## ----env, include=FALSE, echo=FALSE, cache=FALSE------------------------------
suppressPackageStartupMessages(library("qcmetrics"))
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-------------------------------------------------------------
plot(qc)
plot(qc) <- function(object, ... ) boxplot(qcdata(object, "x"), ...)
plot(qc)

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

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

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

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

## ----protqc1------------------------------------------------------------------
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----------------------------------------------------
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------------------------------------------------------------------
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, eval = FALSE-----------------
#  qcReport(protqcm, reportname = "protqc", clean=FALSE, quiet=TRUE)

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

## ---- results='markup', out.width="100%", fig.align="center", fig.cap="Proteomics QC report", echo=FALSE----
knitr::include_graphics("./figs/protqc.png")

## ---- eval = FALSE------------------------------------------------------------
#  browseURL(example_reports("protqc"))

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

## ----qcinc--------------------------------------------------------------------
## 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-------------------------------------------------------------------
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-------------------------------------------------------------------
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----------------------------------------------------------
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--------------------------------------------------------------------
## 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-------------------------------------------------------------------
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---------------------------------------------------------------------
## 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--------------------------------------------------------------------
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-------------------------------------------------------------------
n15qcm <- QcMetrics(qcdata = list(qcinc, qclfc, qcnb))

## ----n15report, eval = FALSE--------------------------------------------------
#  qcReport(n15qcm, reportname = "n15qcreport",
#           title = expinfo(experimentData(psm))["title"],
#           author = expinfo(experimentData(psm))["contact"],
#           clean = FALSE)

## ---- results='markup', out.width="100%", fig.align="center", fig.cap="N15 QC report", echo=FALSE----
knitr::include_graphics("./figs/n15qcreport.png")

## ---- eval = FALSE------------------------------------------------------------
#  browseURL(example_reports("n15qc"))

## ----Qc2Tex-------------------------------------------------------------------
qcmetrics:::Qc2Tex
qcmetrics:::Qc2Tex(n15qcm, 1)

## ----maqcreport4, eval = FALSE------------------------------------------------
#  qcReport(n15qcm, reportname = "report", qcto = Qc2Tex2)
#  qcReport(n15qcm, reportname = "report", qcto = Qc2Tex3) ## for smiley/frowney

## ---- results='markup', out.width="100%", fig.cap="Customised QC report", echo=FALSE----
knitr::include_graphics("./figs/custom.png")

## ---- eval = FALSE------------------------------------------------------------
#  browseURL(example_reports("custom"))

## ----qcpkg0, eval=FALSE-------------------------------------------------------
#  package.skeleton("N15QC", list = "n15qc")

## ----si-----------------------------------------------------------------------
sessionInfo()