Batch effects and GC content
Michael Love
Previously, we looked at distances between samples in high-throughput sequencing experiments, exploring sequencing depth as a technical artifact, the effect of various transformation on stabilizing variance, and PCA and hierarchical clustering as methods for ordination of samples.
In distances, we noticed many genes where we see differences in measurements across the sequencing center, and in hclust, we saw high level clustering of samples by sequencing center, when we focused on a single human population. These technical differences in high-throughput measurements are often referred to as “batch effects”, a term which encompasses any kind of technical artifact due to the lab or time at which a sample was processed. Note that the time of sample preparation, the batch, can have its own unique technical distortion on measurements, just as much as a sample from a different sequencing center.
Two references which give a general outline of batch effects in sequencing data are:
- Tackling the widespread and critical impact of batch effects in high-throughput data
- Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries
It is critical, if you are involved in the design of a high-throughput experiment, that the biological conditions are not confounded with the sample preparation batches, instead to use either a block or randomized design. Block designs are easy, simply ensuring that each sample preparation batch contains each of the biological conditions of interest, so that the batch effects can be isolated apart from the biological differences between samples.
Here I will show the origin of batch effects across lab can sometimes be isolated to batch-specific distortions related to the DNA sequence of the genes. As described in the second link above, a step in nearly all high-throughput experiments is to amplify the DNA fragments using PCR. As PCR is an exponential process of copying molecules, it is very sensitive to slight variations, and results in distortions of measurements which are specific to the particular place and time when the samples were processed: the batch.
GEUVADIS dataset
We start again with the GEUVADIS RNA-seq samples prepared in distances.
library(DESeq2)
load("geuvadis.rda")As typical, though we conceptually have a simple task to perform, many steps are needed just to get the covariates that we need. Here our simple task is to model the counts on covariates like the genes’ lengths and sequence content. However, it will take many steps to obtain and summarize the sequence content into a single number.
A critical summary of DNA sequence is the GC content, which is the fraction of G’s and C’s out of the total sequence. Note that, because G pairs with C, and C with G, the GC content of a sequence and it’s complement is identical. GC content varies along the genome, for biological reasons. Different organisms have very different GC content of their genomes. But we care about GC content also for a technical reason: because the GC content of a piece of DNA can make it difficulty to amplify. We will see this below.
In order to get the GC content of each gene, we start by downloading the information about the genes used by recount2 for quantification, which was version 25 of the Gencode genes. We download this file from the following link (37 Mb).
library(GenomicFeatures)We need to import this gene information from a GTF (gene transfer format) file, and turn it into a TxDb. We then save the database, so we can skip the makeTxDbFromGFF step in the future. Note that a warning message about The "phase" metadata column contains non-NA values, and another one about closing connection are both OK.
txdb <- makeTxDbFromGFF("gencode.v25.annotation.gtf.gz")
saveDb(txdb, file="gencode.sqlite")We load the TxDb and extract all the exons, grouped by gene. We check that we have the same gene names in ebg and in dds.
txdb <- loadDb("gencode.sqlite")
ebg <- exonsBy(txdb, by="gene")
head(names(ebg))## [1] "ENSG00000000003.14" "ENSG00000000005.5" "ENSG00000000419.12" "ENSG00000000457.13"
## [5] "ENSG00000000460.16" "ENSG00000000938.12"table(names(ebg) %in% rownames(dds))##
## TRUE
## 58037table(rownames(dds) %in% names(ebg))##
## TRUE
## 58037Note that the exons in ebg contain redundant sequence. We can see this by plotting the ranges for a given gene. Note that, after running reduce, we remove all redundant sequence in ebg for a given gene.
e <- ebg[[1]]
library(rafalib)
plotRanges <- function(e) {
l <- length(e)
r <- ranges(range(e))
nullplot(start(r), end(r), 0, l+1)
segments(start(e), 1:l, end(e), 1:l, lwd=5)
}
plotRanges(e)
plotRanges(reduce(e))
Now we put the reduced exons in correct order (in this case they are already in correct order), and we store them as rowRanges of the dataset.
exons <- reduce(ebg)
exons <- exons[ rownames(dds) ]
rowRanges(dds) <- exonsCalculate GC content and length of reduced exons
Now we extract the exonic sequence for every gene, using the extractTranscriptSeqs function.
library(BSgenome.Hsapiens.UCSC.hg38)
dna <- extractTranscriptSeqs(Hsapiens, rowRanges(dds))We then calculate the GC content (ratio of G or C to total basepairs) with letterFrequency, and save this as a metadata column gc. We also save the total number of basepairs to a metadata column len.
mcols(dds)$gc <- as.numeric(letterFrequency(dna, "GC", as.prob=TRUE))
mcols(dds)$len <- sum(width(rowRanges(dds)))with(mcols(dds), hist(gc))
with(mcols(dds), hist(log10(len)))
We know have all the covariates we need for modeling how the counts vary by GC content. We can make simple plots to see if we see a dependence. Note that, outside of GC content of .35-.65, we see very few large counts, although there do appear to be genes with this GC content. It is very difficult to amplify the fragments of cDNA from these genes, and so they are often missing from high-throughput sequencing experiments like RNA-seq.
plot(mcols(dds)$gc, log10(counts(dds)[,1]+1), cex=.1)
abline(v=c(.35,.65))
We also plot the counts over the length of the gene. It is expected that, everything else being equal, we should see higher counts from longer genes. Keep in mind though, that we do not expect a line in this scatterplot, due to differences in gene expression of the genes at a given length. In other words, any point on the x-axis, the differences in a vertical band can be explained by gene expression (as well as any other technical covariates, like GC content).
plot(log10(mcols(dds)$len), log10(counts(dds)[,1]+1), cex=.1)
Conditional Quantile Normalization
We will now use a Bioconductor package called cqn to model systematic dependence of counts on gene GC content and length. We provide the cqn function with the counts, the GC content and the gene lengths.
You can ignore a warning about use of 'sig2' is deprecated...
library(cqn)
idx <- dds$population == "TSI"
dds2 <- dds[,idx]
cts <- counts(dds2)
fit <- cqn(cts, mcols(dds2)$gc, mcols(dds2)$len)## fitting ...
##
|
| | 0%
|
|============================== | 33%
|
|============================================================ | 67%
|
|==========================================================================================| 100%## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2)) insteadThe plots show estimated spline dependence of counts on GC (n=1)…
cqnplot(fit, n=1)
…and dependence on length (n=2).
cqnplot(fit, n=2, xlim=c(-2,4.5))
Both of these plots are typical: for GC content, we see an upside-down “U”, where the low and high GC content fragments are systematically underrepresented (although we see lots of sample-sample variability in the estimated splines). And it is typical to have more counts for longer genes due to fragmentation. The tick marks on the x-axis indicate the knots of the splines, by default, these are data quantiles: 0.025, 0.25, 0.5, 0.75, 0.975.
We can draw the lines with the sequencing center as the color. Doing this, we see that the lines cluster by sequencing center. The dependence of counts on the GC content of the features being sequenced is highly batch specific.
library(rafalib)
bigpar()
cqnplot(fit, n=1, col=dds2$center)
legend("bottom", levels(dds2$center), fill=1:nlevels(dds2$center))
There is less variation of counts on gene length across center:
bigpar()
cqnplot(fit, n=2, col=dds2$center, xlim=c(-2,4.5))
legend("topleft", levels(dds2$center), fill=1:nlevels(dds2$center))
Finally, we zoom in to two sequencing centers, to show how different the dependence of counts on GC content can be across batch. You can also cross reference this final plot with Figure 2 here, which goes into more depth on the topic.
idx <- dds$population == "TSI" & dds$center %in% c("CGR","UG")
dds3 <- dds[,idx]
cts <- counts(dds3)
fit <- cqn(cts, mcols(dds3)$gc, mcols(dds3)$len)## fitting ...
##
|
| | 0%
|
|============================== | 33%
|
|============================================================ | 67%
|
|==========================================================================================| 100%## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2)) insteaddds3$center <- droplevels(dds3$center)
bigpar()
cqnplot(fit, n=1, col=dds3$center)
legend("bottom", levels(dds3$center), fill=1:nlevels(dds3$center))
Downstream use of CQN
Above we show plots of the bias of the counts of sequenced reads over aspects of the features (genes) include sequence composition (GC) and the number of basepairs in the reduced exon ranges. The cqn output is also valuable for plotting corrected data, and for providing offsets for use in downstream statistical analysis.
names(fit)## [1] "counts" "lengths" "sizeFactors" "subindex" "y" "x"
## [7] "offset" "offset0" "glm.offset" "func1" "func2" "grid1"
## [13] "grid2" "knots1" "knots2" "call"A natural log scale offset is provided in glm.offset. This can be converted back to count scale (similar to a size factor but across all genes x samples):
cqnOffset <- fit$glm.offset
cqnNormFactors <- exp(cqnOffset)The cqn output also has count values that are corrected for differential bias across samples. We can compare before…
bigpar()
filter <- as.integer(rowSums(counts(dds3)) > 10 & mcols(dds)$len > 1000)
head(order(mcols(dds3)$gc * filter,
decreasing=TRUE), 10)## [1] 15820 32706 5674 53583 16699 15239 14765 57124 57134 57177idx <- 32706
mcols(dds)[idx,]## DataFrame with 1 row and 2 columns
## gc len
## <numeric> <integer>
## ENSG00000234965.2 0.743772 1405boxplot(counts(dds3, normalized=TRUE)[idx,] ~ dds3$center,
ylim=c(0,300), col=1:2,
xlab="center", ylab="scaled counts")
to after:
bigpar()
exprs <- fit$y + fit$offset # note: atypical definition of offset...
boxplot(exprs[idx,] ~ dds3$center, col=1:2,
xlab="center", ylab="log2 scale expression")
While not all of the high GC features are exactly corrected, we see that for a number of them, the differences have been accounted for by modeling on technical covariates.
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] splines parallel stats4 stats graphics grDevices datasets utils methods
## [10] base
##
## other attached packages:
## [1] cqn_1.36.0 quantreg_5.75
## [3] SparseM_1.78 preprocessCore_1.52.0
## [5] nor1mix_1.3-0 mclust_5.4.7
## [7] BSgenome.Hsapiens.UCSC.hg38_1.4.3 BSgenome_1.58.0
## [9] rtracklayer_1.50.0 Biostrings_2.58.0
## [11] XVector_0.30.0 rafalib_1.0.0
## [13] GenomicFeatures_1.42.1 AnnotationDbi_1.52.0
## [15] DESeq2_1.30.0 SummarizedExperiment_1.20.0
## [17] Biobase_2.50.0 MatrixGenerics_1.2.0
## [19] matrixStats_0.57.0 GenomicRanges_1.42.0
## [21] GenomeInfoDb_1.26.2 IRanges_2.24.1
## [23] S4Vectors_0.28.1 BiocGenerics_0.36.0
## [25] testthat_3.0.1 rmarkdown_2.6
## [27] devtools_2.3.2 usethis_2.0.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ellipsis_0.3.1 rprojroot_2.0.2
## [4] fs_1.5.0 MatrixModels_0.4-1 remotes_2.2.0
## [7] bit64_4.0.5 fansi_0.4.1 xml2_1.3.2
## [10] codetools_0.2-18 geneplotter_1.68.0 knitr_1.30
## [13] pkgload_1.1.0 Rsamtools_2.6.0 annotate_1.68.0
## [16] dbplyr_2.0.0 compiler_4.0.3 httr_1.4.2
## [19] assertthat_0.2.1 Matrix_1.3-0 cli_2.2.0
## [22] htmltools_0.5.0 prettyunits_1.1.1 tools_4.0.3
## [25] gtable_0.3.0 glue_1.4.2 GenomeInfoDbData_1.2.4
## [28] dplyr_1.0.2 rappdirs_0.3.1 Rcpp_1.0.5
## [31] vctrs_0.3.6 conquer_1.0.2 xfun_0.19
## [34] stringr_1.4.0 ps_1.5.0 lifecycle_0.2.0
## [37] XML_3.99-0.5 zlibbioc_1.36.0 scales_1.1.1
## [40] hms_0.5.3 RColorBrewer_1.1-2 yaml_2.2.1
## [43] curl_4.3 memoise_1.1.0 ggplot2_3.3.3
## [46] biomaRt_2.46.2 stringi_1.5.3 RSQLite_2.2.1
## [49] genefilter_1.72.0 desc_1.2.0 pkgbuild_1.2.0
## [52] BiocParallel_1.24.1 rlang_0.4.9 pkgconfig_2.0.3
## [55] bitops_1.0-6 evaluate_0.14 lattice_0.20-41
## [58] purrr_0.3.4 GenomicAlignments_1.26.0 bit_4.0.4
## [61] processx_3.4.5 tidyselect_1.1.0 magrittr_2.0.1
## [64] R6_2.5.0 generics_0.1.0 DelayedArray_0.16.0
## [67] DBI_1.1.0 pillar_1.4.7 withr_2.3.0
## [70] survival_3.2-7 RCurl_1.98-1.2 tibble_3.0.4
## [73] crayon_1.3.4 BiocFileCache_1.14.0 progress_1.2.2
## [76] locfit_1.5-9.4 grid_4.0.3 blob_1.2.1
## [79] callr_3.5.1 digest_0.6.27 xtable_1.8-4
## [82] openssl_1.4.3 munsell_0.5.0 sessioninfo_1.1.1
## [85] askpass_1.1