Hi All,

I am using edgeR to analyze data from my RNA-Seq time course (4 time points) experiment.

This is my experimental design. It is a factorial design 2*2*4*2*2 with a total of 64 samples (16 samples with 2 replicates each so a total of 32 samples from year1 and 16 samples with 2 replicates each so a total of 32 samples from year 2) comprising of four factors. They are Genotypes with 2 levels, Tissues with 2 levels, Time points with 4 levels and Year with 2 levels. I have the following code to create a design matrix and fit negative binomial GLM and conduct likelihood ratio test.

The following is my code using EdgeR package:

library(edgeR)

#The following is the count data
year<-read.csv("C:/berryallsamplestotalexonscountyear.csv",row.names="Symbol")

group<-factor(c("FP1","FP1","FP2","FP2","FP3","FP3","FP4","FP4","FS1","FS1","FS2","FS2","FS3","FS3","FS4","FS4","MP1","MP1","MP2","MP2","MP3","MP3","MP4","MP4","MS1","MS1","MS2","MS2","MS3","MS3","MS4","MS4","FPA","FPA","FPB","FPB","FPC","FPC","FPE","FPE","FSA","FSA","FSB","FSB","FSC","FSC","FSE","FSE","MPA","MPA","MPC","MPC","MPD","MPD","MPE","MPE","MSA","MSA","MSC","MSC","MSD","MSD","MSE","MSE"))
# Note: labels FP1 until MS4 denotes year 1 samples,FPA until MSE denotes year 2 samples

y<-DGEList(counts=year,group=group)
y
target<-read.csv("C:/targetyear.csv",header=TRUE)
target
Genotype<-factor(target$Genotype)
Tissue<-factor(target$Tissue)
Timepoint<-factor(target$Timepoint)
Year<-factor(target$Year)

design <- model.matrix(~Genotype+Tissue+Timepoint+Year+Genotype:Tissue+Genotype:Timepoint+Genotype:Year+Tissue:Timepoint+Tissue:Year+Timepoint:Year+Genotype:Tissue:Timepoint+Tissue:Timepoint:Year+Genotype:Timepoint:Year+Genotype:Tissue:Timepoint:Year,data=target)
design

y<-estimateGLMCommonDisp(y,design)
y<-estimateGLMTrendedDisp(y,design)
y<-estimateGLMTagwiseDisp(y,design)
y

fit <- glmFit(y,design)
colnames(fit)

Testing for main effects:

topTags(glmLRT(fit,coef=2:5))

Testing for interaction effects:

topTags(glmLRT(fit,coef=6:16))

I have the following questions:

1. I would like to know if the design matrix is correct to test for significance of main and interaction effects of the four variables (Genotype, Tissue, Timepoints, and Year) in the same model.
2. Also, how to interpret the genewise coefficients of all these variables and their interactions after fitting the model if they are significant or not?
3. Is it possible to do the post hoc test like Tukey’s HSD after edgeR? If not, could you please suggest an alternate way?

Thanks for all your time and the suggestions.

Swami

For your first question; you're setting off on a long and difficult road with your current analysis strategy. The concept of main effects only makes sense when the interaction term is negligible. To illustrate, consider a gene that has a positive log-fold change between genotypes for one tissue (I'll ignore all the other factors, to keep it simple) and a negative log-fold change for another tissue. The main effect of genotype has no meaning for this particular gene, because the genotype effect is conditional on the tissue that you're looking at, i.e., there is a significant interaction term between the genotype and tissue factors.

Simply dropping the coefficients corresponding to each factor is not sufficient to test for a "main effect", due to the presence of interaction terms in your design matrix. Instead, you have to identify genes with non-significant interaction terms involving the genotype; prune the design matrix such that those interaction terms are removed; and then repeat the testing for the main effect of genotype. Needless to say, this is a serious pain in the... neck.

For your analysis, I would recommend a one-way layout. It sounds simple, and it is, and that's why it's so nice:

grouping <- factor(paste(Genotype, Tissue, Timepoint, Year, sep="."))
design <- model.matrix(~0 + grouping)
colnames(design) <- levels(grouping)

Each combination of factors is considered to be a group, and the fitted value of each coefficient in the design matrix represents the average log-expression level for the corresponding group. It is trivial to compare between groups with makeContrasts. If you want something akin to a "main effect", you can test whether average of all groups with one genotype is equal to the average of all groups with another genotype. For example, if you have genotypes A and B in tissues X, Y and Z, you could do something like:

con <- makeContrasts((A.X + A.Y + A.Z)/3 - (B.X + B.Y + B.Z)/3, levels=design)

This can be fed to the contrast argument in glmLRT. Alternatively, you could do something like:

con <- makeContrasts(A.X - B.X, A.Y - B.Y, A.Z - B.Z, levels=design)

You could then select genes that have the same sign of log-fold change across all three comparisons , to identify DE genes that are changing between genotypes in a consistent manner across tissues. I feel that this comparison is more useful than testing for a main effect, as it demonstrates that there is significant DE in a consistent direction across all tissues.

As for your second question; the interpretation of the coefficients in your current model is horrendously complicated. This is why I suggest using the one-way layout instead, where the interpretation of each coefficient is a lot easier. Of course, in the one-way layout, it makes no sense to test each coefficient in isolation for significance - that would be like asking whether a gene is significantly expressed, which doesn't really mean much. Rather, you need to use makeContrasts to test for significant differences between coefficients.

Finally; if I recall correctly, the null hypothesis for Tukey's HSD test states that the means for all groups in an ANOVA are equal. If that's all you want to test, you can use makeContrasts to do that for you. Re-using the example above:

con <- makeContrasts(A.X - A.Y, A.X - A.Z, A.X - B.X, A.X - B.Y, A.X - B.Z, levels=design)

This will test for any difference between any of the groups. There's no need to specify, e.g., A.Z - A.Y, because if A.X is equal to A.Y, and A.X is equal to A.Z, then A.Y and A.Z must implicitly be equal (and so on, for all other pairs of groups). I've arbitrarily chosen A.X as a reference in this contrast, but you can use any group you want.