inst/doc/RW3_BAM_QCs.R

In ramwas: Fast Methylome-Wide Association Study Pipeline for Enrichment Platforms

## ----global_options, include=FALSE--------------------------------------------
#getwd()
#knitr::opts_chunk$set(fig.align="center", fig.retina=1)
knitr::opts_chunk$set(fig.retina=1)
library(ramwas)

## ----loadCgGset---------------------------------------------------------------
library(ramwas)
filename = system.file("extdata", "bigQC.rds", package = "ramwas")
qc = readRDS(filename)$qc

## ----nbams--------------------------------------------------------------------
cat("Number of BAMs:", qc$nbams)

## ----reads.total--------------------------------------------------------------
cat("Reads total:", qc$reads.total)

## ----reads.aligned------------------------------------------------------------
{
cat("Reads aligned:", qc$reads.aligned, "\n")
cat("This is ", qc$reads.aligned / qc$reads.total * 100,
    "% of all reads", sep="")
}

## ----reads.recorded-----------------------------------------------------------
{
cat("Reads recorded:",qc$reads.recorded,"\n")
cat("This is ", qc$reads.recorded / qc$reads.aligned * 100,
    "% of aligned reads", sep="")
}

## ----reads.recorded.no.repeats------------------------------------------------
{
cat("Reads without duplicates:", qc$reads.recorded.no.repeats, "\n")
cat("This is ", qc$reads.recorded.no.repeats / qc$reads.recorded * 100,
    "% of aligned reads", "\n", sep="")
}

## ----frwrevNR-----------------------------------------------------------------
{
    cat("Excluding duplicate reads", "\n")
    cat("Reads on forward strand:", qc$frwrev.no.repeats[1], "\n")
    cat("Reads on reverse strand:", qc$frwrev.no.repeats[2], "\n")
    cat("Fraction of reads on forward strand:",
        qcmean(qc$frwrev.no.repeats), "\n")
}

## ----frwrev-------------------------------------------------------------------
{
    cat("Not excluding duplicate reads", "\n")
    cat("Reads on forward strand:", qc$frwrev[1], "\n")
    cat("Reads on reverse strand:", qc$frwrev[2], "\n")
    cat("Fraction of reads on forward strand (before QC):",
        qcmean(qc$frwrev), "\n")
}

## ----hist.score1, fig.width=8-------------------------------------------------
{
cat("Average alignment score, after filter:", qcmean(qc$hist.score1),    "\n")
cat("Average alignment score, no filter:   ", qcmean(qc$bf.hist.score1), "\n")
par(mfrow=c(1,2))
plot(qc$hist.score1)
plot(qc$bf.hist.score1)
}

## ----hist.length.matched, fig.width=8-----------------------------------------
{
cat("Average aligned length, after filter:", 
    qcmean(qc$hist.length.matched),    "\n")
cat("Average aligned length, no filter:   ",
    qcmean(qc$bf.hist.length.matched), "\n")
par(mfrow = c(1,2))
plot(qc$hist.length.matched)
plot(qc$bf.hist.length.matched)
}

## ----hist.edit.dist1----------------------------------------------------------
{
cat("Average edit distance, after filter:", 
    qcmean(qc$hist.edit.dist1),    "\n")
cat("Average edit distance, no filter:   ", 
    qcmean(qc$bf.hist.edit.dist1), "\n")
par(mfrow = c(1,2))
plot(qc$hist.edit.dist1)
plot(qc$bf.hist.edit.dist1)
}

## ----cnt.nonCpG.reads---------------------------------------------------------
{
cat("Non-CpG reads:", qc$cnt.nonCpG.reads[1], "\n")
cat("This is ", qcmean(qc$cnt.nonCpG.reads)*100, "% of recorded reads", sep="")
}

## ----avg.cpg.coverage---------------------------------------------------------
{
cat("Summed across", qc$nbams, "bams", "\n")
cat("Average     CpG coverage:", qc$avg.cpg.coverage, "\n")
cat("Average non-CpG coverage:", qc$avg.noncpg.coverage, "\n")
cat("Enrichment ratio:", qc$avg.cpg.coverage / qc$avg.noncpg.coverage, "\n")
cat("Noise level:", qc$avg.noncpg.coverage / qc$avg.cpg.coverage)
}

## ----avg.coverage.by.density--------------------------------------------------
{
cat("Highest coverage is observed at CpG density of",
    qcmean(qc$avg.coverage.by.density)^2)
plot(qc$avg.coverage.by.density)
}

## ----hist.isolated.dist1------------------------------------------------------
plot(qc$hist.isolated.dist1)

## ----chrXY--------------------------------------------------------------------
{
cat("ChrX reads: ", qc$chrX.count[1], ", which is ",
    qcmean(qc$chrX.count)*100, "% of total", sep="", "\n")
cat("ChrY reads: ", qc$chrY.count[1], ", which is ",
    qcmean(qc$chrY.count)*100, "% of total", sep="", "\n")
}