proteoQC: an R package for proteomics data quality assessment"

knitr::opts_chunk$set(tidy = FALSE,message = FALSE)
library("BiocStyle")
BiocStyle::markdown()
suppressPackageStartupMessages(library("proteoQC"))
suppressPackageStartupMessages(library("R.utils"))

Introduction

The proteoQC package provides a integrated pipeline for mass spectrometry-based proteomics quality control. It allows to generate a dynamic report starting from a set of mgf or mz[X]ML format peak list files, a protein database file and a description file of the experimental design. It performs an MS/MS search against the protein data base using the X!Tandem search engine [@Craig:2004] and the rTANDEM package [@rTANDEM]. The results are then summarised and compiled into an interactive html report using the Nozzle.R1 package [@Nozzle.R1,@Gehlenborg:2013].

Example data

We are going to use parts a dataset from the ProteomeXchange repository (http://www.proteomexchange.org/). We will use the rpx package to accessed and downloaded the data.

library("rpx")
px <- PXDataset("PXD000864")
px

There are a total of r length(pxfiles(px)) files available from the ProteomeXchange repository, including raw data files (raw), result files (-pride.xml.gz), (compressed) peak list files (.mgf.gz) and, the fasta database file (TTE2010.zip) and one README.txt file.

head(pxfiles(px))
tail(pxfiles(px))

The files, in particular the mgf files that will be used in the rest of this document are named as follows TTE-CC-B-FR-R where CC takes values 55 or 75 and stands for the bacteria culture temperature in degree Celsius, B stands for the biological replicate (only 1 here), FR represents the fraction number from 01 to 12 and the leading R documents one of three technical replicates. (See also http://www.ebi.ac.uk/pride/archive/projects/PXD000864 for details). Here, we will make use of a limited number of samples below. First, we create a vector that stores the file names of interest.

mgfs <- grep("mgf", pxfiles(px), value = TRUE)
mgfs <- grep("-0[5-6]-[1|2]", mgfs, value=TRUE)
mgfs

These files can be downloaded [^1] using the pxget, providing the relevant data object (here px) and file names to be downloaded (see ?pxget for details). We also need to uncompress (using gunzip) the files.

mgffiles <- pxget(px, mgfs)
library("R.utils")
mgffiles <- sapply(mgffiles, gunzip)

To reduce the file size of the demonstration data included for this package, we have trimmed the peak lists to 1/10 of the original number of spectra. All the details are provided in the vignette source.

## Generate the lightweight qc report, 
## trim the mgf files to 1/10 of their size.

trimMgf <- function(f, m = 1/10, overwrite = FALSE) {
    message("Reading ", f)
    x <- readLines(f)
    beg <- grep("BEGIN IONS", x)
    end <- grep("END IONS", x)
    n <- length(beg)
    message("Sub-setting to ", m)
    i <- sort(sample(n, floor(n * m)))
    k <- unlist(mapply(seq, from = beg[i], to = end[i]))
    if (overwrite) {
        unlink(f)
        message("Writing ", f)
        writeLines(x[k], con = f)
        return(f)
    } else {
        g <- sub(".mgf", "_small.mgf", f)
        message("Writing ", g)
        writeLines(x[k], con = g)
        return(g)
    }    
}

set.seed(1)
mgffiles <- sapply(mgffiles, trimMgf, overwrite = TRUE)

Similarly, below we download the database file and unzip it.

fas <- pxget(px, "TTE2010.zip")
fas <- unzip(fas)
fas

Running proteoQC

## code to regenerate the design file
sample <- rep(c("55","75"),each=4)
techrep <- rep(1:2, 4)
biorep <- rep(1, length(mgffiles))
frac <- rep((rep(5:6, each = 2)), 2)
des <- data.frame(file = mgffiles,
                  sample = sample,
                  bioRep = biorep, techRep = techrep,
                  fraction = frac,
                  row.names = NULL)

write.table(des, sep = " ", row.names=FALSE,
            quote = FALSE,
            file = "../inst/extdata/PXD000864-design.txt")

Preparing the QC

The first step in the proteoQC pipeline is the definition of a design file, that provides the mgf file names, sample numbers, biological biocRep) and technical (techRep) replicates and fraction numbers in a simple space-separated tabular format. We provide such a design file for our r length(mgfs) files of interest.

design <- system.file("extdata/PXD000864-design.txt", package = "proteoQC")
design
read.table(design, header = TRUE)

Running the QC

We need to load the proteoQC package and call the msQCpipe function, providing appropriate input parameters, in particular the design file, the fasta protein database, the outdir output directory that will contain the final quality report and various other peptide spectrum matching parameters that will be passed to the rTANDEM package. See ?msQCpipe for a more in-depth description of all its arguments. Please note that if you take mz[X]ML format files as input, you must make sure that you have installed the rTANDEM that the version is greater than 1.5.1.

qcres <- msQCpipe(spectralist = design,
                  fasta = fas, 
                  outdir = "./qc",
                  miss  = 0,
                  enzyme = 1, varmod = 2, fixmod = 1,
                  tol = 10, itol = 0.6, cpu = 2,
                  mode = "identification")

The msQCpipe function will run each mgf input file documented in the design file and search it against the fasta database using the tandem function from the rTANDEM. This might take some time depending on the number of files to be searched and the search parameters. The code chunk above takes about 3 minutes using 2 cores (cpu = 2 above) on a modern laptop.

You can load the pre-computed quality control directory and result data that a shipped with proteoQC as shown below:

zpqc <- system.file("extdata/qc.zip", package = "proteoQC")
unzip(zpqc)
qcres <- loadmsQCres("./qc")
print(qcres)

Set MS/MS searching parameters

When we perform the QC analysis, we need to set several parameters for MS/MS searching. proteoQC provides a table about modifications. Users can select modifications using this table. Please use function showMods to print the available modifications. For the enzyme setting, please use function showEnzyme to print the available enzyme.

showMods()

Generating the QC report

The final quality report can be generated with the reportHTML, passing the qcres object produced by the msQCpipe function above or the directory storing the QC data, as defined as parameter to the msQCpipe.

html <- reportHTML(qcres)

or

html <- reportHTML("./qc")
## Remove these files as they are really big
## but this breaks reportHTML(qcres), though
unlink("./qc/database/target_decoy.fasta")
unlink("./qc/result/*_xtandem.xml")
unlink("../inst/extdata/qc.zip")
zip("../inst/extdata/qc.zip", "./qc")

The report can then be opened by opening the qc/qc_report.html file in a web browser or directly with browseURL(html).

The QC report

The dynamic html report is composed of 3 sections: an introduction, a methods and data section and a result part. The former are purely descriptive and summarise the design matrix and analysis parameters, as passed to msQCpipe.

The respective sections and sub-sections can be expanded and collapsed and each figure in the report can be zoomed in. While the dynamic html report is most useful for interactive inspection, it is also possible to print the complete report for archiving.

The results section provides tables and graphics that summarise

Some useful functions

Protein inference

Protein inference from peptide identifications in shotgun proteomics is a very important task. We provide a function proteinGroup for this purpose. This function is based on the method used in our another package sapFinder [@wen2014sapfinder]. You can use the function as below:

pep.zip <- system.file("extdata/pep.zip", package = "proteoQC")
unzip(pep.zip)
proteinGroup(file = "pep.txt", outfile = "pg.txt")

Isobaric tagging reagent labeling efficiency

The labeling efficiency of the isobaric tag reagents to peptides, such as iTRAQ and TMT, is a very important experiment quality metrics. We provide a function labelRatio to calculate this metrics. You can use the function as below:

mgf.zip <- system.file("extdata/mgf.zip", package = "proteoQC")
unzip(mgf.zip)
a <- labelRatio("test.mgf",reporter = 2)

Precusor charge distribution

Given an MGF file, chargeStat function can be used to get the precusor charge distribution.

library(dplyr)
library(plotly)
mgf.zip <- system.file("extdata/mgf.zip", package = "proteoQC")
unzip(mgf.zip)
charge <- chargeStat("test.mgf")
pp <- plot_ly(charge, labels = ~Charge, values = ~Number, type = 'pie') %>%
        layout(title = 'Charge distribution',
        xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
        yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))
pp

Session information

All software and respective versions used to produce this document are listed below.

sessionInfo()

References

[^1]: In the interest of time, the files are not downloaded when this vignette is compiled and the quality metrics are pre-computed (see details below). These following code chunks can nevertheless be executed to reproduce the complete pipeline.}



Try the proteoQC package in your browser

Any scripts or data that you put into this service are public.

proteoQC documentation built on Oct. 31, 2019, 3:20 a.m.