Integrating bulk RNA-Seq data with different sequencing depths and from different sources

Question

I am attempting to integrate different bulk RNA-Seq datasets. While this is not ideal, I'm trying to reduce the technical variability in these datasets by using data generated by similar protocols (for example, all of them have been sequenced with a Poly-A enrichment strategy and Illumina short read sequencing). I am aware the results I obtain from this analysis will be very limited, but still want to pursue the analysis. I'm looking for advice on two aspects of the analysis, but first some more info about the data:

• There are 5 bulk RNA-Seq datasets and 4 biological groups. Each datasets only comprises one biological group. Only two biological groups have two datasets, and such datasets overlap in PCA.
• In Principal Component Analysis (PCA) and the first PC (~26% of variation) is very correlated with the biological variable (cell differentiation), suggesting there is a biological signal among all the noise/technical bias.

Despite these preliminary results I'm not sure I'm accounting correctly for some aspects: The datasets have been sequenced with different sequencing coverage (in fact, different numbers and length of reads). Some are single read and some paired-end. My questions are:

1. For PCA, I'm normalizing with variance stabilizing transformation (VST). How should I account for the differences in read numbers and type between datasets?

2. For differential expression (DE) analysis I'm using DESeq2. How should I account for the differences in read numbers and type between datasets? There is no batch variable, as all datasets come from different sources (batches, laboratories)

I'm happy to provide more information if necessary. Thanks

Answer 1

Per-sample coverage / count differences are expected in all differential expression analysis, and correction for that is built into DESeq2. You shouldn't need to do any additional correction with respect to proportional counts when looking at a PCA derived from the VST.

For accounting for different sources / types / etc, it's a good idea to add that information into the metadata, and test for differential expression with and without accounting for those variables to see if anything changes. This is recommended in the DESeq2 manual:

If there is unwanted variation present in the data (e.g. batch effects) it is always recommend to correct for this, which can be accommodated in DESeq2 by including in the design any known batch variables or by using functions/packages such as svaseq in sva (Leek 2014) or the RUV functions in RUVSeq (Risso et al. 2014) to estimate variables that capture the unwanted variation. In addition, the ashr developers have a specific method for accounting for unwanted variation in combination with ashr (Gerard and Stephens 2017).

See the multi-factor design section of the manual for more information (which has a "single-read" vs "paired-end" example):

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs

If after accounting for the different variables there is not enough information to properly separate groups, DESeq should tell you. Ideally you should aim for six replicates of each state of each independent variable (easier to do when there are interacting variables), but it's quite rare to see that in differential expression analysis.

Answer 2

In case of different sequencing depth one performs rarefaction, to make the data comparable.

The standard approach is rarefaction by size, when one takes the number of reads in the smallest sample, and samples without replacement the same number of reads from the other samples. Potentially, one could carry out upsizing to a greater number of reads (e.g., if the smallest sample is an outlier in terms of the number of reads, and downsizing may result in losing too much information.) I suggest reading the tutorial for MOMR package for more detailed description (the R code is also available.)

There are more sophisticated methods of rarefaction, such as rarefaction by coverage - these are of significant scientific interests, but tricky to use and not very widespread. See, e.g., this article (it deals with sampling species, but generalization to reads is straightforward... though reads are not species.)

Remark: Of course many good packages (like those mentioned by @gringer) might already perform rarefaction by default - but it is always a good idea to check the documentation.