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