In lgatto/QFeaturesScpWorkshop2021: Reproducing a single-cell proteomics data analysis using QFeatures and scp

library(BiocStyle)
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      collapse = TRUE)

Last modified: r file.info("v01-QFeaturesClass.Rmd")$mtime
Compiled: r date()

Introduction

QFeatures is a package for the manipulation of mass spectrometry (MS)-based quantitative proteomics data. The quantitative data is generally obtained after running a preprocessing algorithm on raw MS profiles called spectra. The algorithm attempts to match each spectra to a peptide database to infer the peptide sequence. A spectrum for which a corresponding peptide could be found is called peptide to spectrum match (PSM). The intensities recorded by the MS are then used to quantify the PSMs. The quantified PSMs are used as input by QFeatures and undergo various data wrangling steps to reconstruct the peptide data and then the protein data, usually of interest for elucidating the biological question at hand.

MS-based quantification data can be represented as a matrix of quantitative values for features (PSMs, peptides, proteins) arranged along the rows, measured for a set of samples, arranged along the columns. The matrix format is a common representation for any quantitative data set. We will be using the SummarizedExperiment [@SE] class:

knitr::include_graphics("figs/SE.png")

• The sample (columns) metadata can be access with the colData() function.
• The features (rows) metadata can be access with the rowData() column.
• The quantiative data can be accessed with assay().
• assays() returns a list of matrix-like assays.

Note that other classes that inherits from SummarizedExperiment can also be used. We will see in the other vignette that for single-cell proteomics data, we will use SingleCellExperiment objects.

QFeatures

While mass spectrometers acquire data for spectra/peptides, the biological entity of interest are the protein. As part of the data processing, we are thus required to aggregate low-level quantitative features into higher level data.

par(mar = c(0, 0, 0, 0))
plot(NA, xlim = c(0, 12), ylim = c(0, 20),
     xaxt = "n", yaxt = "n",
     xlab = "", ylab = "", bty = "n")

for (i in 0:7)
    rect(0, i, 3, i+1, col = "lightgrey", border = "white")


for (i in 8:12)
    rect(0, i, 3, i+1, col = "steelblue", border = "white")

for (i in 13:18)
    rect(0, i, 3, i+1, col = "orange", border = "white")

for (i in 19)
    rect(0, i, 3, i+1, col = "darkgrey", border = "white")


for (i in 5:7)
    rect(5, i, 8, i+1, col = "lightgrey", border = "white")

for (i in 8:10)
    rect(5, i, 8, i+1, col = "steelblue", border = "white")

for (i in 11:13)
    rect(5, i, 8, i+1, col = "orange", border = "white")

for (i in 14)
    rect(5, i, 8, i+1, col = "darkgrey", border = "white")

rect(9, 8, 12, 8+1, col = "lightgrey", border = "white")
rect(9, 9, 12, 9+1, col = "steelblue", border = "white")
rect(9, 10, 12, 10+1, col = "orange", border = "white")
rect(9, 11, 12, 11+1, col = "darkgrey", border = "white")

segments(3, 8, 5, 8, lty = "dashed")
segments(3, 6, 5, 7, lty = "dashed")
segments(3, 4, 5, 6, lty = "dashed")
segments(3, 0, 5, 5, lty = "dashed")

segments(3, 10, 5, 9, lty = "dashed")
segments(3, 11, 5, 10, lty = "dashed")
segments(3, 13, 5, 11, lty = "dashed")

segments(3, 14, 5, 12, lty = "dashed")
segments(3, 16, 5, 13, lty = "dashed")
segments(3, 19, 5, 14, lty = "dashed")

segments(3, 20, 5, 15, lty = "dashed")


segments(8, 5, 9, 8, lty = "dashed")
segments(8, 8, 9, 9, lty = "dashed")
segments(8, 11, 9, 10, lty = "dashed")
segments(8, 14, 9, 11, lty = "dashed")
segments(8, 15, 9, 12, lty = "dashed")

We are going to start to familiarise ourselves with the QFeatures class implemented in the QFeatures package. The class is derived from the Bioconductor MultiAssayExperiment [@MAE] class. Let's start by loading the QFeatures package.

library("QFeatures")

Next, we load the feat1 test data, which is composed of single assay of class SummarizedExperiment composed of 10 rows and 2 columns.

data(feat1)
feat1

Let's perform some simple operations to familiarise ourselves with the QFeatures class:

• Extract the sample metadata using the colData() accessor (like you have previously done with SummarizedExperiment objects).

colData(feat1)

• Extract the first (and only) assay composing this QFeaures data using the [[ operator (as you have done to extract elements of a list) by using the assay's index or name.

feat1[[1]]
feat1[["psms"]]

• Extract the psms assay's row data and quantitative values.

assay(feat1[[1]])
rowData(feat1[[1]])

Feature aggregation

The central functionality of the QFeatures infrastructure is the aggregation of features into higher-level features while retaining the link between the different levels. This can be done with the aggregateFeatures() function.

The call below will

• operate on the psms assay of the feat1 objects;
• aggregate the rows the assay following the grouping defined in the peptides row data variables;
• perform aggregation using the colMeans() function;
• create a new assay named peptides and add it to the feat1 object.

feat1 <- aggregateFeatures(feat1, i = "psms",
                           fcol = "Sequence",
                           name = "peptides",
                           fun = colMeans)
feat1

Note: we do not claim that colMeans is the best aggregation method. There are other more robust methods for aggregating features, for instance MsCoreUtils::robustSummary is a good alternative. See the documentation of aggregateFeatures for more information.

Exercise at home: check that you understand the effect of feature aggregation and repeat the calculations manually. Compare your manual results agains the content of the new assay's row data.

Solution

## Data before aggregation
data.frame(assay(feat1[["psms"]]), 
           Sequence = rowData(feat1[["psms"]])$Sequence)

## SYGFNAAR
colMeans(assay(feat1[["psms"]])[1:3, ])
assay(feat1[[2]])["SYGFNAAR", ]

## ELGNDAYK
colMeans(assay(feat1[["psms"]])[4:6, ])
assay(feat1[[2]])["ELGNDAYK", ]

## IAEESNFPFIK
colMeans(assay(feat1[["psms"]])[7:10, ])
assay(feat1[[2]])["IAEESNFPFIK", ]

rowData(feat1[[2]])

We can now aggregate the peptide-level data into a new protein-level assay using the colMedians() aggregation function.

feat1 <- aggregateFeatures(feat1, i = "peptides",
                           fcol = "Protein",
                           name = "proteins",
                           fun = colMedians)
feat1
assay(feat1[["proteins"]])

You can also get a graphical overview of the QFeatures object thanks to the plot function. Arrows represent the links between the assays.

plot(feat1)

Subsetting and filtering

The subsetting a QFeatures object can be performed using the 3-index subset operator, [. The first index will subset the features of interest. If the feature is linked to features from other assays, those features will also

knitr::include_graphics("figs/QFeatures_subsetting.png")

We provide an example of subsetting features that are associated to features from other assays. When we subset for protein A as shown below or using the subsetByFeature() function, this creates a new instance of class QFeatures containing assays with the expression data for protein, its peptides and their PSMs.

feat1["ProtA", , ]

The filterFeatures() function can be used to filter rows the assays composing a QFeatures object using the row data variables. We can for example retain rows that have a pval < 0.05, which would only keep rows in the psms assay because the pval is only relevant for that assay.

filterFeatures(feat1, ~ pval < 0.05)

On the other hand, if we filter assay rows for those that localise to the mitochondrion, we retain the relevant protein, peptides and PSMs.

filterFeatures(feat1, ~ location == "Mitochondrion")

As an exercise, let's filter rows that do not localise to the mitochondrion.

filterFeatures(feat1, ~ location != "Mitochondrion")

Another useful filtering functionality is filterNA that removes features based on the amount of missing data it contains. To illustrate this, we will load a data set containing missing data

data("ft_na")
ft_na
assay(ft_na[[1]])

We can remove feature that contain more that a given percentage missing values. For instance, let's remove features with more than 25 \% missing data. We expect to keep only 1 features as all the other features contain 33 \% missing values.

filterNA(ft_na, i = "na", pNA = 0.25)

The resulting QFeatures object contains a single feature as we expected.

Data processing

The Qfeatures package provides a few utility functions to process the data:

• Logarithmic transform: logTransform creates a new assay in a QFeatures object that contains the log-transformed quantification of the target assay.

logTransform(feat1, i = "proteins", base = 2, name = "logproteins")

• Normalization: normalize creates a new assay in a QFeatures object that contains the normalized quantification of the target assay. See the ?QFeatures::normalize documentation for a list of available normalization methods. If the method you are looking for is not available, you could also consider QFeatures::sweep.

normalize(feat1, i = "proteins", method = "center.median", name = "normproteins")

• Imputation: imputes predicts missing values in a new assay. See the ?MsCoreUtils::impute_matrix() for a list of available imputation methods. In this example, we impute missing values with zero.

assay(ft_na[[1]])
ft_na <- impute(ft_na, i = 1, method = "zero")
assay(ft_na[[1]])

Further reading

You can refer to the Quantitative features for mass spectrometry data vignette and the QFeature manual page for more details about the class.

Session information {-}

sessionInfo()

References {-}

lgatto/QFeaturesScpWorkshop2021 documentation built on Dec. 21, 2021, 10:44 a.m.