limma voom contrasts
1
0
Entering edit mode
S T ▴ 80
@s-t-4287
Last seen 3 months ago
United Kingdom

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

contrasts limma voom • 230 views
0
Entering edit mode
@gordon-smyth
Last seen 9 hours ago
WEHI, Melbourne, Australia

If you don't wish to define a contrast then you must use

mm <- model.matrix(~group)


without the 0+ term. This is explained in the limma User's Guide Section 9.2.

0
Entering edit mode

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?

mm <- model.matrix(~0 + group + gender )
y <- voom(d, mm, plot = F)
fit <- lmFit(y, mm, method="robust")
contr <- makeContrasts(group1 - group0, levels = colnames(coef(fit)))
temp <- contrasts.fit(fit, contr)
fit2 <- eBayes(temp,robust=T)

0
Entering edit mode

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.