3.2 Multi-group comparison
RNA-seq experiments often contain more than two groups, for example when comparing expression among different treatments or distinct tissues. Here, we load RNA-seq samples from multiple tissues of Oryza sativa Japonica Group E-MTAB-2037 (Sakai et al. 2011) to identify genes whose expression varies among tissues.
atlas <- getAtlasData("E-MTAB-2037")
se <- atlas[[1]]$rnaseq
count_data <- as.matrix(assay(se, "counts"))
keep <- rowSums(cpm(count_data) > 5) > ncol(count_data) / 2
count_data <- count_data[keep, , drop = FALSE]
count_data <- as.matrix(count_data)
head(count_data)## DRR001024 DRR001025 DRR001026 DRR001027 DRR001028 DRR001029 DRR001030 DRR001031 DRR001032 DRR001033 DRR001034 DRR001035 DRR001036 DRR001037 DRR001038 DRR001039 DRR001040 DRR001041 DRR001042 DRR001043 DRR001044 DRR001045 DRR001046 DRR001047 DRR001048 DRR001049 DRR001050 DRR001051
## ENSRNA049464535 9 8 20 23 86 45 31 5 13 21 42 53 38 28 3 2 5 11 31 23 12 4 2 16 12 45 12 10
## ENSRNA049465003 6 10 18 30 39 49 71 5 7 20 21 30 44 34 1 1 2 5 17 7 21 5 1 10 13 28 15 44
## ENSRNA049465125 105 2 130 54 45 3 11 38 3 116 50 32 2 8 24 2 49 13 14 1 3 57 0 128 44 47 1 13
## ENSRNA049468277 27 54 22 52 23 49 51 10 36 25 31 21 25 30 4 9 6 7 11 9 18 0 0 0 0 1 0 4
## ENSRNA049469734 56 24 89 97 170 254 143 20 18 59 81 147 182 136 13 6 16 12 59 76 68 14 4 31 22 95 57 91
## ENSRNA049471102 374 306 330 508 292 341 471 122 221 379 359 232 247 311 112 45 73 101 80 78 177 283 54 313 197 170 79 303col_data <- as.data.frame(colData(se))
exp_design <- data.frame(
group = factor(col_data[, c("organism_part")]))
exp_design## group
## 1 leaf
## 2 root
## 3 shoot
## 4 pre-flowering panicle
## 5 post-flowering panicle
## 6 callus
## 7 seed
## 8 leaf
## 9 root
## 10 shoot
## 11 pre-flowering panicle
## 12 post-flowering panicle
## 13 callus
## 14 seed
## 15 leaf
## 16 root
## 17 shoot
## 18 pre-flowering panicle
## 19 post-flowering panicle
## 20 callus
## 21 seed
## 22 leaf
## 23 root
## 24 shoot
## 25 pre-flowering panicle
## 26 post-flowering panicle
## 27 callus
## 28 seedThe count matrix and sample information are then combined into a SeqCountData object.
By default, calc_nf() uses TMM normalization and a edgeR
quasi-likelihood test during DEG elimination.
In a multi-group design, the default DEGES-based workflow screens candidate DEGs across all groups,
analogous to an omnibus test.
If normalization should remove candidate DEGs for a specific comparison,
such as leaf versus seed, see the Advanced Usage chapter.
Convert the normalized object into a edgeR DGEList,
then perform an omnibus DEG test across tissues.
d <- as_DGEList(x)
design <- model.matrix(~ group, data = exp_design)
d <- estimateDisp(d, design = design)
fit <- glmQLFit(d, design = design)
qlf <- glmQLFTest(fit, coef = 2:ncol(design))
test_stats <- topTags(qlf, n = nrow(d$counts), sort.by = "none")$table
head(test_stats)## logFC.groupleaf logFC.grouppost.flowering.panicle logFC.grouppre.flowering.panicle logFC.grouproot logFC.groupseed logFC.groupshoot logCPM F PValue FDR
## ENSRNA049464535 -1.8106164 -0.1261421 -0.7197434 -1.89328712 -0.9999078 -0.9361496 3.435229 11.548901 2.477801e-07 2.971281e-07
## ENSRNA049465003 -2.0446002 -0.9691048 -1.0131605 -2.21956914 0.1254135 -1.1900287 3.317926 18.947184 4.207533e-10 5.818692e-10
## ENSRNA049465125 5.5948607 3.2434212 4.1755554 0.33022784 1.8109071 5.8879459 4.170530 94.207388 1.070822e-22 4.375462e-22
## ENSRNA049468277 -0.3282526 -1.5110287 -0.1838303 0.60239310 -0.1125906 -0.6267799 3.257092 2.488242 3.958891e-02 4.004107e-02
## ENSRNA049469734 -1.7761631 -1.2270622 -1.7335223 -3.07740234 -0.7886127 -1.5386057 5.087633 47.970518 2.895103e-16 6.148744e-16
## ENSRNA049471102 1.0403364 -0.8880537 0.3994052 0.06660787 0.3651634 0.6209549 6.843600 18.220871 7.563264e-10 1.028933e-09coef = 2:ncol(design) tests whether each gene differs among the tissue groups.
For a specific tissue comparison,
set different values in coef or use the contrast argument.
Refer to the edgeR documentation for detailed usage.
The same multi-group comparison can be performed with DESeq2 using an LRT.
x <- newSeqCountData(count_data, exp_design = exp_design)
x <- calc_nf(x, norm_func = norm_rle, test_func = test_deseq2)
dds <- as_DESeqDataSet(x, design = ~ group)
dds <- estimateDispersions(dds)
dds <- nbinomLRT(dds, full = ~ group, reduced = ~ 1)
test_stats <- results(dds)
head(test_stats)## log2 fold change (MLE): group shoot vs callus
## LRT p-value: '~ group' vs '~ 1'
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSRNA049464535 20.4498 -0.895527 0.303373 56.73906 2.06318e-10 2.48946e-10
## ENSRNA049465003 17.7081 -1.121003 0.300518 83.57475 6.51346e-16 8.67411e-16
## ENSRNA049465125 44.8311 6.096052 0.602161 472.07131 8.70206e-99 2.98878e-98
## ENSRNA049468277 18.4052 -0.385961 1.171369 3.12669 7.92786e-01 7.93205e-01
## ENSRNA049469734 68.3128 -1.486553 0.183893 239.43985 7.39359e-49 1.55523e-48
## ENSRNA049471102 256.2497 0.700061 0.256541 62.75730 1.23688e-11 1.52301e-11