How to normalise Bulk RNA-seq data to account for transcript length, coverage depth, library size and batch effects?

Question

I'm currently working on a project where I'm comparing aggregate measurements (mean, median, etc.) of expression data (RNA-seq) from different groups of genes across various samples with different characteristics (tissue type, health status, etc.). Additionally, the raw counts were collected from several different labs using various techniques.

Since I am conducting between-gene measurements, the data should be normalised to account for differences in transcript length and coverage depth (TPM, RPKM, FPKM). However, I am also interested in comparisons across samples based on tissue type and other factors. Therefore, the data should also be normalised to account for library size (TMM, quantile, etc.), and, as the data were collected from multiple sources, it should be corrected for batch effects.

I have read through many papers but am unsure and confused about how to proceed with the normalisation procedure starting with the raw counts. Can I simply string the methods together, starting with batch effect correction, followed by library size normalisation, and then the within-sample normalisations?

I would appreciate any insights or suggestions on this. Thanks

Answer 1

All factors you mention bar the batch effects can be addressed by some sort of normalization (I use a modified version of UQpgQ2). Batch effects can be addressed using ComBat/SVA/ComBat-seq but my team's experience with them has been bad - they introduce artifacts and skew the data quite a bit especially if the batch effect is not pronounced.

Unfortunately, there is no way to do simple comparisons without huge caveats. You're better off using DESeq2 to perform statistically sound comparisons using the batch variable as a covariate.

Answer 2

Have a look through the DESeq2 manual. It goes into a lot of detail about how it works, and how you can use it to do all the things you're asking about:

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

The standard workflow already has all these things strung together, and works fairly well for the most common differential expression tasks:

Quick start

Here we show the most basic steps for a differential expression analysis. There are a variety of steps upstream of DESeq2 that result in the generation of counts or estimated counts for each sample, which we will discuss in the sections below. This code chunk assumes that you have a count matrix called cts and a table of sample information called coldata. The design indicates how to model the samples, here, that we want to measure the effect of the condition, controlling for batch differences. The two factor variables batch and condition should be columns of coldata.

dds <- DESeqDataSetFromMatrix(countData = cts,
                             colData = coldata,
                             design= ~ batch + condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_trt_vs_untrt")
# or to shrink log fold changes association with condition:
res <- lfcShrink(dds, coef="condition_trt_vs_untrt", type="apeglm")