**60**wrote:

Dear all,

I am having some conceptual issues with the contrasts generated by limma. From my understanding, contrasts acts as an method for defining comparisons between groups. For example, we have expression data of six samples with two samples being controls ("trt1"), two samples being treated with drug A ("trt2") and the remaining two samples being treated with drug B ("trt3"). Here I would like to compare the following:

- control versus drug A (trt1 vs trt2)
- control versus drug B (trt1 vs trt3)
- drug A versus drug B (trt2 vs trt3) - Not as important, but nonetheless interesting

To this end, I thought of generating contrasts using the makeContrasts() function in limma. However, I thought that defining contrasts was conceptually the same as extracting the expression data of the samples trt1 and trt2, trt1 and trt3 and trt2 and trt3 and generating a linear model using lmFit() for each comparison separately. Could anybody tell me what causes the discrepancy or where I could find documentation that describes the difference?

Sample data:

set.seed(1) #Simulate data Mvals <- matrix(rnorm(n = 36), ncol = 6) Mvals[,3:4] <- Mvals[,3:4] + 5 Mvals[,5:6] <- Mvals[,5:6] + 10 colnames(Mvals) <- c("sam1_trt1", "sam2_trt1", "sam3_trt2", "sam4_trt2", "sam5_trt3", "sam6_trt3") #Factor of interest fac_int <- as.factor(rep(c("trt1", "trt2", "trt3"), each = 2)) #Contrasts using limma design <- model.matrix(~0 + fac_int) colnames(design) <- c("trt1", "trt2", "trt3") lmfit <- lmFit(Mvals, design) cont <- makeContrasts(trt1-trt2, trt1-trt3, trt2-trt3, levels = design) lmfit.cont <- contrasts.fit(lmfit, cont) lmfit.cont.ebayes <- eBayes(lmfit.cont) lmfit.cont.ebayes$p.value #Contrasts by extracting the correct Mvals and design #trt1-trt2 M12 <- Mvals[,which(fac_int %in% c("trt1", "trt2"))] design12 <- model.matrix(~as.factor(rep(c("trt1", "trt2"), each = 2))) lmfit12 <- lmFit(M12, design12) lmfit12.ebayes <- eBayes(lmfit12) lmfit12.ebayes$p.value[,2] #trt1-trt3 M13 <- Mvals[,which(fac_int %in% c("trt1", "trt3"))] design13 <- model.matrix(~as.factor(rep(c("trt1", "trt3"), each = 2))) lmfit13 <- lmFit(M13, design13) lmfit13.ebayes <- eBayes(lmfit13) lmfit13.ebayes$p.value[,2] #trt2-trt3 M23 <- Mvals[,which(fac_int %in% c("trt2", "trt3"))] design23 <- model.matrix(~as.factor(rep(c("trt2", "trt3"), each = 2))) lmfit23 <- lmFit(M23, design23) lmfit23.ebayes <- eBayes(lmfit23) lmfit23.ebayes$p.value[,2] #Compare the contrast with the extraction lmfit.cont.ebayes$p.value # trt1 - trt2 trt1 - trt3 trt2 - trt3 # [1,] 0.0001299617 5.130891e-09 3.057964e-05 # [2,] 0.0025873059 5.503494e-08 4.186192e-05 # [3,] 0.0000182689 1.659910e-08 1.207367e-03 # [4,] 0.0003252686 2.228194e-07 1.930707e-03 # [5,] 0.0011712039 4.689758e-07 1.342449e-03 # [6,] 0.0003469217 1.733174e-08 6.331602e-05 data.frame(trt1_trt2 = lmfit12.ebayes$p.value[,2], trt1_trt3 = lmfit13.ebayes$p.value[,2], trt2_trt3 = lmfit23.ebayes$p.value[,2]) # trt1_trt2 trt1_trt3 trt2_trt3 #1 0.0019729784 3.487858e-08 7.756619e-05 #2 0.0149031837 1.962147e-07 1.040529e-04 #3 0.0005273948 8.156824e-08 1.523928e-03 #4 0.0036619967 5.550043e-07 2.320744e-03 #5 0.0087060279 9.729137e-07 1.674585e-03 #6 0.0038247856 8.417320e-08 1.459670e-04 sessionInfo() # R version 3.2.2 (2015-08-14) # Platform: x86_64-apple-darwin13.4.0 (64-bit) # Running under: OS X 10.11 (El Capitan) # # locale: # [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 # # attached base packages: # [1] stats graphics grDevices utils datasets methods base # # other attached packages: # [1] limma_3.24.15 BayesFactor_0.9.12-2 Matrix_1.2-2 coda_0.17-1 # # loaded via a namespace (and not attached): # [1] magrittr_1.5 tools_3.2.2 pbapply_1.1-1 MatrixModels_0.4-1 Rcpp_0.12.1 # [6] mvtnorm_1.0-3 stringi_0.5-5 grid_3.2.2 stringr_1.0.0 gtools_3.5.0 # [11] lattice_0.20-33