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