3.3 Multi-factor comparison
Some RNA-seq experiments include multiple experimental factors, such as treatment, time point, tissue, genotype, and batch. A multi-factor model can account for these factors simultaneously and can also test whether the effect of one factor depends on another.
Common biological questions include:
- Is there a difference between treatment and control?
- Is there a difference between tissues?
- Does the treatment effect depend on time point?
- Does the genotype effect depend on tissue?
- Within a treatment group, does gene expression differ between time points?
The statistical design should be selected according to the biological question.
In this example, we analyze soybean RNA-seq data from control and treated leaves collected at several time points after water-deficit treatment (Rodrigues et al. 2015). The dataset is available from the EBI Expression Atlas under accession E-GEOD-69469.
atlas <- getAtlasData("E-GEOD-69469")
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)## SRR2048167 SRR2048168 SRR2048169 SRR2048170 SRR2048171 SRR2048172 SRR2048173 SRR2048174 SRR2048175 SRR2048176 SRR2048177 SRR2048178 SRR2048179 SRR2048180 SRR2048181 SRR2048182 SRR2048183 SRR2048184 SRR2048185 SRR2048186 SRR2048187 SRR2048188 SRR2048189 SRR2048190 SRR2048191 SRR2048192 SRR2048193 SRR2048194 SRR2048195 SRR2048196 SRR2048197 SRR2048198 SRR2048199 SRR2048200 SRR2048201 SRR2048202
## ENSRNA049760259 46 43 16 50 68 29 11 28 28 40 28 9 30 70 91 11 28 73 37 97 61 21 65 23 58 44 180 34 27 67 78 34 33 54 16 41
## ENSRNA050000373 130 104 104 84 52 48 35 88 111 87 59 25 93 110 207 127 27 72 73 124 204 49 178 47 177 114 387 96 40 144 200 66 60 82 38 109
## ENSRNA050000380 3372 2412 2216 3014 2822 1702 1151 2824 2202 1907 1641 983 2501 2704 6699 1265 1023 3336 1915 2989 3910 990 3840 1145 2698 2381 5155 1587 971 2131 3103 1267 1090 1379 632 1918
## ENSRNA050000423 2971 2103 1634 2091 1876 1647 749 1796 1602 1283 1220 447 662 949 3485 318 562 1336 875 2153 4403 414 1419 695 1574 1268 4420 754 692 1684 2685 909 929 1225 556 1231
## ENSRNA050000426 5143 4175 3686 5639 5197 2034 3378 6924 6188 5320 2876 2233 5349 5616 14272 2169 2437 6220 3927 7912 8149 1735 5653 2003 5786 3694 9436 2461 1668 3805 5656 2270 1932 2545 1140 3016
## ENSRNA050001862 163 133 125 185 185 69 110 198 234 170 93 80 141 173 417 80 80 180 125 231 214 65 177 55 177 135 295 75 49 125 191 52 56 79 53 94col_data <- as.data.frame(colData(se))
exp_design <- data.frame(
group = factor(ifelse(grepl('control', col_data[, c("environmental_stress")]), 'control', 'stress')),
timepoint = factor(gsub("Zeitgeber time ", "", col_data[, c("sampling_time_point")])))
exp_design## group timepoint
## 1 stress 0
## 2 stress 0
## 3 stress 0
## 4 control 0
## 5 control 0
## 6 control 0
## 7 stress 4
## 8 stress 4
## 9 stress 4
## 10 control 4
## 11 control 4
## 12 control 4
## 13 stress 8
## 14 stress 8
## 15 stress 8
## 16 control 8
## 17 control 8
## 18 control 8
## 19 stress 12
## 20 stress 12
## 21 stress 12
## 22 control 12
## 23 control 12
## 24 control 12
## 25 stress 16
## 26 stress 16
## 27 stress 16
## 28 control 16
## 29 control 16
## 30 control 16
## 31 stress 20
## 32 stress 20
## 33 stress 20
## 34 control 20
## 35 control 20
## 36 control 20Create a SeqCountData object from the prepared count matrix and experimental design.
During DEGES-based normalization, candidate DEGs are removed from the data
before the normalization factors are recomputed from the remaining genes.
By default, calc_nf() screens candidate DEGs using the available experimental design,
so genes associated with any modeled effect can be excluded from normalization factor estimation.
If normalization should remove candidate DEGs for a specific comparison, such as control versus treatment after adjusting for time point, see the Advanced Usage chapter.
Then convert the normalized object into a edgeR DGEList
and run a edgeR quasi-likelihood test for the treatment effect while adjusting for time point;
that is, coef = 2 in this design.
d <- as_DGEList(x)
design <- model.matrix(~ group + timepoint, data = exp_design)
d <- estimateDisp(d, design = design)
fit <- glmQLFit(d, design = design)
qlf <- glmQLFTest(fit, coef = 2)Do not assume that a biological effect is always represented by coef = 2.
The coefficient number depends on the factor levels, their reference levels,
the order of terms in the model, and the presence of interaction terms.
Before running a statistical test, inspect the design matrix column names.
After confirming the meaning of each coefficient, test the coefficient of interest.
In this example, coef = 2 corresponds to the treatment effect used above.
## [1] "(Intercept)" "groupstress" "timepoint12" "timepoint16" "timepoint20" "timepoint4" "timepoint8"Extract the results.
## logFC logCPM F PValue FDR
## ENSRNA049760259 -0.3372351845 2.694545 6.639535e+00 0.01425689 0.0620473
## ENSRNA050000373 -0.0958208298 3.861405 3.310878e-01 0.56862490 0.7237248
## ENSRNA050000380 -0.1593806127 8.333032 4.934011e+00 0.03275791 0.1109259
## ENSRNA050000423 -0.0003588716 7.738853 9.163468e-06 0.99760154 0.9988933
## ENSRNA050000426 0.0110208839 9.273367 3.896653e-02 0.84463327 0.9109317
## ENSRNA050001862 -0.0942038960 4.316354 1.342311e+00 0.25431019 0.4369289Here, the analysis compares treatment groups while accounting for differences among time points. When the model includes interaction terms, the coefficient for a main effect describes that effect only at the reference level of the other factor. It should not be interpreted as the average effect across all tissues, time points, or developmental stages. Therefore, the coefficient or contrast must be selected according to the specific biological comparison of interest.
The same multi-factor comparison can be performed with DESeq2.
library(DESeq2)
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 + timepoint)
dds <- estimateDispersions(dds)
dds <- nbinomLRT(dds, full = ~ group + timepoint, reduced = ~ timepoint)
test_stats <- results(dds)
head(test_stats)## log2 fold change (MLE): timepoint 8 vs 0
## LRT p-value: '~ group + timepoint' vs '~ timepoint'
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSRNA049760259 47.8538 -0.419818 0.2282933 5.96967911 0.0145539 0.0598528
## ENSRNA050000373 110.1385 -0.202770 0.2897823 0.27864507 0.5975907 0.7477595
## ENSRNA050000380 2456.5260 -0.600964 0.1196823 4.47284045 0.0344377 0.1107921
## ENSRNA050000423 1625.7716 -1.671562 0.2133009 0.00811839 0.9282061 0.9609125
## ENSRNA050000426 4715.0534 -0.267795 0.0818337 0.24745743 0.6188711 0.7634494
## ENSRNA050001862 152.0347 -0.418654 0.1291651 0.75774444 0.3840353 0.5680314The interpretation of an LRT is determined by the difference between the full and reduced models.
Here, the reduced model removes all terms containing group.
The test therefore asks whether treatment-related terms improve model fit.
It does not represent one simple pairwise contrast.