2.1 Preparing your data

Previous to running the duplication rate analysis, the BAM file with your mapped reads has to be duplicate marked with a like Picard, or the faster BamUtil dedup BioBamBam. The dupRadar package only works with duplicate marked BAM files.

If you do not have/want a duplication marking step in your default pipeline, the dupRadar package includes, for your convenience, wrappers to properly call some of these tools from your R session. Note that you still have to supply the path of the dupmarker installation though:

library(dupRadar)

Now, simply call the wrapper:

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

Simply specify which tool to use, the path where this tool is installed, and the input bam file to be analyzed. After marking duplicates, it’s safe to remove to original BAM file in order to save space.

The dupRadar package currently comes with support for:

• Picard MarkDuplicates
• BamUtil

After the BAM file is marked for duplicates, dupRadar is ready to analyze how the duplication rate is related with the estimated gene expression levels.

2.2 A GTF file

Unless there is any specific reason, dupRadar can use the same GTF file that will be later used to count the reads falling on features.

A valid GTF file can be obtained from UCSC, Ensembl, the iGenomes or other projects.

Note that the resulting duplication rate plots depend on the GTF annotation used. GTF files from the gencode projects result in a less clear picture of duplication rates, as there are many more features and feature types annotated, which overlap heavily as well. In some cases creating the plots only using subsets of gencode annotation files (e.g. just protein coding genes) serve the QC purpose of this tool better.

2.3 AnnotationHub as a source of GTF files

The Bioconductor AnnotationHub package provides an alternative approach to obtain annotation in GTF format from entities such as Ensembl, UCSC, ENCODE, and 1000 Genomes projects.

This is partly outlined in the AnnotationHub ‘HOWTO’ vignette section “Ensembl GTF and FASTA files for TxDb gene models and sequence queries”; for the Takifugu example, the downloaded GTF file is available from the cache.

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"

## 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")

## duprate boxplot
duprateExpBoxplot(DupMat=dm)        # a good looking plot
duprateExpBoxplot(DupMat=dm.bad)    # a dataset with duplication problems

5.1 Fitting a model into the data

To summarize the relationship between duplication rates and gene expression, we propose fitting a logistic regression curve onto the data. With the coefficients of the fitted model, one can get an idea of the initial duplication rate at low read counts (Intercept) and the progression of the duplication rate along with the progression of the read counts (Slope).

duprateExpDensPlot(DupMat=dm)

# or, just to get the fitted model without plot
fit <- duprateExpFit(DupMat=dm)
cat("duprate at low read counts: ",fit$intercept,"\n",
    "progression of the duplication rate: ",fit$slope,fill=TRUE)

## duprate at low read counts:  0.04075061 
##  
## progression of the duplication rate:  3.186793

Our main use case for that function is the condensation of the plots into quality metrics that can be used for automatic flagging of possibly problematic samples in large experiments or aggregation with other quality metrics in large tables to analyse interdependencies.

The duprateExpBoxplot plot shows the range of the duplication rates at 5% bins (default) along the distribution of RPK gene counts. The x-axis displays the quantile of the RPK distribution, and the average RPK of the genes contained in this quantile.

Individual genes can be identified in the plot:

## INTERACTIVE: identify single points on screen (name="ID" column of dm)
duprateExpPlot(DupMat=dm)       # a good looking plot
duprateExpIdentify(DupMat=dm)

One can also call the function duprateExpPlot to get smooth color density representation of the same data:

par(mfrow=c(1,2))
cols <- colorRampPalette(c("black","blue","green","yellow","red"))
duprateExpPlot(DupMat=dm,addLegend=FALSE)
duprateExpPlot(DupMat=dm.bad,addLegend=FALSE,nrpoints=10,nbin=500,colramp=cols)

Any further parm sent to the duprateExpPlot is also sent to the smoothScatter function.

5.2 Comparing the fitted parameters to other datasets

The parameters of the fitted model may mean very little (or just nothing) for many, unless there’s other data to compare with. We provide the pre-computed duplication matrices for all the RNA-Seq experiments publicly available from the ENCODE project, for human and mouse.

With the experience from the ENCODE datasets, we expect from single read experiments little duplication at low RPKM (low intercept) rapidly rising because of natural duplication (high slope). In contrast, paired-end experiments have a more mild rising of the natural duplication (low slope) due to having higher diversity of reads pairs since pairs with same start may still have different end.

The common denominator for both designs is the importance of having a low intercept, suggesting that duplication rate at lowly expressed genes may serve as a quality measure.

All the pre-computed duplication matrices are available in the dupRadar sourceforge site.

5.3 Other plots

CAVEAT: Sometimes in discussions duplicate reads (i.e. two physically different reads are mapped to the exact same position) and multi-mapping reads (i.e. a single read is mapped to two or more locations in the genome) are mixed up. dupRadar’s main focus are PCR duplicates in RNA-Seq. However internally we keep track of unique mappers and multimappers, and we use both in some of the examples, to illustrate use cases of our package beyond the main aim. Multi-mapping reads are completely independent of PCR duplicates.

Apart from the plots relating RPK and duplication rate, the dupRadar package provides some other useful plots to extract information about the gene expression levels.

An interesting quality metric are the fraction of reads taken up by groups of genes binned by 5% expression levels.

readcountExpBoxplot(DupMat=dm)

In the example we see that the 5% of highest expressed genes in our sample data set take up around 60% of all reads.

The distribution of RPK values per genes can be plotted with:

expressionHist(DupMat=dm)

This would help in identifying skewed distributions with unusual amount of lowly expressed genes, or to detect no consensus between replicates.

6.1 Fraction of multimappers per gene

Analogous to per gene duplication rate, the fraction of mutimappers can be easily calculated fom the duplication matrix.

Taking the counts from the column allCountsMulti, and substracting form it the counts from the column allCounts, one can get the total number of multihits. Thus, the fraction of multihits per gene can be calculating then dividing by allCountsMulti.

# calculate the fraction of multimappers per gene
dm$mhRate <- (dm$allCountsMulti - dm$allCounts) / dm$allCountsMulti
# how many genes are exclusively covered by multimappers
sum(dm$mhRate == 1, na.rm=TRUE)

## [1] 295

# and how many have an RPKM (including multimappers) > 5
sum(dm$mhRate==1 & dm$RPKMMulti > 5, na.rm=TRUE)

## [1] 8

# and which are they?
dm[dm$mhRate==1 & dm$RPKMMulti > 5, "ID"]

## [1] GPR89C       LOC728643    LOC606724    TBC1D3C      TBC1D3H     
## [6] MIR650       LOC100133050 HIST1H4J    
## 23228 Levels: 1/2-SBSRNA4 A1BG A1BG-AS1 A1CF A2LD1 A2M A2ML1 A2MP1 ... ZZZ3

We can also generate an overall picture about less extreme cases:

hist(dm$mhRate, 
     breaks=50, 
     col="red", 
     main="Frequency of multimapping rates", 
     xlab="Multimapping rate per gene", 
     ylab="Frequency")

Also the direct comparison of reads per kilobase between uniquely and multimappers is possible.

# comparison of multi-mapping RPK and uniquely-mapping RPK
plot(log2(dm$RPK), 
     log2(dm$RPKMulti), 
     xlab="Reads per kb (uniquely mapping reads only)",
     ylab="Reads per kb (all including multimappers, non-weighted)"
)

Use with identify() to annotate interesting cases interactively.

identify(log2(dm$RPK), 
         log2(dm$RPKMulti),
     labels=dm$ID)

6.2 Connection between possible PCR artefacts and GC content

In some cases we wondered about influence of GC content on PCR artefacts. An easy way to check using our dupRadar package in conjunction with biomaRt is demonstrated in the following. For simplicity we use our demo data here in which we do not see a big influence.

library(dupRadar)
library(biomaRt)

## for detailed explanations on biomaRt, please see the respective
## vignette

## set up biomart connection for mouse (needs internet connection)
ensm <- useMart("ensembl")
ensm <- useDataset("mmusculus_gene_ensembl", mart=ensm)

## get a table which has the gene GC content for the IDs that have been
## used to generate the table (depends on the GTF annotation that you
## use)
tr <- getBM(attributes=c("mgi_symbol", "percentage_gc_content"),
            values=TRUE, mart=ensm)

## create a GC vector with IDs as element names
mgi.gc <- tr$percentage_gc_content
names(mgi.gc) <- tr$mgi_symbol

## using dm demo duplication matrix that comes with the package
## add GC content to our demo data and keep only subset for which we can
## retrieve data
keep <- dm$ID %in% tr$mgi_symbol
dm.gc <- dm[keep,]
dm.gc$gc <- mgi.gc[dm.gc$ID]

## check distribution of annotated gene GC content (in %)
boxplot(dm.gc$gc, main="Gene GC content", ylab="% GC")

Now we can compare the dependence of duplication rate on expression level independently for below and above median GC genes (and to mention again, in this data set we don’t have a big difference).

par(mfrow=c(1,2))

## below median GC genes
duprateExpDensPlot(dm.gc[dm.gc$gc<=45,], main="below median GC genes")

## above median GC genes
duprateExpDensPlot(dm.gc[dm.gc$gc>=45,], main="above median GC genes")

#!/usr/bin/env Rscript

########################################
##
## dupRadar shell script
## call dupRadar R package from the shell for 
## easy integration into pipelines
##
## Holger Klein & Sergi Sayols
##
## https://sourceforge.net/projects/dupradar/
##
## input:
## - _duplicate marked_ bam file
## - gtf file
## - parameters for duplication counting routine:
##   stranded, paired, outdir, threads.
##
########################################

library(dupRadar)

####################
##
## get name patterns from command line
##
args   <- commandArgs(TRUE)

## the bam file to analyse
bam <- args[1]
## usually, same GTF file as used in htseq-count
gtf <- gsub("gtf=","",args[2])
## no|yes|reverse
stranded <- gsub("stranded=","",args[3])
## is a paired end experiment
paired   <- gsub("paired=","",args[4])
## output directory
outdir   <- gsub("outdir=","",args[5])
## number of threads to be used
threads  <- as.integer(gsub("threads=","",args[6]))

if(length(args) != 6) { 
  stop (paste0("Usage: ./dupRadar.sh <file.bam> <genes.gtf> ",
               "<stranded=[no|yes|reverse]> paired=[yes|no] ",
               "outdir=./ threads=1"))
}

if(!file.exists(bam)) { 
  stop(paste("File",bam,"does NOT exist"))
}

if(!file.exists(gtf)) {
  stop(paste("File",gtf,"does NOT exist"))
}

if(!file.exists(outdir)) { 
  stop(paste("Dir",outdir,"does NOT exist"))
}

if(is.na(stranded) | !(grepl("no|yes|reverse",stranded))) {
  stop("Stranded has to be no|yes|reverse")
}

if(is.na(paired) | !(grepl("no|yes",paired))) {
  stop("Paired has to be no|yes")
}

if(is.na(threads)) {
  stop("Threads has to be an integer number")
}

stranded <- if(stranded == "no") 0 else if(stranded == "yes") 1 else 2

## end command line parsing
##
########################################

########################################
##
## analyze duprates and create plots
##
cat("Processing file ", bam, " with GTF ", gtf, "\n")

## calculate duplication rate matrix
dm <- analyzeDuprates(bam,
                      gtf,
                      stranded,
                      (paired == "yes"),
                      threads)

## produce plots

## duprate vs. expression smooth scatter
png(file=paste0(outdir,"/",gsub("(.*)\\.[^.]+","\\1",basename(bam)),"_dupRadar_drescatter.png"),
    width=1000, height=1000)
duprateExpDensPlot(dm, main=basename(bam))
dev.off()

## expression histogram
png(file=paste0(outdir,"/",gsub("(.*)\\.[^.]+","\\1",basename(bam)),"_dupRadar_ehist.png"),
    width=1000, height=1000)
expressionHist(dm)
dev.off()

## duprate vs. expression boxplot
png(file=paste0(outdir,"/",gsub("(.*)\\.[^.]+","\\1",basename(bam)),"_dupRadar_drebp.png"),
    width=1000, height=1000)
par(mar=c(10,4,4,2)+.1)
duprateExpBoxplot(dm, main=basename(bam))
dev.off()

8.1 Citing dupRadar

Please consider citing dupRadar if used in support of your own research:

citation("dupRadar")

## 
##   Sergi SayolsEmail author, Denise Scherzinger and Holger Klein
##   (2016): dupRadar: a Bioconductor package for the assessment of
##   PCR artifacts in RNA-Seq data. BMC Bioinformatics, 17:428,
##   doi:10.1186/s12859-016-1276-2
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     title = {dupRadar: a Bioconductor package for the assessment of PCR artifacts in RNA-Seq data},
##     author = {Sergi Sayols and Denise Scherzinger and Holger Klein},
##     year = {2016},
##     journal = {BMC Bioinformatics},
##     doi = {10.1186/s12859-016-1276-2},
##     url = {http://dx.doi.org/10.1186/s12859-016-1276-2},
##     volume = {17},
##     issue = {1},
##     pages = {428},
##   }

8.2 Reporting problems or bugs

If you run into problems using dupRadar, the Bioconductor Support site is a good first place to ask for help. If you think there is a bug or an unreported feature, you can report it using the dupRadar sourceforge site.