4.3 GLM for multi-factor comparison

For experiments with multiple experimental factors, the statistical model should be chosen to match the biological question. You can supply different models to screen candidate DEGs for specific effects during normalization.

Here, we load and prepare a dataset for a multi-factor comparison.

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)

col_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")])))

x <- newSeqCountData(count_data, exp_design = exp_design)
x

## SeqCountData object
##   genes:   19332
##   samples: 36
##   groups:  control, stress
##   factors: not set

By default, calc_nf() screens candidate DEGs using the available experimental design. This can exclude genes associated with any modeled additive or interaction effects during DEGES-based normalization.

x <- calc_nf(x)
x@meta$nf

## 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 
##  0.8288599  0.9153675  0.7937189  0.9208547  0.8953694  0.9057120  0.8960018  0.8345749  0.8658727  0.8939510  0.9283132  0.7419535  1.0964737  0.9575650  0.9904778  0.6272125  1.0611643  1.0086907  1.0860319  1.2016754  1.1830043  1.1345363  0.9491144  1.1672350  1.0581581  1.1260418  1.2245340  1.1272566  1.1640574  1.2205572  1.0663899  1.0116131  1.0336780  1.0482278  1.0565291  0.9792259

In some analyses, the normalization step should screen candidate DEGs for a specific factor only. In that case, pass additional arguments to the test function used by calc_nf().

For example, using TMM normalization and a edgeR quasi-likelihood test, the workflow below removes candidate DEGs for the treatment effect while adjusting for time point.

design <- model.matrix(~ group + timepoint, data = exp_design)
head(design)

##   (Intercept) groupstress timepoint12 timepoint16 timepoint20 timepoint4 timepoint8
## 1           1           1           0           0           0          0          0
## 2           1           1           0           0           0          0          0
## 3           1           1           0           0           0          0          0
## 4           1           0           0           0           0          0          0
## 5           1           0           0           0           0          0          0
## 6           1           0           0           0           0          0          0

After confirming that the second column of the design matrix represents the treatment effect, pass design and coef = 2 to calc_nf().

x <- calc_nf(x, design = design, coef = 2)
x@meta$nf

## 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 
##  0.8434084  0.9188214  0.8050358  0.9327536  0.9093589  0.9379802  0.8954608  0.8460603  0.8687117  0.9072096  0.9536584  0.7639547  1.0838529  0.9539572  0.9892172  0.6350583  1.0692729  1.0085336  1.0637426  1.1773847  1.1646484  1.0859950  0.9280794  1.1391140  1.0509235  1.1131365  1.2298603  1.1127919  1.1472473  1.2113866  1.0752994  1.0172603  1.0358385  1.0590556  1.0742533  0.9916768

See the edgeR documentation for more detailed guidance on model matrices, coefficients, and contrasts.

DESeq2 uses a different testing interface. For an LRT, specify a full model and a reduced model. The test evaluates whether terms present only in the full model explain additional variation. For example, to test the treatment effect after adjusting for time point, use the following workflow.

x <- newSeqCountData(count_data, exp_design = exp_design)
x <- calc_nf(x, norm_func = norm_rle, test_func = test_deseq2, full = ~ group + timepoint, reduced = ~ timepoint)
x@meta$nf

## 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 
##  0.8460345  0.9180475  0.8075463  0.9322872  0.9213414  0.9346221  0.8984660  0.8520127  0.8734221  0.9087962  0.9495643  0.7713711  1.0869889  0.9499309  0.9853505  0.6302187  1.0554200  1.0068082  1.0582351  1.1693629  1.1539737  1.0873290  0.9236152  1.1308167  1.0466599  1.1148338  1.2231097  1.1145235  1.1482744  1.2093569  1.0782560  1.0267380  1.0462877  1.0642993  1.0765078  0.9995915