Hello,

I am trying to use contrasts with a cell means model in EdgeR to compute log fold changes. I am a little confused by the results I get by modifying weights of the contrasts. For instance, I modified the following piece of code from one the EdgeR examples.

#------------ Defining the data and design
> ngenes <- 100
> n1 <- 3
> n2 <- 3
> n3 <-3
> nlibs <- n1+n2+n3
> mu <- 100
> phi <- 0.1
> group <- c(rep(1,n1), rep(2,n2), rep(3,n3))
> design <- model.matrix(~0+as.factor(group))

> ### 4-fold change for the first 5 genes
> i <- 1:5
> fc <- 4
> mu <- matrix(mu, ngenes, nlibs)
> mu[i, 1:n1] <- mu[i, 1:n1]*fc


> counts <- matrix(rnbinom(ngenes*nlibs, mu=mu, size=1/phi), ngenes, nlibs)
> d <- DGEList(counts=counts,lib.size=rep(1e6, nlibs), group=group)

> gfit <- glmFit(d, design, dispersion=phi)

#-------------------------------------- Contrasts --------------

> tr <- glmTreat(gfit, contrast=c(1,1,-1), lfc=1)
> tr2 <- glmTreat(gfit, contrast=c(.5,.5,-1), lfc=1)

> topTags(tr)
Coefficient:  1*as.factor(group)1 1*as.factor(group)2 -1*as.factor(group)3 
       logFC unshrunk.logFC   logCPM       PValue          FDR
7  -14.49698      -14.50068 6.428122 1.029359e-92 5.793101e-91
32 -14.47646      -14.47995 6.482914 1.158620e-92 5.793101e-91
74 -14.26280      -14.26602 6.621232 6.749246e-91 2.249749e-89
72 -14.11476      -14.11748 6.832884 6.579315e-90 1.537966e-88
79 -14.12208      -14.12489 6.711183 7.689829e-90 1.537966e-88
86 -14.07918      -14.08192 6.735456 1.716826e-89 2.861376e-88
11 -14.09690      -14.09979 6.547746 2.035350e-89 2.907644e-88
13 -13.87883      -13.88147 6.573221 1.673322e-87 1.921928e-86
76 -13.84859      -13.85116 6.785719 1.858638e-87 1.921928e-86
50 -13.91007      -13.91282 6.333895 1.921928e-87 1.921928e-86

> topTags(tr2)
Coefficient:  0.5*as.factor(group)1 0.5*as.factor(group)2 -1*as.factor(group)3 
        logFC unshrunk.logFC   logCPM     PValue      FDR
3   1.2569501      1.2577449 7.994976 0.08410010 0.999649
5   1.2546815      1.2557728 7.541736 0.08618485 0.999649
4   1.1923287      1.1931614 7.749677 0.11551910 0.999649
1   1.1707078      1.1716912 7.535702 0.12723752 0.999649
2   0.9689215      0.9695723 7.570693 0.27704382 0.999649
60  0.7809965      0.7821674 6.599907 0.46447814 0.999649
47  0.7787079      0.7797861 6.698892 0.46783455 0.999649
32 -0.7165637     -0.7175457 6.482914 0.53317722 0.999649
72 -0.7129839     -0.7137471 6.832884 0.54020766 0.999649
68  0.6988273      0.6998642 6.588292 0.55360906 0.999649

> counts[72,]
[1] 133  60  94 104  76  87 134 165 155
> gfit$coefficients[72,]
as.factor(group)1 as.factor(group)2 as.factor(group)3 
        -9.253335         -9.325471         -8.795200

What confuses me here is that I expected the log fold change for the first contrast to be less negative than the second, since the weights for the first two factors are reduced in the second. However, since the coefficients are negative, it appears that the results are the opposite. In particular, for the first contrast, I expected the log fold change to be positive or closer to 0, given the counts, say for row 72 (realizing the contrasts do not directly use counts and incorporates more information such as library sizes and dispersions), I intuitively assumed (factor 1+ factor 2) would be more upregulated than factor 3.

Am I making a misinterpreting the outputs? How do I interpret the results?

Thank you!

For the group-means design matrix that you have defined, the contrasts must consist of weights that add to 0, otherwise they don't correspond to comparisons between groups.

In your code, tr2 is a correct contrast that compares group3 to the average of group1 and group2.

Your first contrast however doesn't add to zero and doesn't correspond to any comparison between the groups. It is actually returning a type of intercept rather than a logFC.