QCs {ramwas} | R Documentation

## Quality Control Measures

### Description

RaMWAS calculates a number of QC measures for each BAM and saves them in R .rds files.

For full description of the QC measures and the ploting options run
`vignette("RW3_BAM_QCs")`

### Usage

    qcmean(x)
    ## S3 method for class 'NULL'
    qcmean(x)
    ## S3 method for class 'qcChrX'
    qcmean(x)
    ## S3 method for class 'qcChrY'
    qcmean(x)
    ## S3 method for class 'qcCoverageByDensity'
    qcmean(x)
    ## S3 method for class 'qcEditDist'
    qcmean(x)
    ## S3 method for class 'qcEditDistBF'
    qcmean(x)
    ## S3 method for class 'qcFrwrev'
    qcmean(x)
    ## S3 method for class 'qcHistScore'
    qcmean(x)
    ## S3 method for class 'qcHistScoreBF'
    qcmean(x)
    ## S3 method for class 'qcIsoDist'
    qcmean(x)
    ## S3 method for class 'qcLengthMatched'
    qcmean(x)
    ## S3 method for class 'qcLengthMatchedBF'
    qcmean(x)
    ## S3 method for class 'qcNonCpGreads'
    qcmean(x)

    ## S3 method for class 'qcHistScore'
    plot(x, samplename="", xstep = 25, ...)
    ## S3 method for class 'qcHistScoreBF'
    plot(x, samplename="", xstep = 25, ...)
    ## S3 method for class 'qcEditDist'
    plot(x, samplename="", xstep = 5, ...)
    ## S3 method for class 'qcEditDistBF'
    plot(x, samplename="", xstep = 5, ...)
    ## S3 method for class 'qcLengthMatched'
    plot(x, samplename="", xstep = 25, ...)
    ## S3 method for class 'qcLengthMatchedBF'
    plot(x, samplename="", xstep = 25, ...)
    ## S3 method for class 'qcIsoDist'
    plot(x, samplename="", xstep = 25, ...)
    ## S3 method for class 'qcCoverageByDensity'
    plot(x, samplename="", ...)

### Arguments

`x` | The QC object. See the examples below.
`samplename` | Name of the sample for plot title.
`xstep` | The distance between x axis ticks.
`...` | Parameters passed to the underlying `plot` or `barplot` function.

### Value

Function `qcmean` returns one value summary of most QC measures. Run `vignette("RW3_BAM_QCs")` for description of values returned by it.

### Author(s)

Andrey A Shabalin andrey.shabalin@gmail.com

### See Also

See vignettes: `browseVignettes("ramwas")`.

### Examples

    # Load QC data from a sample project
    filename = system.file("extdata", "bigQC.rds", package = "ramwas")
    qc = readRDS(filename)$qc

    ## The number of BAM files
    cat("N BAMs:", qc$nbams)

    ## Total number of reads in the BAM file(s)
    cat("Reads total:", qc$reads.total)

    ## Number of reads aligned to the reference genome
    cat("Reads aligned:", qc$reads.aligned, "\n")
    cat("This is ", qc$reads.aligned / qc$reads.total * 100,
        "% of all reads", sep="")

    ## Number of reads that passed minimum score filter and are recorded
    cat("Reads recorded:",qc$reads.recorded,"\n")
    cat("This is ", qc$reads.recorded / qc$reads.aligned * 100,
        "% of aligned reads", sep="")

    ## Number of recorded reads aligned to
    ## the forward and reverse strands respectively
    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:", qcmean(qc$frwrev), "\n")

    ## Distribution of the read scores
    cat("Average alignment score:", 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)

    ## Distribution of the length of the aligned part of the reads
    cat("Average aligned length:", 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)

    ## Distribution of edit distance between
    ## the aligned part of the read and the reference genome
    cat("Average edit distance:", 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)

    ## Number of reads after removal of duplicate reads
    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="")
    cat("Fraction of reads on forward strand (with    duplicates):",
        qcmean(qc$frwrev), "\n")
    cat("Fraction of reads on forward strand (without duplicates):",
        qcmean(qc$frwrev.no.repeats), "\n")

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

    ## Average coverage of CpGs and non-CpGs
    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)

    ## Coverage around isolated CpGs
    plot(qc$hist.isolated.dist1)

    ## Fraction of reads from chrX and chrY
    cat("ChrX reads: ", qc$chrX.count[1], ", which is ",
        qcmean(qc$chrX.count)*100, "% of total", sep="", "\n")
    cat("ChrX reads: ", qc$chrY.count[1], ", which is ",
        qcmean(qc$chrY.count)*100, "% of total", sep="", "\n")

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

[Package ramwas version 1.8.0 Index]