Entering edit mode
Dear all, I am performing a very simple limma voom analysis on 2 groups, but I am getting different results when I use make.contrasts() and when I don't. The aim is simply to find differentially expressed genes between 2 groups. I have read the limma voom guide and edgeR guide but they seem to give different examples (using make.contrasts and not, which gives me different results). I am wondering if you might be able to identify the issue? Thanks very much in advance.
myvec
[1] 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
[38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[75] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1
[149] 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1
[186] 1 1
d <- DGEList(counts=samples2, group=myvec)
d <- calcNormFactors(d)
group <- d$samples$group
mm <- model.matrix(~0 + group)
y <- voom(d, mm, plot = F)
#### here I try to fit without contrasts
fit <- lmFit(y, mm)
fit2 <- eBayes(temp)
#### here I try to fit with contrasts
contr <- makeContrasts(group1 - group0, levels = colnames(coef(fit)))
temp <- contrasts.fit(fit, contr)
fit3 <- eBayes(temp)
#### topTable(fit2) and topTable(fit3) produce very different results, with p values in the former being 10e-300 and those in the latter being 10e-6
> head(fit3$coefficients)
group0 group1
ENSG00000227232 0.9232934 0.7042146
ENSG00000225972 -0.2555306 -1.3402430
ENSG00000225630 6.6519330 6.4041992
ENSG00000237973 4.1737093 4.0783065
ENSG00000229344 1.7111760 1.3761588
ENSG00000248527 8.3435431 8.6294344
> head(fit2$coefficients)
Contrasts
group1 - group0
ENSG00000227232 -0.21907877
ENSG00000225972 -1.08471239
ENSG00000225630 -0.24773378
ENSG00000237973 -0.09540285
ENSG00000229344 -0.33501722
ENSG00000248527 0.28589126
> sessionInfo( )
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /gpfs3/apps/eb/2020b/skylake/software/FlexiBLAS/3.0.4-GCC-11.2.0/lib64/libflexiblas.so.3.0
locale:
[1] C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] edgeR_3.34.1 limma_3.48.3
loaded via a namespace (and not attached):
[1] MASS_7.3-54 compiler_4.1.2 Rcpp_1.0.7 grid_4.1.2
[5] locfit_1.5-9.4 statmod_1.4.36 lattice_0.20-45
Dear Gordon, Thank you very much indeed. Could I just confirm, if I wanted to compare group0 and group1 adjusting for gender, then this would be the way to do it?
Yes, but please do not use `method="robust", which is not recommended in any of the limma workflow examples. See for example simultaneous use of robust and weighting methods in limma. Unless you are an expert able to make your own decisions it is better to stick to the standard recommended options.