Hello everyone,

I am working on RNA-Seq data. I have 2 conditions (C = control, T = treatment), 4 time points (1, 5, 8 and 10) and 2 or 3 replicates (depending on the quality of the libraries) for each {condition, time point}

I am interested in the genes whose the expression patterns are different in the control and the treatment over time. It means that I want the genes whose the Fold change is different from 1 in at least one time point.

To this end, I build the model :

~ condi + time + time:condi

where condi is my condition (C or T) and time is the time point (1, 5, 8 and 10) so that my design matrix is :

Intercept  condiT  time5  time8  time10  condiT.time5  condiT.time8  condiT.time10
1          0       0      0      0       0             0             0
1          0       0      0      0       0             0             0
1          0       0      0      0       0             0             0
1          1       0      0      0       0             0             0
1          1       0      0      0       0             0             0
1          1       0      0      0       0             0             0
1          0       1      0      0       0             0             0
1          0       1      0      0       0             0             0
1          0       1      0      0       0             0             0
1          1       1      0      0       1             0             0
1          1       1      0      0       1             0             0
1          1       1      0      0       1             0             0
1          0       0      1      0       0             0             0
1          0       0      1      0       0             0             0
1          0       0      1      0       0             0             0
1          1       0      1      0       0             1             0
1          1       0      1      0       0             1             0
1          0       0      0      1       0             0             0
1          0       0      0      1       0             0             0
1          0       0      0      1       0             0             0
1          1       0      0      1       0             0             1
1          1       0      0      1       0             0             1
1          1       0      0      1       0             0             1

Note that at time point {1}, my control (C) and my treatments (T) are biologically identical (they are actual replicates). It means that the differences between C and T at a given time point (5,8 or 10) is modeled by the the interaction term condiT.time(5, 8, or 10), as the condiT term should be 0 for almost every gene. So, I don't need to build a contrast matrix.

For sanity checks, I look at the importance of the condiT term, It gives me only one gene for which this condiT term is significantly different from 0 (pval<0.05), and I am super happy with this. code : topTable(efit, coef=c(2), n=Inf, p.value=0.05) # the coeff '2' is the condiT term

If I check the importance of the condi.Time terms (which is what I am interested in), It gives me 16 genes. code : topTable(efit, coef=c(6:8), n=Inf, p.value=0.05) # the coeff '(6:8)' are the interaction terms

As this is a very few genes, I tried to see what was going on, so I added the condiT tem to the test :

code : topTable(efit, coef=c(2,6:8), n=Inf, p.value=0.05) # same line of code as the one above but with '(2,6:8)' instead of '(6:8)' And this gives me 1531 genes for which this four terms are important (pval<0.05) to model the gene expression...

I just don't get how 16 genes become 1531 genes by adding a term that is not significant taken independently.

Would someone be able to help me ? I am very confused and I have no idea what to trust. Thanks a lot in advance.

Here is my code, if you want to have a look:

design <- model.matrix(  ~  condi + time + time:condi , data=coldata)
colnames(design) <- str_replace_all(colnames(design),":",".") %>% str_replace_all(.,"[()]","") # remove parenthesis from Intercept and change ':' to '.'
idx <- filterByExpr(gse, design) # gse is my expression matrix
dds <- gse[idx,,keep.lib.size=FALSE]
dds <- calcNormFactors(dds)
v <- voom(dds, design)
fit <- lmFit(v, design)
efit <- eBayes(fit)

Best,

Arnaud Stigliani

University of Copenhagen

When you do

topTable(efit, 2, Inf, p.value = 0.05)

You are testing for any gene that is different between the two conditions.

when you do

topTable(efit, 6:8, Inf, p.value = 0.05)

You are now doing an F-test, where you are testing for any differences between those three coefficients, each of which is computing the interaction between T and C, between a given time and time zero. This makes some sense; you want to know if any genes are affected differently at one of the time points. I probably wouldn't do that, but instead would do the three t-tests for each interaction, but it's not my analysis.

If you then add in coefficient 2 to your F-test, you are testing for any gene that has a coefficient that is different between the three interaction terms and the treatment main effect. This doesn't make sense! Comparing four things, one of which is a completely different thing, isn't a reasonable thing to do. So yeah, you get more genes, but you should expect that because you are making a nonsensical comparison.