Some RNA-seq data might suffer from genomic DNA (gDNA) contamination due to carrying over of residual gDNA in RNA preparation into sequencing library. The R package CleanUpRNAseq provides a set of functionalities to detect and correct for gDNA contamination, thus facilitate better quantitation of gene expression and differential expression analysis.
CleanUpRNAseq 0.99.15
RNA-seq has become a state-of-art technology for studying gene expression (Wang, Gerstein, and Snyder 2009). However, due to improper RNA preparation and choice of some library preparation protocols, such as rRNA-depletion based RNA-seq protocol (O’Neil, Glowatz, and Schlumpberger 2013) and the SMART-Seq protocol, RNA-seq data might suffer from genomic DNA (gDNA) contamination, which skews quantitation of gene expression and hinders differential gene expression analysis (Li et al. 2022; Verwilt et al. 2020; Markou et al. 2021). Some quality control tools have been developed to check the quality of RNA-seq data at the raw read and post-alignment levels, including FastQC, RSeQC (Wang, Wang, and Li 2012), Qualimap (GarcÃa-Alcalde et al. 2012), RNA-SeQC/RNA-SeQC 2 (DeLuca et al. 2012; Graubert et al. 2021), QoRTs (Hartley and Mullikin 2015), RNA-QC-Chain (Zhou et al. 2018), and RNAseqQC. Those post-alignment tools can report percentages of reads mapping to different genomic features, such as genes, exons, introns, intergenic regions, and rRNA exons. Thus, they can be used to roughly detect gDNA contamination. So far, SeqMonk and the gDNAx package are the only tool which can be used to both detect and correct for gDNA contamination in RNA-seq data. However, SeqMonk is a Java-based GUI tool which makes it not available in most high performance computing cluster. More importantly, seqMonk assumes a uniform distribution of reads derived from gDNA contamination and its performance on correcting for gDNA contamination in RNA-seq data is not fully evaluated and peer reviewed. On the other hand, gDNAx is an R/Bioconductor package and can be incorporated into automated RNA-seq data analysis pipeline easily. However, gDNAx only simply removes reads mapping to intergenic regions and introns, but not those mapping to exons, to mitigate gDNA contamination. Thus, gDNAx’s algorithm for correcting for gDNA contamination is not sophisticated. To this end, we developed the R pacakge CleanUpRNAseq, which provides a full set of functions for detecting and correcting for gDNA contamination in RNA-seq data across all computing platforms.
As for any other sequencing data analysis, users should first check the quality of raw RNA-seq sequencing data using a tool like FastQC followed by MultiQC (Ewels et al. 2016). Depending on the quality of the raw data, trimming low quality bases and/or adpator sequences might be necessary to increase mapping quality and rate. Subsequently, reads are mapped to the reference genome of interest using tools, such as STAR (Dobin et al. 2012), HISAT2 (Kim et al. 2019), GMAP (Wu et al. 2016), and Subread (Liao, Smyth, and Shi 2013). The resulting alignment files in the BAM format are used for post-alignment quality control, including detection of gDNA contamination. Meanwhile, if the RNA-seq library is unstranded, reads are mapped to the reference transcriptome using Salmon (Patro et al. 2017) to get transcript-level abundance, counts, and effective lengths, since Salmon can resolve reads mapped to multiple transcripts, which are imported into R using the tximport to get gene-level counts, abundance and length information. However, if the library is stranded, reads are mapped to the reference transcriptome using Salmon twice: one using the designed strandedness for the ‘–libType’ option, the other using the opposite strandedness. See Salmon library type discussion https://salmon.readthedocs.io/en/latest/library_type.html. Then gene-level counts, abundance and length information are imported into R as above.
First, load the required packages, including CleanUpRNAseq. Then users can perform a step-by-step analysis of the RNA-seq data as below for more flexibility. For users’ convenience, CleanUpRNAseq also offers two wrapper functions, create_diagnostic_plots and correct_for_contamination, for detecting and correcting for gDNA contamination in RNA-seq data. As for how to use these wrappers, please refer to their function documentation.
suppressPackageStartupMessages({
library("CleanUpRNAseq")
#devtools::load_all("../../CleanUpRNAseq")
library("ggplotify")
library("patchwork")
library("ensembldb")
library("utils")
})
Users can list all current EnsDb packages from AnnotationHub and load the package of choice, if available, as follows. Here, an EnsDb package for the human genome is loaded. It is crucial to use an EnsDb pacakge which represents the genome annotation file (GTF) used for RNA-seq read alignment.
suppressPackageStartupMessages({
library("EnsDb.Hsapiens.v86")
})
hs_ensdb_sqlite <- EnsDb.Hsapiens.v86
Otherwise, users can easily prepare an EnsDb database from an Ensembl GTF file. For all following steps, this option is adopted because the latest version of human transcriptome has been used for read mapping by STAR and Salmon to quantify gene expression of the example RNA-seq data.
options(timeout = max(3000, getOption("timeout")))
tmp_dir <- tempdir()
gtf <- system.file("extdata", "example.gtf.gz",
package = "CleanUpRNAseq")
hs_ensdb_sqlite <-
ensDbFromGtf(
gtf = gtf,
outfile = file.path(tmp_dir, "EnsDb.hs.v110.sqlite"),
organism = "Homo_Sapiens",
genomeVersion = "GRCh38",
version = 110
)
Potential DNA contamination exists if a significantly high portion of RNA-seq reads mapped to intergenic regions. CleanUpRNAseq uses the featureCounts function from the Rsubread package to quantify the percentage of reads mapped to different genomic features: genes, exons, introns, intergenic regions, rRNA exons, and organelle genome(s). So users need to prepare SAF files for these genomic features. In addition, a BAM file is needed to provide the lengths of the chromosomes/scaffolds.
bam_file <- system.file("extdata", "K084CD7PCD1N.srt.bam",
package = "CleanUpRNAseq"
)
saf_list <- get_saf(
ensdb_sqlite = hs_ensdb_sqlite,
bamfile = bam_file,
mitochondrial_genome = "MT"
)
Reads mapped to different genomic features is summarized by using the featureCounts function with the SAF files generated above as annotation. CleanUpRNAseq also quantifies the reads spanning exon-exon junctions and the reads mapped to exons of each gene using the GTF file as annotation. The junction read count table is used to determine the unexpressed, multiexonic genes, while the gene-level read count table is used for comparing samples at the gene level. Here two downsampled BAM files are used for speeding up the demonstration, while a precomputed R object in the .RData format is used for actual downstream analysis.
tmp_dir <- tempdir()
in_dir <- system.file("extdata", package = "CleanUpRNAseq")
gtf.gz <- dir(in_dir, ".gtf.gz$", full.name = TRUE)
gtf <- file.path(tmp_dir, gsub("\\.gz", "", basename(gtf.gz)))
R.utils::gunzip(gtf.gz, destname= gtf,
overwrite = TRUE, remove = FALSE)
in_dir <- system.file("extdata", package = "CleanUpRNAseq")
BAM_file <- dir(in_dir, ".bam$", full.name = TRUE)
salmon_quant_file <- dir(in_dir, ".sf$", full.name = TRUE)
sample_name = gsub(".+/(.+?).srt.bam", "\\1", BAM_file)
salmon_quant_file_opposite_strand <- salmon_quant_file
col_data <- data.frame(sample_name = sample_name,
BAM_file = BAM_file,
salmon_quant_file = salmon_quant_file,
salmon_quant_file_opposite_strand =
salmon_quant_file_opposite_strand,
group = c("CD1N", "CD1P"))
sc <- create_summarizedcounts(lib_strand = 0, colData = col_data)
## featurecounts
capture.output({counts_list <- summarize_reads(
SummarizedCounts = sc,
saf_list = saf_list,
gtf = gtf,
threads = 1,
verbose = TRUE
)}, file = tempfile())
## load salmon quant
salmon_counts <- salmon_tximport(
SummarizedCounts = sc,
ensdb_sqlite = hs_ensdb_sqlite
)
Precomputed R object in the .RData format (feature_counts_list.rda, and
salmon_quant.rda)
containing the featureCounts output and Salmon quantification output
imported by using the tximport
function for a RNA-seq dataset of 8 samples
from two treatment groups are used for demo. GC-content for genes and
intergenic regions (intergenic_GC.rda, and gene_GC.rda) were also
precomputed and used for demo.
data("feature_counts_list")
data("salmon_quant")
data("intergenic_GC")
data("gene_GC")
lib_strand <- 0
col_data_f <- system.file("extdata", "example.colData.txt",
package = "CleanUpRNAseq")
col_data <- read.delim(col_data_f, as.is = TRUE)
## create fake bam files
tmp_dir <- tempdir()
bamfiles <- gsub(".+/", "", col_data$BAM_file)
null <- lapply(file.path(tmp_dir, bamfiles), file.create)
## create fake quant.sf files
quant_sf <- file.path(tmp_dir, gsub(".srt.bam$",
"quant.sf",
bamfiles))
null <- lapply(quant_sf, file.create)
col_data$BAM_file <- file.path(tmp_dir, bamfiles)
col_data$salmon_quant_file <- quant_sf
## pretend this is stranded RA=NA-seq data
col_data$salmon_quant_file_opposite_strand <- quant_sf
sc <- create_summarizedcounts(lib_strand, col_data)
sc$set_feature_counts(feature_counts_list)
sc$set_salmon_quant(salmon_quant)
sc$set_salmon_quant_opposite(salmon_quant)
p_assignment_stat <-plot_assignment_stat(SummarizedCounts = sc)
wrap_plots(p_assignment_stat)
assigned_per_region <- get_region_stat(SummarizedCounts = sc)
p <- plot_read_distr(assigned_per_region)
p