In wenbostar/proteoQC: An R package for proteomics data quality control
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
- Summaries of identification results for individual files as well
as technical and biological replicates at the protein, peptide and
spectrum levels.
- Summary overview charts that describe number of missed
cleavages, peptide charge distributions, peptide length, precursor
and fragment ion mass deviations, number of unique spectra/peptides
per proteins and protein mass distributions for each sample.
- A contamination summary table generated using the common
Repository of Adventitious Proteins (\textit{cRAP}).
- Reproducibility summaries that compare fractions, replicates and
samples, representing total number of spectra, number of identified
spectra, number of peptides and proteins and overlap of peptides and
proteins across replicates.
- Summary histograms of mass accuracies for fragment and precursor
ions.
- A summary of the separation efficiency showing the effect of
accumulating fractions for all samples.
- A summary of identification-independent QC metrics.
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.}
wenbostar/proteoQC documentation built on Nov. 20, 2020, 10:15 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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
- Summaries of identification results for individual files as well as technical and biological replicates at the protein, peptide and spectrum levels.
- Summary overview charts that describe number of missed cleavages, peptide charge distributions, peptide length, precursor and fragment ion mass deviations, number of unique spectra/peptides per proteins and protein mass distributions for each sample.
- A contamination summary table generated using the common Repository of Adventitious Proteins (\textit{cRAP}).
- Reproducibility summaries that compare fractions, replicates and samples, representing total number of spectra, number of identified spectra, number of peptides and proteins and overlap of peptides and proteins across replicates.
- Summary histograms of mass accuracies for fragment and precursor ions.
- A summary of the separation efficiency showing the effect of accumulating fractions for all samples.
- A summary of identification-independent QC metrics.
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.}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.