creating contrast matrix (limma) for two factorial in R

Question

I am attempting to construct a contrast matrix that I can run in R, using the limma bioconductor package, but I am not sure that I have coded the contrast matrix correctly. A previous post and the limma guide were helpful, but my two factorial design is more complicated than what is illustrated there.

The first factor is the treatment, with two levels (control=c and stress=s), and the second factor is the genotype, with five levels (g1, g2, g3, g4, g5). Each genotype/treatment consists of 3-biological replicates (30xsamples total). My dataset has already been normalized and log2 transformed. It consists of 1208 proteins (based upon spectral counting for those that care) that measures protein abundance differences in the five genotypes and two treatments. The dataset is complete, meaning each sample/condition has a datapoint.

Subset of the data:

proteinID   g1.s1   g1.s2   g1.s3   g1.c1   g1.c2   g1.c3   g2.s1   g2.s2   g2.s3   g2.c1   g2.c2   g2.c3   g3.s1   g3.s2   g3.s3   g3.c1   g3.c2   g3.c3   g4.s1   g4.s2   g4.s3   g4.c1   g4.c2   g4.c3   g5.s1   g5.s2   g5.s3   g5.c1   g5.c2   g5.c3
prot1   -9.70583694 -9.940059478    -9.764489183    -9.691937821    -9.547306096    -9.668928704    -9.821333234    -10.00376839    -9.843380585    -10.0789111 -9.958506961    -9.791583706    -10.04996359    -10.10279896    -10.0689715 -9.989303332    -10.05414639    -10.00619809    -9.907032795    -10.09700113    -10.00902876    -10.05603575    -10.26218387    -10.15527373    -9.88009858 -9.748974338    -9.730010667    -9.899956956    -9.773955101    -9.957684691
prot2   -9.810354967    -9.844319231    -9.896748977    -9.777040294    -9.821308434    -9.906798728    -9.832236541    -9.876359355    -9.935535795    -10.05991278    -9.831098077    -9.789738587    -10.08470861    -10.18515166    -10.10371621    -10.01971224    -9.977142493    -10.09055782    -9.739831978    -9.586647999    -9.949407778    -9.800183583    -9.83900565 -9.943521592    -9.99229056 -9.744850134    -9.794814509    -9.98542989 -9.766324886    -9.95430439
prot3   -11.70842601    -11.72521838    -11.90389475    -11.98273998    -11.915401  -11.88620205    -11.91603643    -11.96029519    -12.14926486    -12.23846499    -12.26650985    -11.84300821    -12.64562082    -12.41471031    -12.66462278    -12.577619  -12.90001898    -12.31577711    -11.66323243    -11.50283992    -11.4844068 -11.60402491    -11.95270942    -11.68245512    -12.32380181    -12.24294758    -12.23990879    -12.21563403    -12.33730369    -12.437377
prot4   -10.88942769    -11.16906693    -11.13942576    -11.31332257    -11.04718433    -11.11811122    -11.17687812    -11.12503828    -10.9724186 -11.16837945    -11.19642214    -10.96468249    -11.3975887 -11.28808753    -11.32778647    -11.34124725    -11.30972182    -11.29564372    -10.74370929    -10.92223539    -10.97733154    -11.40528844    -11.1238659 -11.15938598    -11.24937805    -10.8691392 -11.12478375    -10.75566728    -10.99485703    -11.09493115
prot5   -10.0102959 -9.936796529    -9.964629149    -9.842835973    -9.791578592    -9.773380518    -9.72290866 -9.715837804    -9.79028651 -9.951486129    -9.636225505    -9.820715987    -10.41899204    -10.25269382    -10.26949484    -10.02644184    -10.13120897    -10.20756299    -9.752087376    -9.687001368    -10.07111473    -9.815279198    -9.995624174    -9.993526894    -9.722360141    -9.551502595    -9.551929198    -9.724500546    -9.502769792    -9.65324573
prot6   -10.34051005    -10.27571947    -10.14968761    -10.17419023    -10.47812301    -10.11019796    -10.40447672    -10.15885481    -10.22900798    -10.26612428    -10.21920493    -10.17186677    -10.66125689    -10.95438025    -10.63751536    -10.65825783    -10.60857688    -10.78516027    -10.33890785    -10.49726978    -10.47100414    -10.64742463    -10.78932619    -10.5318634 -10.26494688    -9.975182247    -10.24870036    -10.2356165 -10.26689552    -10.13061368
prot7   -10.24930429    -10.37307132    -10.03573128    -10.29985129    -9.991216794    -10.05854902    -10.1958704 -10.30549818    -10.2078462 -10.28795766    -10.23314344    -10.23897922    -9.997472306    -10.27461285    -10.20805608    -10.06261332    -10.24876706    -10.12643737    -9.906088449    -10.07316322    -10.23545822    -10.30970717    -10.40745591    -10.36432166    -10.22423532    -10.25703553    -10.44925268    -9.902554721    -9.891163766    -10.0695915
prot8   -10.98782595    -10.84184533    -10.76496107    -10.68290092    -10.55763113    -10.91736394    -10.87505278    -10.76474268    -10.58319007    -10.87547281    -10.71948079    -10.95011831    -10.99753277    -11.061728  -10.8852958 -10.86371208    -10.96638746    -11.24112703    -10.46809937    -10.78446288    -10.71240489    -10.80931259    -10.6598091 -10.54801115    -10.70612733    -10.7339808 -10.8184854 -10.53370359    -10.47323989    -10.62675183
prot9   -8.83857166 -8.736344638    -8.743339515    -8.8152675  -8.743086044    -8.719612156    -8.898093257    -8.902781886    -9.071574958    -8.945970659    -8.862394746    -8.825061244    -8.82313363 -9.161452294    -8.905846232    -8.940119002    -9.024995852    -8.943721201    -8.768488159    -8.802155458    -8.721187011    -8.84850416 -8.931513624    -8.86743278 -8.856904592    -8.675257846    -8.900833162    -8.676117406    -8.758661701    -8.925717389
prot10  -10.65297508    -10.74532307    -10.65940071    -10.36671791    -10.50431649    -10.54915637    -11.07154003    -10.79884265    -10.97164196    -11.1201714 -11.14821342    -10.9254445 -10.92875918    -10.90806369    -10.77581175    -11.2324716 -11.31360896    -11.01070959    -11.04450945    -10.89694291    -10.76865867    -10.92983387    -11.07365287    -11.43888216    -11.14948441    -10.69611194    -10.85827316    -10.64470128    -10.79046792    -10.86048168

Code that I am attempting to utilize:

proteins.mat <- as.matrix(proteins.df)
treat = c("g1.s","g1.c","g2.s","g2.c","g3.s","g3.c","g4.s","g4.c","g5.s","g5.c")
factors = gl(10,3,labels=treat)
design <- model.matrix(~0+factors)
colnames(design) <- treat

Here is the design for my model:

> design
   g1.s g1.c g2.s g2.c g3.s g3.c g4.s g4.c g5.s g5.c
1     1    0    0    0    0    0    0    0    0    0
2     1    0    0    0    0    0    0    0    0    0
3     1    0    0    0    0    0    0    0    0    0
4     0    1    0    0    0    0    0    0    0    0
5     0    1    0    0    0    0    0    0    0    0
6     0    1    0    0    0    0    0    0    0    0
7     0    0    1    0    0    0    0    0    0    0
8     0    0    1    0    0    0    0    0    0    0
9     0    0    1    0    0    0    0    0    0    0
10    0    0    0    1    0    0    0    0    0    0
11    0    0    0    1    0    0    0    0    0    0
12    0    0    0    1    0    0    0    0    0    0
13    0    0    0    0    1    0    0    0    0    0
14    0    0    0    0    1    0    0    0    0    0
15    0    0    0    0    1    0    0    0    0    0
16    0    0    0    0    0    1    0    0    0    0
17    0    0    0    0    0    1    0    0    0    0
18    0    0    0    0    0    1    0    0    0    0
19    0    0    0    0    0    0    1    0    0    0
20    0    0    0    0    0    0    1    0    0    0
21    0    0    0    0    0    0    1    0    0    0
22    0    0    0    0    0    0    0    1    0    0
23    0    0    0    0    0    0    0    1    0    0
24    0    0    0    0    0    0    0    1    0    0
25    0    0    0    0    0    0    0    0    1    0
26    0    0    0    0    0    0    0    0    1    0
27    0    0    0    0    0    0    0    0    1    0
28    0    0    0    0    0    0    0    0    0    1
29    0    0    0    0    0    0    0    0    0    1
30    0    0    0    0    0    0    0    0    0    1
attr(,"assign")
[1] 1 1 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$factors
[1] "contr.treatment"

My contrast model. I want to test for interaction, differences between genotypes, and to see if specific genotypes respond differently to the treatment from one another:

cmtx <- makeContrasts(
  GenotypevsTreatment=(g1.s-g1.c)-(g2.s-g2.c)-(g3.s-g3.c)-(g4.s-g4.c)-(g5.s-g5.c),
  genotype=(g1.s+g1.c)-(g2.s+g2.c)-(g3.s+g3.c)-(g4.s+g4.c)-(g5.s+g5.c),
  Treatment=(g1.s+g2.s+g3.s+g4.s+g5.s)-(g1.c+g2.c+g3.c+g4.c+g5.c),
  levels=design)

What my contrast model looks like, but I don't think this is correct:

> cmtx
      Contrasts
Levels GenotypevsTreatment Genotype Treatment
  g1.s                   1        1         1
  g1.c                  -1        1        -1
  g2.s                  -1       -1         1
  g2.c                   1       -1        -1
  g3.s                  -1       -1         1
  g3.c                   1       -1        -1
  g4.s                  -1       -1         1
  g4.c                   1       -1        -1
  g5.s                  -1       -1         1
  g5.c                   1       -1        -1

Fitting the linear model by empirical bayes statistics for differential expression:

fit <- eBayes(contrasts.fit(lmFit(proteins.mat, design), cmtx))
topTable(fit, adjust.method="BH")

The below topTable proteins are the same as the subset of data from above:

> topTable(fit, adjust.method="BH")
       GenotypevsTreatment Genotype    Treatment    AveExpr        F      P.Value    adj.P.Val
prot1        -0.40786338 60.30918  0.073054723  -9.918822 17308.55 1.124646e-39 1.232079e-36
prot2        -0.09255219 59.60864  0.061701713  -9.897968 15801.43 3.304533e-39 1.232079e-36
prot3        -0.23880357 73.48557  0.536672827 -12.090016 15650.65 3.701463e-39 1.232079e-36
prot4        -0.11834000 66.76931  0.305471823 -11.122034 15522.46 4.079731e-39 1.232079e-36
prot5        -0.15210172 59.21509 -0.183849274  -9.876144 14734.51 7.556112e-39 1.423908e-36
prot6        -0.15761118 62.87467  0.155340561 -10.389362 14565.87 8.658504e-39 1.423908e-36
prot7        -0.03886438 61.15652 -0.166795475 -10.182834 14551.88 8.757515e-39 1.423908e-36
prot8        -0.10425341 64.63523 -0.186904167 -10.780359 14461.18 9.429854e-39 1.423908e-36
prot9        -0.03426380 53.48057  0.007403722  -8.854471 13713.49 1.767090e-38 2.021378e-36
prot10       -0.75250251 66.62646  0.327497120 -10.894506 13480.51 2.164184e-38 2.021378e-36

Aside from thinking that I didn’t do this correctly, the result for Genotype looks incorrect to me. Any input would be much appreciated.

Answer 1

The contrast are mean to select what do you want to compare. Your contrast matrix doesn't mean what its names make you think:

cmtx <- makeContrasts(
  GenotypevsTreatment=(g1.s-g1.c)-(g2.s-g2.c)-(g3.s-g3.c)-(g4.s-g4.c)-(g5.s-g5.c),
  genotype=(g1.s+g1.c)-(g2.s+g2.c)-(g3.s+g3.c)-(g4.s+g4.c)-(g5.s+g5.c),
  Treatment=(g1.s+g2.s+g3.s+g4.s+g5.s)-(g1.c+g2.c+g3.c+g4.c+g5.c),
  levels=design)

The GenotypeVsTreatment compares the expression on genome 1 treatment s (g1.s) the genome 2 treatment c, genome 3 treatment c and genome 4 treatment c against everything else. The genotype compares the genome 1 expression against all the other samples. The Treatment compares what it means, the treatment s minus treatment c.

cmtx <- makeContrasts(
      g1vsg2 = (g1.s+g1.c)/2 - (g2.s+g2.c)/2
      g2vsg3 = (g2.s+g2.c)/2 - (g3.s+g3.c)/2
      g3vsg4 = (g3.s+g3.c)/2 - (g4.s+g4.c)/2
      g4vsg5 = (g4.s+g4.c)/2 - (g5.s+g5.c)/2,
      Treatment=(g1.s+g2.s+g3.s+g4.s+g5.s)-(g1.c+g2.c+g3.c+g4.c+g5.c),
      g1 = g1.s-g1.c
      g2 = g2.s-g5.c
      g3 = g3.s-g5.c
      g4 = g4.s-g5.c
      g5 = g5.s-g5.c
      levels=design)

The givsgj are the different comparisons of one genome against the other, to see differences between genomes. With your design you can't test it in general if they are different. The gi compares the effect of the treatment in a genome. If the result of g1, g2, g3, g4, and g5 is different then the genomes differ on the response.