A combined test for feature selection on sparse metaproteomics data—an alternative to missing value imputation

Experimental data sets

Analysis of replicates

Consider a technical replication of the analysis of a biological sample i. The probability that a feature j is missing in the technical replicate, given that the observed intensity Y_i,j is equal to y in the original analysis, is defined as follows:

$$p_{i,j}^0(x)=\mathbb{P}({Y_{i,j}}^{\prime }=NA|Y_{i,j}=y)$$

with Y′_i,j the observed intensity in the replicate. The replicate data sets were used to infer the missingness probability function under the assumption that the probability p⁰_i,j(x) is independent of the biological sample and of the feature: p⁰_i,j(x) = p⁰(x). Then, observed intensities were stratified in 5% quantiles: $(y_0,\dots ,y_{20})$ and p⁰ was approximated by:

$$p\left({\displaystyle \frac{x_{\ell }+x_{\ell +1}}{2}}\right)=\mathbb{P}({Y_{i,j}}^{\prime }=NA|Y_{i,j}=[x_{\ell },y_{\ell +1}))={\displaystyle \frac{\mathbb{P}({Y_{i,j}}^{\prime }=NA,Y_{i,j}\in [y_{\ell },y_{\ell +1}))}{\mathbb{P}(Y_{i,j}\in [y_{\ell },y_{\ell +1}))}}$$

which was estimated by its empirical counterpart:

$${\displaystyle \frac{\sum _{j=1}^J\sum _{i=1}^8\sum _{r,r^{\mathrm{\prime }}=1,\dots ,7,r\ne r^{\mathrm{\prime }}}1{\mathrm{I}}_{Y_{i,j}^r=x,Y_{i,j}^{r^{\mathrm{\prime }}}=NA}}{J\times 8\times 7\times 6}\times {\displaystyle \frac{J\times 8\times 7}{\sum _{j=1}^J\sum _{i=1}^8\sum _{r=1,\dots ,7}1{\mathrm{I}}_{Y_{i,j}^r=x}}}}$$

with $Y_{i,j}^r$ the intensity of the protein j in the replicate r of the biological sample i.

Simulation study

A simulation study was conducted to illustrate the impact of both the nature of the biological difference between classes and the missingness mechanism. A realistic full data set was generated by kNN imputation of missing values on the log-transformed intensities from the ProteoCardis-cyto data set, after filtering of proteins with less than 10 non-missing values. Two classes of size 49 and 50 were randomly sampled among the 99 samples, so that no proteins are differentially expressed between the two classes except by chance. Then, difference of expression between classes was generated. First of all, 2,000 proteins were randomly selected to be truly differentially expressed, assuming either a difference in intensity (fold change), or a difference in probability of presence. Then, missing values were picked up assuming either MAR or MNAR mechanism, followed by filtering of proteins with less than 20 non-missing values. Details can be found in the Supplemental Material.

Combined test

We propose a protein level univariate combined test that accounts for both missing and non-missing data. Consider m features (here, proteins) observed in n samples belonging to two or more classes. For each feature j, the difference of intensity between the two classes on non-missing observations, and the association between class and missingness are tested via the following linear mixed model (lmm) and generalised linear mixed model (glmm):

$$\{\begin{array}{cc}{Y}_{i,j}=X_i{\beta }_j+Z_i{u}_j+{\epsilon }_{i,j},i=1,\dots ,n\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}{Y}_{i,j}\ne NA& (\mathrm{m}\mathrm{o}\mathrm{d}-\mathrm{l}\mathrm{m}\mathrm{m})\\ \log {\displaystyle \frac{\mathbb{P}[Y_{i,j}\ne NA]}{\mathbb{P}[Y_{i,j}=NA]}={X^{\mathrm{\prime }}}_i{{\beta }^{\mathrm{\prime }}}_j+{Z^{\mathrm{\prime }}}_i{u^{\mathrm{\prime }}}_j+{{\epsilon }^{\mathrm{\prime }}}_{i,j}i=1,\dots ,n}& (\mathrm{m}\mathrm{o}\mathrm{d}-\mathrm{g}\mathrm{l}\mathrm{m}\mathrm{m})\end{array}$$

with Y_i,j the log-transformed observed intensity of feature j in sample i, (X, X′) design matrices of fixed effect including the class effect, (Z, Z′) design matrices of random effects, (ɛ_i,j)_i and (ɛ′_i,j)_i i.i.d. centered gaussian. Let p₁ and p₂ be the p-values of the F-test for class effect in models (mod-lmm) and (mod-glmm) respectively.

Let S be the Fisher combined statistic defined as:

(1) $$S=-{\displaystyle \frac{1}{2}(\log p_1+\log p_2).}$$

S is large if at least one of the two p-values p₁ and p₂ is small; moreover, if only one of the two p-values is small, S is weakly affected by the value of the largest one. If the statistics of the two tests were independent, S would be χ²-distributed under the global null hypothesis, but this assumption may be violated, especially under MNAR assumption, since low protein abundances could simultaneously lead to low observed intensities and high probability to be missing. Therefore, the distribution under the global null hypothesis is obtained by repeated permutations (N^perm) of the classes. In order to reduce the computing time and to increase the number of distinct values that can be taken by S, the distribution under the null hypothesis is assumed to be identical for all variables with the same proportion of missing values. Mathematical details are provided in the Supplemental Material.

FSMs based on NA imputation methods

The combined test was compared with feature selection procedures based on missing value imputation. We considered the three imputation methods proposed by Tang et al. (2020a) (package metaFS), namely: (i) single value imputation, where all missing value are replaced by the smallest intensity observed in the data set; (ii) k-nearest neighbours (kNN); (iii) singular value decomposition (SVD). Following Wang et al. (2020), kNN was implemented using the R function SeqKNN with a number of neighbours k = 10 and SVD was implemented using the function pca (package pcaMethods) with two components. Then, the linear mixed model (mod-lmm) was applied on the vectors of observed and imputed intensities. The choice of a lmm testing procedure after imputation guarantees that the differences observed between the imputation-based FSMs and the combined test are exclusively due to the treatment of missing values. The imputation-based FSMs are denoted as follows.

• Single-lmm: log-transformation + single value imputation + lmm

• KNN-lmm: log-transformation + kNN imputation + lmm

• SVD-lmm: log-transformation + SVD imputation + lmm

Hurdle model

Goeminne et al. (2020) proposed a peptide level model, the Hurdle model, that presents similarities with our approach, by combining MSqRob, a mixed model applied on all peptides observed intensities of each protein including a random sample effect, and a quasi-binomial model on the number of observed unique peptides. The p-values are combined as in (1), but independence is assumed between both statistics. Functions to implement the tests are available at https://github.com/statOmics/MSqRobHurdlePaper.

Following Goeminne et al. (2020), additional filterings at the peptide level were applied to ProteoCardis and Pigs datasets prior to the Hurdle test implementation. Peptides with only one identification were deleted, as the model would be perfectly confounded. Then, proteins identified by only one peptide were removed.

Resampling-based procedure for control of false discovery rate (FDR)

To account for correlation between variables, the false discovery rate (FDR) was controled using the resampling-based procedure proposed by Reiner, Yekutieli & Benjamini (2003): p-values were reestimated by resampling (100 times) from the marginal distribution prior to p-value adjustment (Benjamini & Hochberg, 1995). Even if most FDR procedures only guarantee an upper-bound control and are subject to assumptions on dependence between variables, the number of selected variables for a given FDR threshold is an indicator of the FSM’s power. Moreover, the resampling-based FDR procedure considered here enables the bias due to complex dependence structures to be circumvented. Mathematical details are provided in the Supplemental Material.

Criteria for method comparison

Missingness mechanisms: both MAR and MNAR

The replicate data sets, including seven technical replicates for eight biological samples, allow for the assessment of the technical variability and the analysis of the MAR/MNAR hypotheses. Figure 1 and Figure S2 (left) display the average observed intensity of a protein as a function of the number of times it is missing in the replicate samples. The observed intensity decreased as the number of missing values increased, which suggests that the probability to be missing is higher when the protein abundance is lower, so missingness mechanisms is at least partially MNAR. In particular, we observed a pronounced difference of intensity between proteins with no missing value and protein with one, or a fortiori more than one missing values. Nevertheless, even when the protein was missing in a large proportion of replicates, the average observed intensity could still be high, indicating that missingness mechanisms are not exclusively MNAR. These observations were confirmed by the probability of being missing, that decreased when the observed intensity increased, but remained non-negligible even for consistent observed intensities (Fig. 1 and Fig. S2). For example, for an intensity equal to the median of the observed values on the data set, the probability of being missing was 0.23–0.25, and even for an intensity equal to the 0.9 quantile, the probability was still 0.05.

Comparison of FSM’s performances on simulated data

Figure 2 illustrates the ability to discriminate differentially and non-differentially expressed proteins of three FSMs, assuming two missingness mechanisms (MAR and MNAR) and under two types of difference of expression between proteins (differential presence and differential abundance). We considered two imputation-based FSMs: SVD based on global structure similarity, and single value imputation by the smallest observed intensity, followed by a linear (mixed) model, as well as the combined test. kNN-lmm was not considered since kNN imputation was used to generate the simulated data set. FSM’s performances were highly impacted by both the nature of the difference of expression and the missingness mechanism. Each imputation based FSM failed under one scenario: single-lmm did not achieve to detect differentially abundant proteins under MAR, while the AUC with SVD-lmm was close to 0.5 when proteins were differentially present under MNAR assumption. These observations are coherent with the underlying assumptions under each imputation methods: MAR for SVD and MNAR for single value imputation. In all scenarios, the combined test remains competitive, with AUCs close to the ones of the best of the two imputation-based FSMs.

Poor concordance between imputation-based FSMs based on MAR and MNAR

Comparison of the three imputation-based FSMs in terms of correlation between log-transformed p-values and proportion of common selected features (Fig. 3 and Fig. S3) indicated a strong agreement between SVD-lmm and KNN-lmm, but both methods showed a poor concordance with single-lmm. These observations were common to the ProteoCardis and Pigs data sets, but the concordance between methods was globally higher for Pigs due to a weaker proportion of missing values, which reduced the impact of the imputation method.

Imputation-based FSMs are concordant with either glmm on probability of missingness or lmm on observed intensities

In addition, the imputation-based FSM p-values were compared with the two tests involved in the combined test (Rows 2 and 3 of Figs. 4, 5 and Fig. S4). We observed a very strong correlation between the p-values of mod-glmm and single-lmm. Indeed, as the smallest intensity used for imputation was far from most of the observed intensities (Fig. S9), the proportion of imputed values among a class strongly impacted the average intensity after single value imputation (a large proportion of missing values automatically leads to a small average intensity after single value imputation). Therefore, testing the difference in fold-change and in missingness leads to consistent p-values. Single-lmm p-values were weakly consistent with mod-lmm on observed intensities. Conversely, the p-values from KNN-lmm and SVD-lmm correlated well with mod-lmm, but weakly with mod-glmm.

The combined test reaches a compromise between imputation-based FSMs

The combined test displays a strong agreement with single-lmm and a moderate agreement with KNN-lmm and SVD-lmm in terms of correlation between log-transformed p-values and proportion of common selected features (Fig. 3 and Fig. S3). This observation was confirmed by the scatterplots between log-transformed p-values of the combined test and each imputation-based FSM (first row of Figs. 4, 5 and Fig. S4). Indeed, features found highly significant by any of the imputation based FSMs were at least moderately significant with the combined test, while the opposite was not true. More precisely, for ProteoCardis data sets, the proteins with very low p-values (p < 5.10⁻⁴) with the imputation-based FSMs also had low p-values with the combined test (p < 5.10⁻²); Conversely for each imputation-based FSM, some non-significant proteins displayed very low p-values with the combined test. A more nuanced but similar assessment holds on Pigs, since most of the proteins that were found significant with SVD-lmm or KNN-lmm and non-significant with the combined test were very sparse (Fig. S5) and thus included a large proportion of imputed values, which indicates that imputation-based analysis on these variables is weakly reliable. On the contrary, the variables significant with the combined test but non-significant with SVD-lmm and KNN-lmm have a sparsity level that is either low or high.

Impact of the filtering threshold

The impact of the imputation method is expected to decrease with the proportion of missing values, which itself depends on the filtering threshold. Therefore, we repeated our analyses with higher filtering thresholds: 30, 40 and 50 for ProteoCardis, and 20 and 30 for Pigs, to examine to what extent the comparisons between FSMs were impacted. Figures S6–S8 display the concordance between FSMs for several filtering thresholds. The comparisons between methods still holds when threshold varies: SVD-lmm/KNN-lmm on the one hand, and single-lmm/combined test on the other hand were strongly concordant, while KNN-lmm/SVD-lmm were moderately concordant with the combined test, and poorly concordant with single-lmm. As expected, the agreement between all pairs of FSMs globally increased with the threshold.

A harder thresholding also resulted in a higher number of discarded variables, potentially leading to a loss of biological information. Indeed, Fig. S10 displays the distribution of the p-values as a function of the protein sparsity; we observed that very sparse proteins may exhibit very small p-values, in particular for ProteoCardis data sets. Moreoever, Table S1 indicates that the sets of most significant proteins for each FSM include 60–82% (average 70%) of features with more than half of missing values for ProteoCardis, and 17–57% (average 41%) for Pigs. Note that this proportion of sparse proteins among the selected ones is only slightly lower than in the entire data sets (79% for ProteoCardis and 43% for Pigs). This indicates that sparse proteins are selected almost as frequently as less sparse ones.

Similar FSM’s quantitative performances

Table 1 and Table S2 display the prediction accuracy with SVM and RF classifiers applied on selected proteins. Even though the four FSMs selected very different sets of features, their performances in terms of prediction were not significantly different.

Figure 6 and Figure S11 display the concordance between variable selection performed on independent data sets. On ProteoCardis-cyto, reproducibility of feature selection was similar with the four FSMs, while on ProteoCardis-env and Pigs the combined test and single-lmm outperformed imputation-based FSMs. Therefore variable selection based on the combined test and the single imputation FSM were equally reproducible and tended to be more reproducible than FSMs using structure-based imputation.

Figure S1 displays the number of selected variables for various values of FDR. The methods ranking varied with the data set and the FDR threshold, but the combined test remained competitive in terms of number of selected variables.

Comparison of the combined test with the peptide-level hurdle model

Similarly to the combined test, the Hurdle model proposed by Goeminne et al. (2020) targets simultaneously difference of intensity on observed values and difference in probability of detection, but with a peptide level model. As expected, both procedures lead to consistent p-values (Fig. 7), with Pearson correlation of log-transformed p-values of 0.637–0.786, but far from identical. In terms of feature selection, some proteins highly significant with one method exhibit non-significant p-values with the other, in particular for Pigs.

Regarding prediction, both procedures led to similar accuracy on ProteoCardis data sets, and on Pigs with RF, but the combined test had lower performances with SVM than RF, while the hurdle test maintained a similar accuracy (Table 1 and Table S2). No significant differences in terms of replicability were observed between the protein-level combined test and the peptide-level hurdle test (Fig. 8 and Fig. S12). More precisely, feature selection with the combined test was slightly more replicable than the hurdle test on Pigs, notably for the largest number of selected features, while the opposite was observed on the ProteoCardis data sets.

Missingness blends MAR and MNAR mechanisms and our method addresses both assumptions

Metaproteomics by LC-MS/MS generates a large proportion of missing values, usually imputed prior to statistical analysis. Several categories of imputation methods are routinely considered. Methods based on local similarity (e.g., kNN) or global structure (e.g., SVD) implicitely assume that missingness occurs independently of the true feature concentration (MAR). But the analysis of the replicate data sets clearly indicates that missingness is more likely to occur when the feature has a low abundance. On the other hand, the single imputation method relies on the assumption of a left censoring mechanism, but the distribution of observed intensities as well as the analysis of the replicate data set are not consistent with this assumption. Even if the co-existence of MAR and MNAR mechanisms is admitted in LC-MS/MS proteomics, exploration of their prevalence is often based on biological replicates (e.g., Karpievitch et al. (2009)), assuming similar protein abundances in all samples from a given class. This assumption is questionable, in particular in human gut metaproteomics characterised by a strong individual specificity. On the contrary, our analysis based on technical replicates guarantees that the true protein abundances are identical.

Limits of missing value imputation in metaproteomics

Missing value imputation is the standard method to account for missing values in metaproteomics. This flexible approach enables to address any type of statistical questions (e.g., prediction, network inference,…) using methods developed for data with no missing values. But the downside is the risk to “forget” which values were imputed and to treat them equally to observed values, regardless of implicit assumptions underlying imputation that can strongly impact biological findings when a large proportion of values are missing (O’Brien et al., 2018; Karpievitch, Dabney & Smith, 2012; Lazar et al., 2016). In particular, we observed that global and local structure imputations specifically led to selection of features with a large proportion of imputed and thus weakly reliable values. Therefore, despite its easiness of use, imputation has a cost in terms of reliability and should be limited to moderately sparse data sets. On sparse metaproteomics data, this condition would require to filter out a large part of the features, which may be harmful since we demonstrated that a large part of the potentially interesting proteins have more than half of missing values.

As an alternative to missing value imputation, censored statistical models developed for proteomics data can account simultaneously for MAR and MNAR mechanisms (Karpievitch et al., 2009; Luo et al., 2009; O’Brien et al., 2018). Moreover, Berg et al. (2019) proposed a multiple imputation model which handles MAR and MNAR assumptions, but their method suffers from a methodological bias since imputation is performed within each class, which may artificially increase significance of inter-classes differences by enhancing intra-class similarities. Similarly, Gianetto et al. (2020) proposed a combination of MAR and MNAR imputation for proteomics data, but their method relies on a log-normal distribution of the true abundances, unrealistic in the metaproteomics framework where a large part of proteins are truly missing. Even though these models may address the complexity of missingness mechanisms more acutely than simple imputation methods, they also heavily rely on assumptions regarding missingness mechanisms (e.g., hard thresholding) as well as signal distribution (e.g., additive effect, Gaussian distribution). Therefore, they can not be directly applied to metaproteomics data whose structure and characteristics strongly differ from proteomics.

The combined test reaches a compromise between imputation-based methods

Beyond differences in underlying assumptions, distinct imputation-based FSMs lead to very different sets of selected variables. While selection with the two MAR imputation methods (SVD and kNN) were consistent, these two methods displayed almost no agreement with the FSM based on MNAR assumption (single value imputation) for the highly sparse data sets ProteoCardis, and a low agreement for the moderately sparse data set Pigs. Therefore, the choice of an imputation method can strongly impact the biological conclusions. The combined test, which addresses the two types of missingness by combining a glmm on probability of presence (relevant under MNAR assumption) and a lmm on observed intensities after removal of missing values (amounting to assume MAR mechanisms) displays a correct agreement with each imputation-based FSM. In greater details, the features detected using single value imputation were recovered by the glmm, while the features selected using kNN or SVD were recovered by the lmm on observed values, which is consistent with the assumptions on missingness mechanisms.

Prediction accuracy is a classic criterion to compare FSMs (Tang et al., 2020b), but its relevance is questionable when the main interest is the biological interpretation of selected features. Indeed, as enlightened in our analysis, methods with similar classification abilities can lead to totally different selected sets and feature ranking, so choosing a FSM based on a slightly higher prediction accuracy that - may vary with the chosen classifier - seems hazardous.

Furthermore, the simulation study conducted on archetypical scenarios illustrated that the comparative performances in terms of discrimination between differentially and non-differentially expressed proteins strongly depends on the nature of the protein difference (associated to the biological context) and the missingness mechanisms (mostly dependent on the technology). First of all, this simulation study indicates that quantitative comparison of FSMs’ performances in a specific context has to be taken cautiously in another biological context, in which the proportion of each scenario, notably regarding biological differences, could vary. In particular, transposition of methods comparisons performed on proteomics data to the metaproteomics context can be hazardous, as both the biological context and the technical artefacts can differ. Secondly, while each imputation-based FSM fails in at least one scenario, the combined test remains efficient in all scenarios.

Thus, the combined test realises a compromise in terms of feature ranking and selection between inconsistent methods whose performances are based on unverifiable assumptions, and remains efficient in diverse scenarios. Therefore, it represents a robust solution, notably in the context of few prior knowledge regarding missingness mechanisms and the nature of biological differences.

Protein vs peptide level analysis

In proteomics, analyses can either be realised at the protein or peptide level. In univariate testing procedures on proteomics data, peptide level analyses including a run effect are usually considered as more robust and with higher power (Clough et al., 2012). Nevertheless, in metaproteomics, analysis are routinely performed at the protein level (see review by Tang et al. (2020b)). This choice enables simplifications in the analysis, but biological arguments could be considered: the complexity of metaproteome may lead to a higher proportion of peptide misidentifications; moreover, the individual specificity of the metaproteome generates a very high sparsity at the peptide level (e.g., 95–97% of missing values at the peptide level in the Proteocardis data sets) which is detrimental to the robustness of mixed model.

We compared the combined test with the hurdle model by Goeminne et al. (2020), which displays similarities but is implemented at the peptide level. Interestingly, contrary to what was demonstrated by the authors on proteomics data, this peptide level analysis did not demonstrate superior performances compared to the protein-level combined test on our metaproteomics data sets. Further analyses of the biological findings brought by each method would enrich the comparison.

Beyond metaproteomics, the generic aspect of the combined test enables to use it in other omics or non-omics contexts with data missing both at random and not at random.