1 Introduction to dupRadar

RNA-Seq experiments are a common strategy nowadays to quantify the gene expression levels in cells. Due to its higher sensitivity compared to arrays and the ability to discover novel features, makes them the choice for most modern experiments.

In RNA-Seq - as in all other NGS applications - quality control is essential. Current NGS workflows usually involve PCR steps at some point, which involves the danger of PCR artefacts and over-amplification of reads. For common DNA-based assays PCR duplicates are often simply removed before further analysis and their overall fraction or the read multiplicity taken as quality metrics. In RNA-Seq however, apart from the technical PCR duplicates, there is a strong source for biological read duplicates: for a given gene length and sequencing depth there exists an expression level beyond which it is not possible to add place more reads inside a gene locus without placing a second read exactly on the position of an already existing read. For this reason the overall duplication rate is not a useful measure for RNA-Seq.

As in the NGS/RNA-Seq QC ecosystem there were no suitable tools to address this question, we set out to develop a new tool. The dupRadar package gives an insight into the duplication problem by graphically relating the gene expression level and the duplication rate present on it. Thus, failed experiments can be easily identified at a glance.

Note: by now RNA-Seq protocols have matured so that for bulk RNA protocols data rarely suffer from technical duplicates. With newer low-input or single cell RNA-Seq protocols technical duplicates possible problems are worth to be checked for by default, especially if protocols are pushed to or beyond their boundaries. Paired-end libraries make the distinction between duplicates due to highly expressed genes and PCR duplicates easier, but the problem itself is not completely solved, especially at higher sequencing depths.

2 Getting started using dupRadar

library(dupRadar)

# call the duplicate marker and analyze the reads
bamDuprm <- markDuplicates(dupremover="bamutil",
               bam="test.bam",
               path="/opt/bamUtil-master/bin",
               rminput=FALSE)

if(suppressWarnings(require(AnnotationHub))) {
    ah = AnnotationHub()
    query(ah, c("ensembl", "80", "Takifugu", "gtf"))    # discovery
    cache(ah["AH47101"])                                # retrieve / file path
}

##                                 AH47101 
## "/home/biocbuild//.AnnotationHub/52579"

3 dupRate demo data

In the package we include two precomputed duplication matrices for two RNASeq experiments used as examples of a good and a failed (in terms of high redundancy of reads) experiments. The experiments come from the ENCODE project, as a source of a wide variety of protocols, library types and sequencing facilities.

Load the example dataset with:

attach(dupRadar_examples)

4 The duplication rate analysis

The analysis requires some info about the sequencing run and library preparation protocol:

• The strandess information about the reads: is the sequencing strand specific? if so, are the reads reversely sequenced?
• Are the reads paired, or single?

Due to its phenomenal performance, internally we use the featureCounts() function from the Bioconductor Rsubread package, which also supports multiple threads.

# The call parameters: 
bamDuprm <- "test_duprm.bam"    # the duplicate marked bam file
gtf <- "genes.gtf"  # the gene model
stranded <- 2       # '0' (unstranded), '1' (stranded) and '2' (reversely stranded)
paired   <- FALSE   # is the library paired end?
threads  <- 4       # number of threads to be used

# Duplication rate analysis
dm <- analyzeDuprates(bamDuprm,gtf,stranded,paired,threads)

The duplication matrix contains read counts in different scenarios and RPK and RPKM values for every gene.

5 Plotting and interpretation

The number of reads per base assigned to a gene in an ideal RNA-Seq data set is expected to be proportional to the abundance of its transcripts in the sample. For lowly expressed genes we expect read duplication to happen rarely by chance, while for highly expressed genes - depending on the total sequencing depth - we expect read duplication to happen often.

A good way to learn if a dataset is following this trend is by relating the normalized number of counts per gene (RPK, as a quantification of the gene expression) and the fraction represented by duplicated reads.

The dupRadar offers several functions to assess this relationship. The most prominent may be the 2d density scatter plot:

## make a duprate plot (blue cloud)
par(mfrow=c(1,2))
duprateExpDensPlot(DupMat=dm)       # a good looking plot
title("good experiment")
duprateExpDensPlot(DupMat=dm.bad)   # a dataset with duplication problems
title("duplication problems")