High throughput data such as RNAseq is often complicated by technical errors that are not explicitly modeled in the experimental design and can confound the detection of differential expression. Differences in counts from samples run on different days or different locations or on different machines are an example of a technical error that is very often present and which we should try to model in our experimental design. An approach to this is to build a surrogate variable into our experimental design that explains the batch effect and take it into account in the modeling and differential expression analysis stages. We'll use the SVA package to do this.

We'll need the SVA package from Bioconductor and our Arabidopsis count data in datasets/ch1/arabidopsis.RDS.

For estimating batch effects using SVA, perform the following steps:

1. Load the libraries and data:

library(sva)
arab <- readRDS(file.path(getwd(), "datasets", "ch1", "arabidopsis.RDS"))

2. Filter out rows with too few counts in some experiments:

keep <- apply(arab, 1, function(x) { length(x[x>3])>=2 } )
arab_filtered <- arab[keep,]

3. Create the initial design:

groups <- as.factor(rep(c("mock", "hrcc"), each=3))

4. Set up the test and null models and run SVA:

test_model <- model.matrix(~groups)
null_model <- test_model[,1]
svar <- svaseq(arab_filtered, test_model, null_model, n.sv=1)

5. Extract the surrogate variables to a new design for downstream use:

design <- cbind(test_model, svar$sv)

In step 1, the script begins by loading in a count matrix of the Arabidopsis RNAseq data, which you will recall is a simple experiment with three replicates of mock and three of hrcc treatment.

In step 2, we create a vector of row indices that we wish to retain, which we do by testing whether the row has at least two columns with a count of over 3—this is done by using apply() and an anonymous function over the rows of the count matrix.

With step 3, we create a groups factor describing the experiment sample types.

Step 4 is the one that does the work; we use the groups factor in model.matrix() to create the model design we want to use. We also need a null model, which, in this experimental design, is equivalent to the first column, so we extract that from the test_model design object. We can then use the key svaseq() function to estimate the surrogate variable to add to our design. We add in test_model and null_model and select the number of surrogate variables to use with n.sv, which should be one for a simple design like this.

The final bit, step 5, is to add the surrogate variable to the design model, which we do by binding test_model and the sv column of svar (svsar$sv) together. The final design object can then be used in packages such as edgeR and DESeq2 as any other and those methods will use the surrogate variable to deal with batch effects.