I'm using glmLRT() from edgeR to analyze some RNA-seq data. My experimental design has one Treatment factor with three classes (A, B and C), one continuous Age covariate and one Batch factor with 19 levels.
> colnames(design)
[1] "(Intercept)" "ClassB" "ClassC" "Age"
[5] "Batch1" "Batch2" "Batch3" "Batch4"
[9] "Batch5" "Batch6" "Batch7" "Batch8"
[13] "Batch9" "Batch10" "Batch11" "Batch12"
[17] "Batch13" "Batch14" "Batch15" "Batch16"
[21] "Batch17" "Batch18" "Batch19” This is an ANCOVA design, so when I make my comparisons between the three Treatment classes (A, B, C), I have been taught that I need to pick a specific covariate value to make the contrasts. Typically you might pick Age = 0, Age = min(Age), Age = mean(Age or Age = max(Age) for your contrasts. I've been estimating the contrasts for the mean Age. Let's suppose the mean age is 7.3 years for this example. I'm using the code
> lrt1 = glmLRT(fit,contrast=c(0,1,0,mean(Age),rep(0,19)))
... to compute my contrasts. The code runs and the output makes sense, but my collaborators would like some assurance that glmLRT() is computing the correct contrast. Unfortunately, there are no explanations of covariate effects and contrasts in the edgeR manuals or tutorials, and I can't find any information on Google either. Can glmLRT compute contrasts adjusted for a continuous covariate like this? Or am I computing nonsense here?
First, reply to answers using the "Add Comment" button.
Second, the contrast as posted makes no sense. The null hypothesis corresponding to this contrast is:
ClassB + m*Age = 0 # or m*Age = -ClassB... where
mis the mean age,classBis the log-fold change of B over A, andAgeis the slope w.r.t. age.That is, the log-fold change of class A over B is equal to the gradient of expression w.r.t. time scaled by the mean age. This has no obvious scientific meaning to me. If you want to test for DE between classes, just drop the
Classcoefficients. Think ofAgeas a blocking factor; required for the model fit but not involved in any comparisons (at least, not between classes).Your advice to always use Age = 0 in my contrasts implies that I can only make contrasts at the Y-intercept. That's problematic if there is a significant Age*Class interaction, because the fold changes could change from minAge to meanAge to maxAge. Obviously, I should have specified the design as Y ~ Class*Age + Batch in my original email. That was my mistake. Now, I've coded an example to demonstrate what I was thinking with my contrasts. EdgeR uses a neg binomial family generalized linear model with eBayes "shrinking", but I'm going to approximate the edgeR model using Poisson regression from glm().
# Sample Data Y = c(7,8,8,10,12,11,13,15,16,18,15,13,10,9,6) lY = log(Y) Class = factor(c(rep("A", 5),rep("B", 5),rep("C",5)),levels=c("A","B","C")) Age = c(0.6,2.1,3.3,4.2,7.7,2.9,5.4,6.5,8.0,12.8,1.4,4.7,6.8,8.5,13.1) # Basic Plot plot(x=Age,y=lY,pch=c(rep(17,5),rep(16,5),rep(15,5)), col=c(rep("orange",5),rep("purple",5),rep("green",5)), cex=2,ylim=c(1.6,3.0),yaxt="n",xlim=c(0,14),ylab="Count") axis(side=2,at=log(c(5,10,15,20)),labels=c(5,10,15,20)) legend("bottomleft",c("A","B","C"),pch=17:15,col=c("orange","purple","green"), title="Class") # Fit Poisson regression model with glm() poisfit = glm(Y ~ Class*Age,family=poisson) poiscoef = poisfit$coef summary(poisfit) # Compute regression equations for A, B and C A.reg = poisfit$coef[c(1,4)] B.reg = c(sum(poisfit$coef[1:2]),sum(poisfit$coef[4:5])) C.reg = c(sum(poisfit$coef[c(1,3)]),sum(poisfit$coef[c(4,6)])) # Add regression lines to plot abline(A.reg,col="orange") abline(B.reg,col="purple") abline(C.reg,col="green") # Suppose I compute contrasts at Age = 1.8, Age = 5.4 and Age = 12.2 abline(v=c(1.8,5.4,12.2),lty=3,col="red") A0 = sum(poiscoef*c(1,0,0,0,0,0)) B0 = sum(poiscoef*c(1,1,0,0,0,0)) C0 = sum(poiscoef*c(1,0,1,0,0,0)) A1.8 = sum(poiscoef*c(1,0,0,1.8,0,0)) B1.8 = sum(poiscoef*c(1,1,0,1.8,1.8,0)) C1.8 = sum(poiscoef*c(1,0,1,1.8,0,1.8)) A5.4 = sum(poiscoef*c(1,0,0,5.4,0,0)) B5.4 = sum(poiscoef*c(1,1,0,5.4,5.4,0)) C5.4 = sum(poiscoef*c(1,0,1,5.4,0,5.4)) A12.2 = sum(poiscoef*c(1,0,0,12.2,0,0)) B12.2 = sum(poiscoef*c(1,1,0,12.2,12.2,0)) C12.2 = sum(poiscoef*c(1,0,1,12.2,0,12.2)) # Check mean responses against regression lineThis is a simplified example of my edgeR analysis for one gene. ClassA and ClassB are mostly parallel, while ClassC (green) is negative and almost perpendicular to A and B. This is an obvious interaction effect. I fit the Poisson regression and plotted the regression lines. I compute 9 different predicted values for Classes A, B and C at Ages 1.8, 5.4 and 12.2, then I compute the contrasts of interest two ways. Again, if you only compute contrasts using Age = 0, then you only compare the lines at Y-int. That would only tell you part of the story where A < B < C. I specified the design and contrasts incorrectly in my initial post, but I don't think my overall idea was flawed. However, the edgeR model is more complicated than glm(), so I was hoping someone could confirm that the contrasts would also work in this way in edgeR.
For simplicity, I'm going to reparametrize your example in edgeR terms:
design <- model.matrix(~0 + Class + Class:Age) colnames(design) <- sub(":", "_", colnames(design))The renaming just keeps
makeContrastshappy later. Looking at the column names gives us:The first three terms are the intercepts for each group, while the last three terms are the gradients with respect to age. If you want to test for DE between classes at a specific timepoint
X, you would do:con <- makeContrasts((ClassA + X * ClassA_Age) - (ClassB + X * ClassB_Age), levels=design)... which may or may not be interesting, depending on the importance of
X. I presume you wanted to usemean(Age)asXin your original post, but there's no statistical reason that the mean age would be of particular importance.Usually, the contrast of greater interest/relevance is whether there is any DE between classes at any time point:
The null hypothesis here would be that the slope and intercept are equal for both classes.
Thank you. This is what I was looking for. You are right that mean(X) isn't especially important, but it is a helpful way to code "compare the groups somewhere in the middle of the range of X" without specifying an exact X value. Usually the 3 places you would want to look are
1. Somewhere near min X or zero
2. Somewhere near mean X or median X or any where in the middle
3. Somewhere near max X
Looking at the contrasts near min X and max X will allow you to look at the two higher leverage areas of your experiment. If the fold change difference are going to change direction, they are most likely to do it between min X and max X. Looking somewhere in the middle will get you a comparison away from the high leverage areas, in an area where the regression lines are most likely to cross.
What's the null hypothesis here? I don't see the statistical necessity of testing a whole bunch of points in the covariate range - at least, not with respect to the leverages. If the scientific question involves those specific ages, then I could understand; but if you want to make a statement about the differences in slopes/intercepts between classes, it would be better/simpler to do that directly. I usually just check the expression profiles manually for any genes of interest to make sure that they're not being driven by a few influential observations (or, in extreme cases, I repeat the DE analysis after tossing those observations out, to ensure that the results are robust).