inst/doc/qcmetrics.R

In qcmetrics: A Framework for Quality Control

## ----style, eval=TRUE, echo=FALSE, results="asis"--------------------------
BiocStyle::latex()

## ----env, include=FALSE, echo=FALSE, cache=FALSE-----
library("knitr")
opts_chunk$set(fig.align = 'center', 
               fig.show = 'hold', 
               par = TRUE,
               prompt = FALSE,
               tidy = FALSE,
               eval = TRUE,
               stop_on_error = 1L,
               comment = "##")
options(replace.assign = TRUE, 
        width = 55)

suppressPackageStartupMessages(library("qcmetrics"))
suppressPackageStartupMessages(library("MAQCsubsetAFX"))
suppressPackageStartupMessages(library("yaqcaffy"))
suppressPackageStartupMessages(library("affy"))
suppressPackageStartupMessages(library("AnnotationDbi"))
suppressPackageStartupMessages(library("RforProteomics"))
suppressPackageStartupMessages(library("mzR"))
suppressPackageStartupMessages(library("MSnbase"))
set.seed(1)

## ----qcmetric----------------------------------------
library("qcmetrics")
qc <- QcMetric(name = "A test metric")
qcdata(qc, "x") <- rnorm(100)
qcdata(qc) ## all available qcdata
summary(qcdata(qc, "x")) ## get x
show(qc) ## or just qc
status(qc) <- TRUE
qc

## ----qcmetricplot, dev='pdf', fig.width = 4, fig.height = 4, tidy = FALSE----
plot(qc)
plot(qc) <-
    function(object, ... ) boxplot(qcdata(object, "x"), ...)
plot(qc)

## ----qcmetrics---------------------------------------
qcm <- QcMetrics(qcdata = list(qc))
qcm
metadata(qcm) <- list(author = "Prof. Who",
                      lab = "Big lab")
qcm
mdata(qcm)

## ----qcmdataupdate-----------------------------------
metadata(qcm) <- list(author = "Prof. Who",
                      lab = "Cabin lab",
                      University = "Universe-ity")
mdata(qcm)

## ----maqcdata, eval=FALSE----------------------------
#  library("MAQCsubsetAFX")
#  data(refA)
#  library("affy")
#  deg <- AffyRNAdeg(refA)
#  library("yaqcaffy")
#  yqc <- yaqc(refA)

## ----maqcdata0, echo=FALSE---------------------------
load(system.file("extdata/deg.rda", package = "qcmetrics"))
load(system.file("extdata/yqc.rda", package = "qcmetrics"))

## ----maqc1-------------------------------------------
qc1 <- QcMetric(name = "Affy RNA degradation slopes")
qcdata(qc1, "deg") <- deg
plot(qc1) <- function(object, ...) {
    x <- qcdata(object, "deg")
    nms <- x$sample.names
    plotAffyRNAdeg(x, col = 1:length(nms), ...)
    legend("topleft", nms, lty = 1, cex = 0.8, 
           col = 1:length(nms), bty = "n")
}
status(qc1) <- TRUE
qc1

## ----maqc2-------------------------------------------
qc2 <- QcMetric(name = "Affy RNA degradation ratios")
qcdata(qc2, "yqc") <- yqc
plot(qc2) <- function(object, ...) {
    par(mfrow = c(1, 2))
    yaqcaffy:::.plotQCRatios(qcdata(object, "yqc"), "all", ...)
}
status(qc2) <- FALSE
qc2

## ----maqcm-------------------------------------------
maqcm <- QcMetrics(qcdata = list(qc1, qc2))
maqcm

## ----maqcreport0, echo = FALSE, message = FALSE------
qcReport(maqcm, reportname = "rnadeg", clean = FALSE)

## ----maqcreport, eval = FALSE------------------------
#  qcReport(maqcm, reportname = "rnadeg", type = "pdf")

## ----maqcwrap, tidy=FALSE----------------------------
rnadeg

## ----qcwrap2, eval = FALSE---------------------------
#  maqcm <- rnadeg(refA)

## ----qcwrapstatus0, echo=FALSE-----------------------
status(maqcm) <- c(NA, NA)

## ----qcwrapstatus------------------------------------
status(maqcm)
## check the QC data 
(status(maqcm) <- c(TRUE, FALSE))

## ----qcwrap3, eval = FALSE---------------------------
#  maqcm <- rnadeg(refA, type = "pdf")

## ----protdata0, echo=FALSE---------------------------
load(system.file("extdata/exp.rda", package = "qcmetrics"))

## ----protdata, eval=FALSE----------------------------
#  library("RforProteomics")
#  msfile <- getPXD000001mzXML()
#  library("MSnbase")
#  exp <- readMSData(msfile, verbose = FALSE)

## ----protqc1, cache=TRUE, tidy=FALSE-----------------
qc1 <- QcMetric(name = "Chromatogram")
x <- rtime(exp)
y <- precursorIntensity(exp)
o <- order(x)
qcdata(qc1, "x") <- x[o]
qcdata(qc1, "y") <- y[o]
plot(qc1) <- function(object, ...)
    plot(qcdata(object, "x"),
         qcdata(object, "y"),
         col = "darkgrey", type ="l",
         xlab = "retention time",
         ylab = "precursor intensity")

## ----protqc2, cache=TRUE-----------------------------
qc2 <- QcMetric(name = "MS space")
qcdata(qc2, "p2d") <- plot2d(exp, z = "charge", plot = FALSE)
plot(qc2) <- function(object) { 
    require("ggplot2")
    print(qcdata(object, "p2d"))
}

## ----protqc3, cache=TRUE, messages=FALSE, tidy=FALSE, warnings=FALSE----
qc3 <- QcMetric(name = "m/z delta plot")
qcdata(qc3, "pmz") <- plotMzDelta(exp, plot = FALSE,
                                  verbose = FALSE)
plot(qc3) <- function(object) 
    suppressWarnings(print(qcdata(object, "pmz")))

## ----protqcm, tidy=FALSE-----------------------------
protqcm <- QcMetrics(qcdata = list(qc1, qc2, qc3))
metadata(protqcm) <- list(
    data = "PXD000001",
    instrument = experimentData(exp)@instrumentModel,
    source = experimentData(exp)@ionSource,
    analyser = experimentData(exp)@analyser,
    detector = experimentData(exp)@detectorType,
    manufacurer = experimentData(exp)@instrumentManufacturer)

## ----protreport0, echo = FALSE, message = FALSE------
qcReport(protqcm, reportname = "protqc", clean=FALSE, quiet=TRUE)

## ----protreport, eval=FALSE--------------------------
#  qcReport(protqcm, reportname = "protqc")

## ----n15ex-------------------------------------------
library("ggplot2")
library("MSnbase")
data(n15psm)
psm

## ----qcinc, tidy=FALSE-------------------------------
## incorporation rate QC metric
qcinc <- QcMetric(name = "15N incorporation rate")
qcdata(qcinc, "inc") <- fData(psm)$inc
qcdata(qcinc, "tr") <- 97.5
status(qcinc) <- median(qcdata(qcinc, "inc")) > qcdata(qcinc, "tr")

## ----qcinc2, tidy=FALSE------------------------------
show(qcinc) <- function(object) {
    qcshow(object, qcdata = FALSE) 
    cat(" QC threshold:", qcdata(object, "tr"), "\n")
    cat(" Incorporation rate\n")
    print(summary(qcdata(object, "inc")))
    invisible(NULL)
}

## ----qcinc3, tidy=FALSE------------------------------
plot(qcinc) <- function(object) {
    inc <- qcdata(object, "inc")
    tr <- qcdata(object, "tr")
    lab <- "Incorporation rate"
    dd <- data.frame(inc = qcdata(qcinc, "inc"))
    p <- ggplot(dd, aes(factor(""), inc)) +
        geom_jitter(colour = "#4582B370", size = 3) + 
    geom_boxplot(fill = "#FFFFFFD0", colour = "#000000",
                 outlier.size = 0) +
    geom_hline(yintercept = tr, colour = "red",
               linetype = "dotted", size = 1) +
    labs(x = "", y = "Incorporation rate") 
    p
}

## ----combinefeatures, tidy = FALSE-------------------
fData(psm)$modseq <- ## pep seq + PTM
    paste(fData(psm)$Peptide_Sequence, 
          fData(psm)$Variable_Modifications, sep = "+")
pep <- combineFeatures(psm,
                       as.character(fData(psm)$Peptide_Sequence), 
                       "median", verbose = FALSE)
modpep <- combineFeatures(psm,
                          fData(psm)$modseq,
                          "median", verbose = FALSE)
prot <- combineFeatures(psm,
                        as.character(fData(psm)$Protein_Accession), 
                        "median", verbose = FALSE)

## ----qclfc, tidy=FALSE-------------------------------
## calculate log fold-change
qclfc <- QcMetric(name = "Log2 fold-changes")
qcdata(qclfc, "lfc.psm") <-
    log2(exprs(psm)[,"unlabelled"] / exprs(psm)[, "N15"])
qcdata(qclfc, "lfc.pep") <-
    log2(exprs(pep)[,"unlabelled"] / exprs(pep)[, "N15"])
qcdata(qclfc, "lfc.modpep") <-
    log2(exprs(modpep)[,"unlabelled"] / exprs(modpep)[, "N15"])
qcdata(qclfc, "lfc.prot") <-
    log2(exprs(prot)[,"unlabelled"] / exprs(prot)[, "N15"])
qcdata(qclfc, "explfc") <- c(-0.5, 0.5)

status(qclfc) <-
    median(qcdata(qclfc, "lfc.psm")) > qcdata(qclfc, "explfc")[1] &
    median(qcdata(qclfc, "lfc.psm")) < qcdata(qclfc, "explfc")[2]

## ----qclfc2, tidy=FALSE------------------------------
show(qclfc) <- function(object) {
    qcshow(object, qcdata = FALSE) ## default
    cat(" QC thresholds:", qcdata(object, "explfc"), "\n")
    cat(" * PSM log2 fold-changes\n")
    print(summary(qcdata(object, "lfc.psm")))
    cat(" * Modified peptide log2 fold-changes\n")
    print(summary(qcdata(object, "lfc.modpep")))
    cat(" * Peptide log2 fold-changes\n")
    print(summary(qcdata(object, "lfc.pep")))
    cat(" * Protein log2 fold-changes\n")
    print(summary(qcdata(object, "lfc.prot")))    
    invisible(NULL)
}
plot(qclfc) <- function(object) {
    x <- qcdata(object, "explfc")
    plot(density(qcdata(object, "lfc.psm")),
         main = "", sub = "", col = "red",
         ylab = "", lwd = 2,
         xlab = expression(log[2]~fold-change))    
    lines(density(qcdata(object, "lfc.modpep")),
          col = "steelblue", lwd = 2)
    lines(density(qcdata(object, "lfc.pep")),
          col = "blue", lwd = 2)
    lines(density(qcdata(object, "lfc.prot")),
          col = "orange")
    abline(h = 0, col = "grey")
    abline(v = 0, lty = "dotted")    
    rect(x[1], -1, x[2], 1, col = "#EE000030",
         border = NA)
    abline(v = median(qcdata(object, "lfc.psm")),
           lty = "dashed", col = "blue")
    legend("topright",
           c("PSM", "Peptides", "Modified peptides", "Proteins"),
           col = c("red", "steelblue", "blue", "orange"), lwd = 2,
           bty = "n")
}

## ----qcnb, tidy=FALSE--------------------------------
## number of features
qcnb <- QcMetric(name = "Number of features")
qcdata(qcnb, "count") <- c(
    PSM = nrow(psm),
    ModPep = nrow(modpep),
    Pep = nrow(pep),
    Prot = nrow(prot))
qcdata(qcnb, "peptab") <-
    table(fData(psm)$Peptide_Sequence)
qcdata(qcnb, "modpeptab") <-
    table(fData(psm)$modseq)
qcdata(qcnb, "upep.per.prot") <- 
    fData(psm)$Number_Of_Unique_Peptides

## ----qcnb2, tidy=FALSE-------------------------------
show(qcnb) <- function(object) {
    qcshow(object, qcdata = FALSE)
    print(qcdata(object, "count"))
}
plot(qcnb) <- function(object) {
    par(mar = c(5, 4, 2, 1))
    layout(matrix(c(1, 2, 1, 3, 1, 4), ncol = 3))
    barplot(qcdata(object, "count"), horiz = TRUE, las = 2)
    barplot(table(qcdata(object, "modpeptab")),
            xlab = "Modified peptides")
    barplot(table(qcdata(object, "peptab")),
            xlab = "Peptides")
    barplot(table(qcdata(object, "upep.per.prot")),
            xlab = "Unique peptides per protein ")
}

## ----n15qcm, tidy=FALSE------------------------------
n15qcm <- QcMetrics(qcdata = list(qcinc, qclfc, qcnb))
qcReport(n15qcm, reportname = "n15qcreport",
         title = expinfo(experimentData(psm))["title"],
         author = expinfo(experimentData(psm))["contact"],
         clean = FALSE)

## ----Qc2Tex, tidy=FALSE------------------------------
qcmetrics:::Qc2Tex
qcmetrics:::Qc2Tex(maqcm, 1)

## ----Qc2Tex2, tidy=FALSE-----------------------------
Qc2Tex2

## ----maqcreport2, echo = FALSE, message = FALSE------
qcReport(maqcm, reportname = "rnadeg2", clean = FALSE, qcto = Qc2Tex2)

## ----maqcreport3, echo = FALSE, message = FALSE------
qcReport(maqcm, reportname = "rnadeg3", clean = FALSE, qcto = Qc2Tex3)

## ----maqcreport4, eval = FALSE-----------------------
#  qcReport(maqcm, reportname = "rnadeg2", qcto = Qc2Tex2)

## ----qcpkg0, eval=FALSE------------------------------
#  package.skeleton("RnaDegQC", list = "rnadeg")

## ----sessioninfo, results='asis', echo=FALSE---------
toLatex(sessionInfo())