Covariates and Batch correction

In this module, we will learn:

• How to identify possible “confounding” factors
• How batches or other covariates may impact your data

Differential Expression Workflow

Although we’ve successfully generated DE results, there are instances where patterns in our PCA plots or additional discussions with our collaborators cause us to revise our DE model(s).

In this case, we weren’t aware of any covariates that should be considered in our comparisons based on the information available publicly for these data. However, let’s look back at our initial PCA plot and see if there are any concerning patterns:

pca_plot = plotPCA(rld, intgroup = c('condition'), ntop = 500)
pca_plot

We do see a high amount of variance (58%) along PC1 that doesn’t seem to correspond to the experimental treatment. If we saw this pattern for a “real” dataset, we’d request additional information about the samples

How to model batch effects with DESeq2

Differences between samples can be due to biological covariates such as sex or patient of origin. Differences between samples can also be due to technical reasons, such as collection on different days or different sequencing runs. Groups of samples with different technical handling are called batches and differences due to the handling of samples are called batch effects.

Any relevant technical batches or biological covariates should be added as additional columns in the sample information table and added to a design formula. We can include batch effects and covariates in the same way in our design model, as long as the groups do not overlap, or confound, the biological treatment groups.

Let’s assume that after reviewing the above PCA plot, we reached out to our collaborators and they indicated that the samples were collected on three different days. We can add that information to our samplesheet and then create another DESeq2 object that includes those labels and a model with an additional term to account for possible batch effects due to date of collection.

samplesheet_batch = samplesheet_ready
samplesheet_batch$batch = factor(c(rep(c("Day1", "Day2", "Day3"), 2)), 
                                 levels = c("Day1", "Day2", "Day3"))

dds_batch = DESeqDataSetFromMatrix(countData = count_table,
                      colData = samplesheet_batch,
                      design = ~ batch + condition)

Note that we first created a new column called “batch” and added date of collection labels to our sample sheet and then initialized the DESeq2 object with that sample information and a model that includes an additional term that matches the “batch” column name.

While we don’t have time to delve into model options in detail, more complex model designs including adding “interaction terms” between multiple group labels, are helpfully described in this support thread as well as in the DESeq2 vignette.

Click for example of model design check

When including multiple terms in our model it’s helpful to check the corresponding design matrix to ensure that our batches are not confounded, which would cause DESeq2 to return a model not full rank error when attempting to fit the model. We can do that using the model.matrix function, providing our intended model and our sample information.

model.matrix(~ condition + batch, samplesheet_batch)

         (Intercept) conditiondeficient batchDay2 batchDay3
sample_A           1                  0         0         0
sample_B           1                  0         1         0
sample_C           1                  0         0         1
sample_D           1                  1         0         0
sample_E           1                  1         1         0
sample_F           1                  1         0         1
attr(,"assign")
[1] 0 1 2 2
attr(,"contrasts")
attr(,"contrasts")$condition
[1] "contr.treatment"

attr(,"contrasts")$batch
[1] "contr.treatment"

When we look at the outputs, we can see that all of the returned columns have values (1) included. If our model was not full rank, then we would see a columns with no values (all 0) returned.

Just like in our initial analysis, we’ll also want to filter our data.

keep = rowSums(counts(dds_batch)) >= 10
dds_batch_filtered = dds_batch[keep,]

Now that we’ve initialized a DESeq2 object that includes our batch labels and filtered the data, we can generate the rlog normalized count data and have that meta-data included:

rld_batch = rlog(dds_batch_filtered, blind = TRUE)
head(rld_batch)

class: DESeqTransform 
dim: 6 6 
metadata(1): version
assays(1): ''
rownames(6): ENSMUSG00000000001 ENSMUSG00000000028 ... ENSMUSG00000000049 ENSMUSG00000000056
rowData names(7): baseMean baseVar ... dispFit rlogIntercept
colnames(6): sample_A sample_B ... sample_E sample_F
colData names(4): genotype condition batch sizeFactor

Now that we have the additional labels added and normalized our data, we can generate a new PCA plot to see if the the batches explain any of the variance along PC1:

pca_plot_batch = plotPCA(rld_batch, intgroup = c('batch'), ntop = 500) + ggtitle("Batch labeled - iron deficiency data")
pca_plot_batch

From this plot, we can see that the dates of collection are primarily separated along PC1, suggesting that this variation might be batch effects. For data with multiple technical or biological covariates, we might need to generate PCA plots for each of the additional labels to determine what labels might be relevant to that dataset. However, based on the PCA plot for these data, we can proceed with fitting our model and generating DE results that account for these batches.

Fit a model that includes a covariate

To fit our updated model, just like for our initial analysis, we use the DESeq function:

dds_batch_fitted = DESeq(dds_batch_filtered)
resultsNames(dds_batch_fitted) 

[1] "Intercept"                      "batch_Day2_vs_Day1"             "batch_Day3_vs_Day1"            
[4] "condition_deficient_vs_control"

Notice that we have two additional sets of results that include the batch labels. While this indicates that the term we included for batch was fit, these results would not be biologically interesting.

Let’s look at these results from our comparison of interest:

resultsBatch_deficient_vs_control = results(dds_batch_fitted, name = 'condition_deficient_vs_control')
head(resultsBatch_deficient_vs_control) # results from model with term for batch

log2 fold change (MLE): condition deficient vs control 
Wald test p-value: condition deficient vs control 
DataFrame with 6 rows and 6 columns
                      baseMean log2FoldChange     lfcSE      stat    pvalue      padj
                     <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
ENSMUSG00000000001  1489.83039       0.274966  0.165703  1.659394 0.0970364  0.538526
ENSMUSG00000000028  1748.93544       0.217805  0.152329  1.429833 0.1527651  0.647029
ENSMUSG00000000031  2151.87715       0.138522  0.289999  0.477665 0.6328885  0.936072
ENSMUSG00000000037    24.91672       0.580793  0.709198  0.818944 0.4128186        NA
ENSMUSG00000000049     7.78377      -1.180919  1.360812 -0.867805 0.3855013        NA
ENSMUSG00000000056 19653.54030      -0.204334  0.180015 -1.135096 0.2563350  0.761406

head(results_deficient_vs_control) # results from model with no term

log2 fold change (MLE): condition deficient vs control 
Wald test p-value: condition deficient vs control 
DataFrame with 6 rows and 7 columns
                      baseMean log2FoldChange     lfcSE      stat    pvalue      padj     call
                     <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric> <factor>
ENSMUSG00000000001  1489.83039       0.297760  0.210310  1.415815  0.156830  0.868573       NS
ENSMUSG00000000028  1748.93544       0.226421  0.176795  1.280695  0.200301  0.902900       NS
ENSMUSG00000000031  2151.87715       0.457635  0.764579  0.598545  0.549476  0.995391       NS
ENSMUSG00000000037    24.91672       0.579130  0.561259  1.031840  0.302147  0.950613       NS
ENSMUSG00000000049     7.78377      -0.899483  1.553063 -0.579167  0.562476  0.998043       NS
ENSMUSG00000000056 19653.54030      -0.174048  0.203529 -0.855151  0.392468  0.982479       NS

We can see that while the structure of the results table is the same, the returned statistics are slightly different. However, since we didn’t actually talk to a collaborator to identify if batch could explain the % of variance observed within our treatment groups, we’ll save these results to file but will proceed with visualizing the results of our initial model.

save(dds_batch_fitted,
          file="outputs/Robjs/dds_batch_fitted.Robj")

---

Summary

In this section, we:

• Discussed technical batches and biological covariates
• Fitted a DESeq2 model that includes batch
• Generated tables of differential expression results for our batch model - i.e. fold changes and adjusted pvalues for each gene in dataset
• Saved our fitted model and results with batch to file

Now that we’ve generated multiple differential comparisons, we can determine how many genes are differentially expressed between our conditions and how to visualize our results.

---

---

Sources

• HBC DGE training module, part 1: https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html
• https://biohpc.cornell.edu/doc/RNA-Seq-2020-exercise2.html

---

These materials have been adapted and extended from materials listed above. These are open access materials distributed under the terms of the Creative Commons Attribution license (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

---

Previous lesson	Top of this lesson	Next lesson