Hello,

I am  quite new to limma and I am looking for some advice regarding the matrix disign and contrast matrix.

I have RNA-seq data from 3 genotypes and 2 time points (not a time course), 2 replicates each

The targets file looks as follow:

gen    Read1      Read2        out           zt

N        NAR1.fa   NAR2.fa    NA.bam     0

N         NBR1.fa  NBR2.fa     NB.bam    12 

C....

C  ......

D    ......

D  ..........  each line repeated twice

I am interested in DEGs in N-C and D-C and that go in opposite direction in N and D compared to C. So genes that are high in N (compared to C) and low in D or low in N and high in D. I would like to use the zt information to identify gen zt interactions in the subset of genes identified above.

The command I used are:

> gen <- factor(targets$gen, levels=c("N","C","D"))
> design <- model.matrix(~0+gen)
> colnames(design) <- c("N","C","D")

> v <- voom(dge,design,plot=TRUE)

> fit <- lmFit(v,design)
> cont.matrix <- makeContrasts(NC=N-C, DC=D-C, Diff=(N-C)-(D-C) ,levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)

In this design I am non considering zt.

How can I include zt in order to find gen~zt interactions?

Thanks,

Mirko

The concept of interactions is often unnecessary for DE analyses. Instead, it is usually more intuitive to parametrize the design in a one-way layout:

grouping <- factor(paste0(targets$gen, targets$zt))
design <- model.matrix(~0 + grouping)
colnames(design) <- levels(grouping)

You can then specify contrasts between two or more groups, using makeContrasts as above. If you want to mimic the contrasts you already have, you could compare between the average of the two timepoints for each genotype:

con.CvD <- makeContrasts((C0+C12)/2 - (D0 + D12)/2, levels=design)

You'll have to decide what you mean when you want to look at "gen~zt interactions". My interpretation would be something like, "is the magnitude of DE between genotypes affected by time?". If so, you could do something like:

makeContrasts((C0 - D0) - (C12 - D12), levels=design)

This will test for any genes that differ in the size of the C - D log-fold change between timepoints.

Note that this is equivalent to dropping the interaction term in a classical interaction model. Let's describe the effect of genotype in the C-to-D comparison as D0 - C0, and the effect of time as C12 - C0. Under the null hypothesis of no interaction, the combined effect of time and genotype is additive, so we get:

(D0 - C0) + (C12 - C0) = (D12 - C0)

which, upon rearranging, is the same as the contrast above.