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.
