Processing quantitative proteomics data with QFeatures
1 November 2022
Abstract
This vignette describes how to process quantitative mass spectrometry data with QFeatures: cleaning up unneeded feature variables, adding an experimental design, filtering out contaminants and reverse hits, managing missing values, log-transforming, normalising and aggregating data. This vignette is distributed under a CC BY-SA license.Package
QFeatures 1.8.0
1 Reading data as QFeatures
We are going to use a subset of the CPTAC study 6 containing
conditions A and B (Paulovich et al. 2010). The peptide-level data, as
processed by MaxQuant (Cox and Mann 2008) is available in the msdata
package:
basename(f <- msdata::quant(pattern = "cptac", full.names = TRUE))## [1] "cptac_a_b_peptides.txt"From the names of the columns, we see that the quantitative columns,
starting with "Intensity." (note the dot!) are at positions 56 to
61.
names(read.delim(f))## [1] "Sequence" "N.term.cleavage.window"
## [3] "C.term.cleavage.window" "Amino.acid.before"
## [5] "First.amino.acid" "Second.amino.acid"
## [7] "Second.last.amino.acid" "Last.amino.acid"
## [9] "Amino.acid.after" "A.Count"
## [11] "R.Count" "N.Count"
## [13] "D.Count" "C.Count"
## [15] "Q.Count" "E.Count"
## [17] "G.Count" "H.Count"
## [19] "I.Count" "L.Count"
## [21] "K.Count" "M.Count"
## [23] "F.Count" "P.Count"
## [25] "S.Count" "T.Count"
## [27] "W.Count" "Y.Count"
## [29] "V.Count" "U.Count"
## [31] "Length" "Missed.cleavages"
## [33] "Mass" "Proteins"
## [35] "Leading.razor.protein" "Start.position"
## [37] "End.position" "Unique..Groups."
## [39] "Unique..Proteins." "Charges"
## [41] "PEP" "Score"
## [43] "Identification.type.6A_7" "Identification.type.6A_8"
## [45] "Identification.type.6A_9" "Identification.type.6B_7"
## [47] "Identification.type.6B_8" "Identification.type.6B_9"
## [49] "Experiment.6A_7" "Experiment.6A_8"
## [51] "Experiment.6A_9" "Experiment.6B_7"
## [53] "Experiment.6B_8" "Experiment.6B_9"
## [55] "Intensity" "Intensity.6A_7"
## [57] "Intensity.6A_8" "Intensity.6A_9"
## [59] "Intensity.6B_7" "Intensity.6B_8"
## [61] "Intensity.6B_9" "Reverse"
## [63] "Potential.contaminant" "id"
## [65] "Protein.group.IDs" "Mod..peptide.IDs"
## [67] "Evidence.IDs" "MS.MS.IDs"
## [69] "Best.MS.MS" "Oxidation..M..site.IDs"
## [71] "MS.MS.Count"(i <- grep("Intensity\\.", names(read.delim(f))))## [1] 56 57 58 59 60 61We now read these data using the readQFeatures function. The peptide
level expression data will be imported into R as an instance of class
QFeatures named cptac with an assay named peptides. We also use
the fnames argument to set the row-names of the peptides assay to
the peptide sequences.
library("QFeatures")
cptac <- readQFeatures(f, ecol = i, sep = "\t", name = "peptides", fnames = "Sequence")2 Encoding the experimental design
Below we update the sample (column) annotations to encode the two groups, 6A and 6B, and the original sample numbers.
cptac$group <- rep(c("6A", "6B"), each = 3)
cptac$sample <- rep(7:9, 2)
colData(cptac)## DataFrame with 6 rows and 2 columns
## group sample
## <character> <integer>
## Intensity.6A_7 6A 7
## Intensity.6A_8 6A 8
## Intensity.6A_9 6A 9
## Intensity.6B_7 6B 7
## Intensity.6B_8 6B 8
## Intensity.6B_9 6B 93 Filtering out contaminants and reverse hits
filterFeatures(cptac, ~ Reverse == "")## 'Reverse' found in 1 out of 1 assay(s)## An instance of class QFeatures containing 1 assays:
## [1] peptides: SummarizedExperiment with 11436 rows and 6 columnsfilterFeatures(cptac, ~ Potential.contaminant == "")## 'Potential.contaminant' found in 1 out of 1 assay(s)## An instance of class QFeatures containing 1 assays:
## [1] peptides: SummarizedExperiment with 11385 rows and 6 columnslibrary("magrittr")
cptac <- cptac %>%
filterFeatures(~ Reverse == "") %>%
filterFeatures(~ Potential.contaminant == "")## 'Reverse' found in 1 out of 1 assay(s)## 'Potential.contaminant' found in 1 out of 1 assay(s)4 Removing up unneeded feature variables
The spreadsheet that was read above contained numerous variables that are returned by MaxQuant, but not necessarily necessary in the frame of a downstream statistical analysis.
rowDataNames(cptac)## CharacterList of length 1
## [["peptides"]] Sequence N.term.cleavage.window ... MS.MS.CountThe only ones that we will be needing below are the peptides sequences
and the protein identifiers. Below, we store these variables of
interest and filter them using the selectRowData function.
rowvars <- c("Sequence", "Proteins", "Leading.razor.protein")
cptac <- selectRowData(cptac, rowvars)
rowDataNames(cptac)## CharacterList of length 1
## [["peptides"]] Sequence Proteins Leading.razor.protein5 Managing missing values
Missing values can be very numerous in certain proteomics experiments
and need to be dealt with carefully. The first step is to assess their
presence across samples and features. But before being able to do so,
we need to replace 0 by NA, given that MaxQuant encodes missing data
with a 0 using the zeroIsNA function.
cptac <- zeroIsNA(cptac, i = seq_along(cptac))
nNA(cptac, i = seq_along(cptac))## $nNA
## DataFrame with 1 row and 3 columns
## assay nNA pNA
## <character> <integer> <numeric>
## 1 peptides 30609 44.9194
##
## $nNArows
## DataFrame with 11357 rows and 4 columns
## assay name nNA pNA
## <character> <character> <integer> <numeric>
## 1 peptides AAAAGAGGAG... 4 66.6667
## 2 peptides AAAALAGGK 0 0.0000
## 3 peptides AAAALAGGKK 0 0.0000
## 4 peptides AAADALSDLE... 0 0.0000
## 5 peptides AAADALSDLE... 0 0.0000
## ... ... ... ... ...
## 11353 peptides YYSIYDLGNN... 6 100.0000
## 11354 peptides YYTFNGPNYN... 3 50.0000
## 11355 peptides YYTITEVATR 4 66.6667
## 11356 peptides YYTVFDRDNN... 6 100.0000
## 11357 peptides YYTVFDRDNN... 6 100.0000
##
## $nNAcols
## DataFrame with 6 rows and 4 columns
## assay name nNA pNA
## <character> <character> <integer> <numeric>
## 1 peptides Intensity.... 4669 41.1112
## 2 peptides Intensity.... 5388 47.4421
## 3 peptides Intensity.... 5224 45.9981
## 4 peptides Intensity.... 4651 40.9527
## 5 peptides Intensity.... 5470 48.1641
## 6 peptides Intensity.... 5207 45.8484The output of the nNA function tells us that
- there are currently close to 50% is missing values in the data;
- there are 4051 peptides with 0 missing values, 989 with a single missing values, … and 3014 peptides composed of only missing values;
- the range of missing values in the 6 samples is comparable and ranges between 4651 and 5470.
In this dataset, we have such a high number of peptides without any
data because the 6 samples are a subset of a larger dataset, and these
peptides happened to be absent in groups A and B. Below, we use
filterNA to remove all the peptides that contain one or more missing
values by using pNA = 0 (which also is the default value).
cptac <- filterNA(cptac, i = seq_along(cptac), pNA = 0)
cptac## An instance of class QFeatures containing 1 assays:
## [1] peptides: SummarizedExperiment with 4051 rows and 6 columnsI we wanted to keep peptides that have up to 90% of missing values,
corresponsing in this case to those that have only one value (i.e 5/6
percent of missing values), we could have set pNA to 0.9.
6 Counting unique features
Counting the number of unique features across samples can be used for
quality control or for assessing the identification efficiency between
different conditions or experimental set-ups. countUniqueFeatures
can be used to count the number of features that are contained in each
sample of an assay from a QFeatures object. For instance, we can
count the number of (non-missing) peptides per sample from the
peptides assay. Note that the counts are automatically stored in the
colData of cptac, under peptide_counts:
cptac <- countUniqueFeatures(cptac,
i = "peptides",
colDataName = "peptide_counts")
colData(cptac)## DataFrame with 6 rows and 3 columns
## group sample peptide_counts
## <character> <integer> <integer>
## Intensity.6A_7 6A 7 4051
## Intensity.6A_8 6A 8 4051
## Intensity.6A_9 6A 9 4051
## Intensity.6B_7 6B 7 4051
## Intensity.6B_8 6B 8 4051
## Intensity.6B_9 6B 9 4051We can also count the number of unique proteins. We therefore need to
tell countUniqueFeatures that we need to group by protein (the
protein name is stored in the rowData under Proteins):
cptac <- countUniqueFeatures(cptac,
i = "peptides",
groupBy = "Proteins",
colDataName = "protein_counts")
colData(cptac)## DataFrame with 6 rows and 4 columns
## group sample peptide_counts protein_counts
## <character> <integer> <integer> <integer>
## Intensity.6A_7 6A 7 4051 1125
## Intensity.6A_8 6A 8 4051 1125
## Intensity.6A_9 6A 9 4051 1125
## Intensity.6B_7 6B 7 4051 1125
## Intensity.6B_8 6B 8 4051 1125
## Intensity.6B_9 6B 9 4051 11257 Imputation
The impute method can be used to perform missing value imputation
using a variety of imputation methods. The method takes an instance of
class QFeatures (or a SummarizedExperiment) as input, an a
character naming the desired method (see ?impute for the complete
list with details) and returns a new instance of class QFeatures (or
SummarizedExperiment) with imputed data.
As described in more details in (Lazar et al. 2016), there are two types of mechanisms resulting in missing values in LC/MSMS experiments.
Missing values resulting from absence of detection of a feature, despite ions being present at detectable concentrations. For example in the case of ion suppression or as a result from the stochastic, data-dependent nature of the MS acquisition method. These missing value are expected to be randomly distributed in the data and are defined as missing at random (MAR) or missing completely at random (MCAR).
Biologically relevant missing values, resulting from the absence of the low abundance of ions (below the limit of detection of the instrument). These missing values are not expected to be randomly distributed in the data and are defined as missing not at random (MNAR).
MAR and MCAR values can be reasonably well tackled by many imputation methods. MNAR data, however, requires some knowledge about the underlying mechanism that generates the missing data, to be able to attempt data imputation. MNAR features should ideally be imputed with a left-censor (for example using a deterministic or probabilistic minimum value) method. Conversely, it is recommended to use hot deck methods (for example nearest neighbour, maximum likelihood, etc) when data are missing at random.
Figure 1: Mixed imputation method
Black cells represent presence of quantitation values and light grey corresponds to missing data. The two groups of interest are depicted in green and blue along the heatmap columns. Two classes of proteins are annotated on the left: yellow are proteins with randomly occurring missing values (if any) while proteins in brown are candidates for non-random missing value imputation.
It is anticipated that the identification of both classes of missing values will depend on various factors, such as feature intensities and experimental design. Below, we use perform mixed imputation, applying nearest neighbour imputation on the 654 features that are assumed to contain randomly distributed missing values (if any) (yellow on figure 1) and a deterministic minimum value imputation on the 35 proteins that display a non-random pattern of missing values (brown on figure 1).
8 Data transformation
When analysing continuous data using parametric methods (such as t-test or linear models), it is often necessary to log-transform the data. The figure below (left) show that how our data is mainly composed of small values with a long tail of larger ones, which is a typical pattern of quantitative omics data.
Below, we use the logTransform function to log2-transform our
data. This time, instead of overwriting the peptides assay, we are
going to create a new one to contain the log2-transformed data.
cptac <- addAssay(cptac,
logTransform(cptac[[1]]),
name = "peptides_log")
cptac## An instance of class QFeatures containing 2 assays:
## [1] peptides: SummarizedExperiment with 4051 rows and 6 columns
## [2] peptides_log: SummarizedExperiment with 4051 rows and 6 columnsThe addAssay() function is the general function that adds new assays
to a QFeatures object. The step above could also be fun with the
following syntax, that implicitly returns an updated QFeatures
object.
logTransform(cptac,
i = "peptides",
name = "log_peptides")par(mfrow = c(1, 2))
limma::plotDensities(assay(cptac[[1]]))
limma::plotDensities(assay(cptac[[2]]))
Figure 2: Quantitative data in its original scale (left) and log2-transformed (right)
9 Normalisation
Assays in QFeatures objects can be normalised with the normalize
function. The type of normalisation is defined by the method
argument; below, we use quantile normalisation, store the normalised
data into a new experiment, and visualise the resulting data.
cptac <- addAssay(cptac,
normalize(cptac[["peptides_log"]], method = "center.median"),
name = "peptides_norm")
cptac## An instance of class QFeatures containing 3 assays:
## [1] peptides: SummarizedExperiment with 4051 rows and 6 columns
## [2] peptides_log: SummarizedExperiment with 4051 rows and 6 columns
## [3] peptides_norm: SummarizedExperiment with 4051 rows and 6 columnsAs above, the normalize() function can also be firectly applied to
the QFeatures object.
normalize(cptac,
i = "log_peptides",
name = "lognorm_peptides",
method = "center.median")par(mfrow = c(1, 2))
limma::plotDensities(assay(cptac[["peptides_log"]]))
limma::plotDensities(assay(cptac[["peptides_norm"]]))