A Meta-proteogenomic Approach to Peptide Identification Incorporating Assembly Uncertainty and Genomic Variation

Matching metagenomic and/or metatranscriptomic data, currently often under-used, can be useful reference for metaproteomic tandem mass spectra (MS/MS) data analysis. Here we developed a software pipeline for identification of peptides and proteins from metaproteomic MS/MS data using proteins derived from matching metagenomic (and metatranscriptomic) data as the search database, based on two novel approaches Graph2Pro (published) and Var2Pep (new). Graph2Pro retains and uses uncertainties of metagenome assembly for refer-ence-based MS/MS data analysis. Var2Pep considers the variations found in metagenomic/metatranscriptomic sequencing reads that are not retained in the assemblies (contigs). The new software pipeline provides one stop application of both tools, and it supports the use of met-agenome assembly from commonly used assemblers including MegaHit and metaSPAdes. When tested on two collections of multi-omic microbiome data sets, our pipe-line significantly improved the identification rate of the metaproteomic MS/MS spectra by about two folds, comparing to conventional contig-or read-based approaches (the Var2Pep alone identified 5.6% to 24.1% more unique peptides, depending on the data set). We also showed that identified variant peptides are important for functional profiling of microbiomes. All results suggested that it is important to take into consideration of the assembly uncertainties and genomic variants to facilitate metaproteomic MS/MS data interpretation.

A Meta-proteogenomic Approach to Peptide Identification Incorporating Assembly Uncertainty and Genomic Variation* □ S Sujun Li, Haixu Tang, and Yuzhen Ye ‡ Matching metagenomic and/or metatranscriptomic data, currently often under-used, can be useful reference for metaproteomic tandem mass spectra (MS/MS) data analysis. Here we developed a software pipeline for identification of peptides and proteins from metaproteomic MS/MS data using proteins derived from matching metagenomic (and metatranscriptomic) data as the search database, based on two novel approaches Graph2Pro (published) and Var2Pep (new). Graph2Pro retains and uses uncertainties of metagenome assembly for reference-based MS/MS data analysis. Var2Pep considers the variations found in metagenomic/metatranscriptomic sequencing reads that are not retained in the assemblies (contigs). The new software pipeline provides one stop application of both tools, and it supports the use of metagenome assembly from commonly used assemblers including MegaHit and metaSPAdes. When tested on two collections of multi-omic microbiome data sets, our pipeline significantly improved the identification rate of the metaproteomic MS/MS spectra by about two folds, comparing to conventional contig-or read-based approaches (the Var2Pep alone identified 5.6% to 24.1% more unique peptides, depending on the data set). We also showed that identified variant peptides are important for functional profiling of microbiomes. All results suggested that it is important to take into consideration of the assembly uncertainties and genomic variants to facilitate metaproteomic MS/MS data interpretation. Molecular & Cellular Proteomics 18: S183-S192, 2019. DOI: 10.1074/mcp. TIR118.001233.
Microbiome research has been applied to various studies of microbial organisms associated with different ecosystems, habitats and hosts (1)(2)(3)(4)(5)(6)(7)(8). Recent studies have further shown the broader impacts of microbiota on human health and diseases, including the influence of microbiota on the efficacy of cancer immunotherapy (9). Researchers started to investigate the therapeutic applications of microbiome, which are very promising as demonstrated in recent case studies: Zhu and colleagues edited gut microbiota on ameliorate colitis (10), and in another case, engineered commensal microbiomes were used for diet-mediated colorectal cancer chemoprevention (11). In microbiome research, metagenomic and metatranscriptomic (12,13) data often tell about the taxonomic distribution and potential functions, whereas metaproteomic data provides more direct information about the functionality of microbial communities (14,15,16,17). In a recent study of ocean microbiome (18), single-cell genomics and community metagenomics revealed that Nitrospinae are the most abundant and globally distributed nitrite-oxidizing bacteria in the ocean, whereas metaproteomic and metatranscriptomic analyses suggest that nitrite oxidation is the main pathway of energy production in Nitrospinae. And in another work, Kleiner and colleagues (19) showed that metaproteomics could be used to assess species biomass contributions in microbial communities, which is less prone to the biases in sequencing-based methods. Studies have shown the impact of microbiomes on human health and diseases. Metaproteomics provides a new opportunity of studying the functionality of microbial communities.
Although metaproteomics revealed complementary information to metagenomic and metatranscriptomic studies, it often applies conventional proteomics techniques using mass spectrometry (MS). Metaproteomic data analysis is challenging (20), which has motivated the development of recent algorithms and tools including MetaLab (21), cascaded search (22), Unipept (23), two-steps search (24), MetaPro-teomeAnalyzer (25) and ProteoStorm (26). Successful proteomic database search largely relies on the completeness and specificity of the target protein database. There are two broad categories of approaches for metaproteomic database search: approaches relying on a generic collection of reference proteins (for example, the human/mouse gut bacterial genes (27,21)), and approaches that use specific proteins/ peptides inferred from matching metagenomic and metatranscriptomic data sets of the same microbial community. Different from typical proteomics data analysis, peptide/protein identification from metaproteomics data generally lacks a good reference database of proteins for database searches, except for well-studied microbiomes. On the other hand, although matching metagenomic or metatranscriptomic data sets can provide more adequate reference database, they are often under-used for metaproteomic data analysis. Microbiomes are normally composed of many species, many of those species are closely related, as demonstrated in studies of microbial communities empowered by long sequencing reads (28,29). Despite of the recent advances in the development of assembly algorithms for metagenomes, it remains a challenge to assemble metagenomes from shotgun metagenomic sequences because of the nature of microbial communities: a large number of species co-exist in a sample; these species are of various abundances; and many species are closely related containing similar genome sequences with variations so that assemblers won't be able to untangle them. As a result, metagenome assemblies often provided limited references for peptide/protein identification from metaproteomic data sets.
Approaches have been proposed to address the abovementioned limitation of using metagenome assemblies as reference for peptide/protein identification from MS/MS data. One successful effort is from May and colleagues (30), who developed a method to derive metapeptides (short amino acid sequences that may be encoded by multiple organisms) from shotgun metagenomic sequencing reads, and then use the metapeptides as the reference for peptide identification from metaproteomic MS/MS data. They showed that by constructing site-specific metapeptide databases, they were able to detect more than one and a half times as many peptides as by searching against predicted genes from an assembled metagenome and roughly three times as many peptides as by searching against the NCBI environmental proteome database. Also some benchmarking experiments have been optimized to explore better experimental and computational strategies for metaproteomics (31,32) as well. We have also recognized the limitation of conventional approach of using assembled metagenomes (e.g. using only the contigs) for peptide/protein identification from MS/MS data and developed a graph-centric approach (Graph2Pro) (33) that uses the assembly graph to take into consideration of the assembly uncertainties. Our approach achieved significant improvement of peptide/protein identification from MS/MS data for microbiome research.
Comparing to (meta)genomic sequencing, metaproteomics (and proteomics in general) typically has relatively limited throughput. It is therefore important to develop methods that can extensively use the sequencing reads (in addition to assemblies) for peptide/protein identification from MS/MS data (sequencing reads from low abundant species are unlikely to be assembled into contigs). Further, microbial communities often contain closely related species and strains with genomic variations, many of which are of low abundances (28,29). Short sequencing reads containing variations are often "ignored" by the assemblers, so as a result, our graph-based approach that uses assembly graphs for MS/MS identification is unable to use these reads for MS/MS identification. The open database search approach (34) in principle can be applied to detect variations, which however will be impractical for metaproteomic identification because of the unknown and extremely large search space. The goal of our new variantaware approach is to incorporate the potential peptides encoded by the short sequencing reads that are not even assembled into the assembly graph, to improve the identification from MS/MS data. Instead of using variant callers (such as MIDAS (35) and metaSNV (36)) to detect potential variants from metagenomic and/or metatranscriptomic sequencing data (8), our variant-aware approach retains all potential variants containing peptides, and use them as reference for spectral database search. The rationale is that the variants are often of low abundance, so the typical variant detection tools that rely on sequencing coverage to distinguish variants from sequencing errors may not work. If a potential variant is supported by only few reads, but nevertheless is supported by the MS/MS data, it is a confident variant. We also note that, different from the previous method (30), which includes all potential short peptides from short sequencing reads, our variant-aware approach only includes the peptides that are similar to the proteins that have already been identified from MS/MS spectra. We showed that our approach effectively reduced the search space for MS/MS spectral database searching, and identified more peptides from metaproteomic MS/MS data.

EXPERIMENTAL PROCEDURES
Overview-We developed a pipeline for peptide/protein identification from metaproteomic mass spectrometry data, to optimize the use of matching genomic/transcriptomic data of the same microbial community. The pipeline is empowered by two novel approaches we developed, Graph2Pro (33) and Var2Pep proposed here. Graph2Pro uses the uncertainties of metagenome assembly captured in assembly graphs. Var2Pep uses genomic variations that are not retained in metagenome assembly. Specially, our pipeline uses two search databases for mass spectral identification: (1) the Graph2Pro database, which contains putative proteins and peptides inferred from the assembly graph (including the edges, i.e. the contigs, and the branching structures); and (2) the Var2Pep database, which contains putative peptides predicted from sequencing reads that potentially contain genomic variants as predicted by our new approach (see Fig. 1).
Rationale of the Variant-aware Approach (Var2Pep)-Our variantaware approach Var2Pep aims to retain peptides that can be predicted from unassembled short sequencing reads in the search database to improve MS/MS spectral identification. Unassembled short reads containing sequence variations can be mapped to assembled contigs as long as the contigs contain similar sequences (so the variant-bearing reads can be "used" for downstream analysis based on genomic sequences). However, for MS/MS analysis, a peptide with only a single amino acid difference (caused by a single nucleotide variation) may be overlooked if the peptide was not included in the target database, because two peptides with a single amino acid difference may result in very different mass spectra.
In principle, keeping all putative peptides in short sequencing reads in the search database for MS/MS analysis will address the abovementioned peptide variant problem. However, by doing so, we ex-pand the search database tremendously, resulting in not only the increase of searching time, but also the decrease of identification performance because of the extremely large database containing unbalanced target and decoy entries. By contrast, Var2Pep only considers reads-derived putative peptides that share significant sequence similarity with the proteins/peptides in the Graph2Pro database.
Var2Pep and Its Integration with Graph2Pro-As illustrated in Fig.  1, our pipeline for protein/peptide identification from metaproteomic data integrates the Graph2Pro approach (33) and the new Var2Pep approach. The Var2Pep approach is composed of several steps, as shown in Fig. 1. The first step is to extract short sequencing reads that cannot be mapped to contigs, and the reads that can be mapped to contigs but with mismatches, based on bowtie2 (37) mapping results. In the second step, putative peptides are predicted from the reads with mismatches and unmapped reads using FragGeneScan (38). The third step is to prepare a new search database (called Var2Pep) that contains the putative peptides sharing high similarity (e.g. Ն 70% identity by default, based on RAPSearch2 (39) search results) with proteins/peptides in the Graph2Pro database. The 70% sequence identity is chosen, considering that it allows up to 3 substitutions in identified variant peptides (the average length of identified peptide is 10). This identity cutoff works well in practice, but users can customize this identity threshold when they apply our pipeline.
The MS/MS database searching is applied to Var2Pep database, controlled by false discovery rate (FDR Յ 1%). When combining the Graph2Pro and Var2Pep search results to derive a non-redundant collection of identified peptides, we chose the peptide identification with the higher score for a spectrum when the two searches result in different peptides for the same spectrum.
Data Sets-We tested our approaches using two publicly available collections of meta-omics data sets. The first collection (called the ocean data set) contains meta-omics data sets derived from ocean water samples (30), which were collected from Bering Strait (BSt) chlorphyll maximum layer and from the more norther Chukchi Sea (CS) bottom water. This collection contains matching metagenomic and metaproteomic data. The LC-MS/MS spectra from triplicate acquisitions of peptides from the BSt and CS samples are designated as BSt 45-47, and CS 51-53. We downloaded the raw sequencing data and metaproteomic data from the following website: http://noble.gs. washington.edu/proj/metapeptide/. For the meta-proteogenomic analysis, the same metagenomic data set (BSt) was used for analyzing the three metaproteomics data sets (BSt45, BSt46 and BSt47), and similarly the same metagenomic data set (CS) was used for analyzing metaproteomic data sets CS51-53.
Metagenomic/Metatranscriptomic Assembly-Raw reads were preprocessed using Trimmomatic (version 0.32) (43) and only reads of at least 80 bps were used in downstream analyses. We used SOAP-denovo2 (44) and MegaHit (45) for the assembly of metagenomic and/or metatranscriptomic sequences. The default parameters of SOAPdenovo2, and the k-mer size of 31 were used. For MegaHit, we used -k-list 21,29,39,59,79,99 and default setting for other options; assembly graphs from the last k-mer size (i.e. 99) were used as inputs to our pipeline. Considering metaSPAdes (46) is becoming a popular assembler for metagenomes (a recent study showed that meta-SPAdes produced longer scaffold, but had the longest run times (47)), we extended our Graph2Pro such that it can exploit assembly graphs produced by MetaSPAdes.
Peptide/Protein Identification from MS/MS Data Using Database Search-We used the MS-GFϩ search engine (Version v10089) (48) in our pipeline. We used the following parameters for the MS-GFϩ database searching: (1) instrument type was set as high-resolution LTQ; (2) precursor mass tolerance was set as 15 ppm; (3) Isotope error range was selected as Ϫ1,2; (4) modification was set to include at most 3 modification on the sequence, including variable oxidation on Methionine and fixed carboamidomethy on Cysteine; and (5) maximum charge was set as 7 and minimum charge was set as 1, (6) Number of tolerable termini is semi-tryptic. The false discovery rate (FDR) was estimated by using a target-decoy search approach (49), in which the decoy proteins were generated by reversing the target sequences. Our pipeline outputs two sets of identification results after imposing FDR at the spectrum level and at the peptide level, respectively. RESULTS We tested our integrated pipeline for peptide/protein identification using the ocean water and wastewater multi-omics data sets. We compared our approaches with contig-based (in which proteins predicted from contigs are used as the reference), read-based approach (using metapeptides derived from reads as the reference), and a simple combination of read-and contig-based approaches (in which peptides were derived using the read-and contig-based approaches separately, and then merged as the final results). Our results showed that by considering variations and uncertainties of assembly, we can significantly improve the MS/MS spectral identification rate, and therefore identify more peptides (and proteins) for downstream analyses including functional profiling. The Var2Pep approach alone identified 5.6% to 24.1% more unique peptides, and our pipeline integrating the Graph2Pro approach and the Var2Pep approach improved the overall identification rate of the metaproteomic MS/MS spectra by about two folds, as compared with conventional contig-or read-based approaches.
Performance of the Meta-proteogenomic Approach on the Seawater Microbiome Data Set-We first tested the pipeline for peptide and protein identification from metaproteomic MS/MS data using the ocean data set. Table I lists the size of  Table II summarizes the peptide identification for selected data sets (BSt45 and CS51; see results for all data sets in the supplementary Data File S1) using different approaches. Fig.  2 shows the comparison results for all ocean water data sets. Briefly, our approach approximately doubled the number of identified peptides as compared with the contig-based and read-based approaches, when the same FDR (Յ 0.01, at spectrum level) estimated using the target-decoy search approach was applied. When FDR Յ 0.01 at peptide level was applied, all methods unsurprisingly identified fewer peptide, however, we observed similar trends of improvements by Graph2Pro and Var2Pep: supplemental Fig. S1 shows the comparison of the different approaches on all ocean water data sets, using FDR Յ 0.01 at peptide level, and supplemental Fig. S2 shows the side-by-side comparison of the results by using FDR at spectrum or peptide level for the CS51 data set, respectively. It was shown in reference (30) that readbased approach using short amino acid sequences predicted from shotgun sequencing reads as the search database identified significantly more peptides as compared with the contig-based approach using genes predicted from contigs assembled by SOAPdenovo version 1.06 (30)). Our results confirmed the observation that when SOAPdenovo2 was used as the assembler, contig-based approach resulted in worse MS/MS identification (see Table II) than the read-based approach. However, when MegaHit was used as the assembler, the trend reversed and the contig-based approach identified more peptides than the read-based approach (see Fig.  2), whereas in both cases adding Graph2Pro and Var2Pep steps further improved the peptide identification.
Our pipeline significantly outperformed the simple combination of contig-based and read-based approaches (red bars versus blue bars in Fig. 2). For instance, for the Bst45 data set, reads-based method identified 7795 (8.6% of total) spectra, and the contig-based (using MegaHit) approach identified 8631 (9.5%) of total spectra. Combining contig-based search and reads-based search resulted in the characterization of a total of 16,426 spectra/3813 peptides. By comparison, our pipeline reported the identification of comparable number of spectra (15,913), but many more unique peptides (6240).
Comparison of Our Approach with Grouped and Cascaded Searches-Our pipeline uses separate database searches on the Graph2Pro and Var2Pep databases, and then combines the search results (we called it separate search for comparison purpose). We compared this approach with the other two approaches: a single spectral search on a grouped database combining both databases (called grouped approach), and cascaded search using only rejected spectra from the spectral search against the Graph2Pro database for the subsequent search against Var2Pep database (i.e. the cascaded approach, similar to the approach developed by Kertesz-Farkas and colleagues (22)). Table III summarizes the comparison of separate and cascaded searches at the same FDR level. For example, for CS51, our pipeline based on separate search against the Var2Pep database identified additional 2073 spectra (917 unique peptides); by comparison, the cascaded search against the Var2Pep database identified slightly fewer 1598 spectra (629 unique peptides). Notably, this result appears different from the previous report that cascaded search identified more spectra and unique peptides (22) than the separate search. One possible reason is that removing some high qualify spectra in the cascaded searches may alter the overall score distribution, and thus the estimation of false discovery rate for different searching results. Our approach also outperformed (marginally) the search using combined database. For example, for the CS51 data set, combined database search only identified a total of 20,046 spectra and 7662 unique peptides; by contrast, our approach based on separate searches identified 31,145 spectra and 12,380 peptides, and the cascaded search identified 30,670 spectra and 12,092 peptides from the same spectra data set. Based on our tests, we recommend using either separate or cascaded search, but not using the combined database for spectrum search. Our pipeline supports both separate search (default) and cascaded-search, and the results reported below were based on the separate searches.
Performance of the Meta-proteogenomic Approach on the Wastewater Microbiome Data Set-Next we tested our approach using the wastewater data sets. Fig. 3 illustrates the In all cases, the contig-based approach achieved significantly better performance than the read-based approach, and our approach integrating Graph2Pro and Var2pep significantly outperformed other approaches.
As shown in Fig. 3, significant improvement of peptide identification was achieved by using both metatranscriptomic and metagenomic data sets as the reference as compared with using only metagenomic sequencing data (see supplementary Data File S2 for details). This result is expected, because metatranscriptomic data may provide sequence information for otherwise low-abundant species that may have been missed by metagenomic sequencing. We note that among the wastewater data sets, the Var2Pep approach resulted in the most significant improvements of the peptide Note: S* represents SOAPdenovo2, M* represents MegaHit, PSM stands for peptide spectrum match. Bst stands for Bering strait, and CS stands for Chukchi sea. Graph2Pro (S*) and Graph2Pro (M*) represent using assembly graph from SOAPdenovo2 (S*) and MegaHit (M*) as the reference in Graph2Pro, respectively. In all cases, FDR (false discovery rate) was estimated using a target-decoy search approach, and a cutoff of 1% at spectrum level was applied. This table only shows the results for two datasets. See Fig. 2 and supplementary Data File S1 for results of all data sets.

Meta-proteogenomic Approach to Peptide Identification
Molecular & Cellular Proteomics 18.14 S187 identification for the SD3 data set, with an additional 24.1% and 19.5% more identified unique peptides when using matched metagenomic data set alone (SD3-MG), and both the metagenomic and metatranscriptomic data sets (SD3-MGMT).
Considering the increasing popularity of metagenome assembler MetaSPAdes, we also used the wastewater data sets to test if more contigs assembled from metagenome resulted in better identification of metaproteomics data. Table IV shows the comparison of the results for SD3-MG data set using assemblies from MegaHit and MetaSPAdes as the reference in our pipeline. Strikingly, MetaSPAdes produced about three times more contigs (more than doubled the total bases in the contigs) as compared with MegaHit. However, the increasing of contigs did not lead to more peptide identifications.
Abundant Known (and Unknown) Proteins Revealed by Metaproteomics-We examined proteins that were identified using our pipeline. For example, for CS51 data set, its first most abundant protein was supported by 270 spectra (0.8% of the total 31,145 identified spectra for this data set). A total of 663 proteins or protein fragments (of minimal length of 60 aa) were each supported by at least 10 spectra (see supplementary Data File S3 for details). Functional annotation of these proteins by myRAST (http://blog.theseed.org/down loads/myRAST-Intel.dmg) (50) revealed 90 ribosomal proteins, 61 heat shock proteins (HSP60), 59 elongation factors (GTP_EFTU), 38 chaperone proteins (DnaK), 30 glutamate aspartate periplasmic binding protein precursors, 26 branched-chain amino acid ABC transporters, 10 TonB-dependent receptors, among other functions (see supplementary Data File S3 for details of the annotations). Our result is consistent with the fact that many of these functions (e.g. ribosomal proteins) are typically the most abundant and highly conserved in prokaryotic genomes.
Among the CS51's 663 proteins each with support of at least 10 spectra, a significant fraction (93, 14%) have no functional annotation based on myRAST annotation, including a protein (ID: NODE_17410_length_722_cov_188.0000_ ID_34819_1_456_ϩ) that has the third most spectra support with 257 spectra and many others with strong spectral supports. See supplementary Data File S3 for the complete list of these highly expressed proteins that lack any functional an-notations (and their sequences can be found at our supplementary website). We argue that these proteins are most wanted proteins that wait to be further studied experimentally for their functions. Fig. 4 illustrates two examples of proteins supported by a significant number of variant peptides (details of the variant peptides are available at our supplementary website). Protein366 was predicted from the assembly graph by Graph2Pro and was shown to be supported by a total of 27 peptides, including six variant peptides. The myRAST annotation did not reveal any significant hits for this protein. We also did sequence similarity search of this protein using ncbi blast (http://www.ncbi.nlm.nih.gov/blast) against the nr database (as of Nov 18th, 2018), which returned no hits either. Although no significant sequential similarity was detected between Protein366 and proteins in the public databases, a structural comparison (by FATCAT (51)) of this protein's predicted structure (using QUARK (52)) revealed that it contains a Cystatin-like fold composed of helical segments packed against coiled antiparallel beta-sheet. Protein2811 is another example, which was annotated as a putative peptide/opine/ nickel uptake transporter. This protein has a total of 17 supporting peptides, among which nine are variant peptides.
For the wastewater data sets, an interesting example is a protein (Protein264) identified from SD6 data set. This protein has the most spectra support (2266 spectra), and the domain prediction shows that it contains three SLH (S-layer homology) domains and another larger domain that has unknown function (DUF2815). The presence of SLH domains in this protein suggest that it is likely a cell-surface protein, as studies have shown the utilization of SLH domains as vehicles for surface location of function proteins (53,54); however, the precise function of this protein remains to be determined. Fig.  4 (bottom) shows the locations of the domains and the distribution of identified peptides and variant peptides supporting this protein.
Low Abundant Proteins with Supports Augmented by Variant Peptides-As shown above, metaproteomics is likely limited to revealing highly abundant proteins expressed from highly abundant bacterial species, so metaproteomics-based functional profiling of microbiomes will only provide a shallow glimpse into the functionality of the underlining microbial community. A hope to alleviate this limitation is to augment the identification of proteins using variant peptides. Here we use the wastewater data set (SD3-MG) as an example to illustrate this approach. We used BlastKOALA (https:// www.kegg.jp/blastkoala/) (55) to annotate the functions of identified proteins from this data set. A total of 1401 KEGG families (i.e. the KO families) were identified using proteins with at least one peptide support; however, this number increased to 1938 if variant peptides were also considered (a 38% improvement). We acknowledge that the number of families is still small, because of the shallow functional profiling of microbiomes based on metaproteomics. We note that the variant peptides identified by our pipeline have both nucleotide level (from metagenomic and/or metatranscriptomic sequences) and mass spectral level supports. Here, we examined the types of substitutions found in the variant peptides. Not surprisingly, we found that most variations are the substitutions likely to be found in homologous proteins, with positive BLOSUM scores (BLOSUM matrices contain log-odd scores of all possible substitution pairs of amino acids, where positive scores represent favorable substitutions observed in homologous proteins such as Asp to Asn substitution and negative scores represent unfavorable substitutions such as Arg to Asp) (56). As shown in Fig. 5, the average BLOSUM (BLOSUM62) score between two amino acids is Ϫ2 (see "all" box), whereas the average BLOSUM score of substitutions in variant peptides is positive for Pro-tein366. This result indicates that Var2Pep likely detected biologically meaningful variations.

DISCUSSION AND CONCLUSIONS
We developed the Var2Pep algorithm, which combined with Graph2pro achieved drastic improvement of peptide identification from metaproteomics data. Var2Pep makes use of shotgun sequencing reads for MS/MS identification, and it reduces the search space imposed by the large number of sequencing reads by only keeping peptides that share similarity with proteins already identified by Graph2Pro. As shown in Table I, there were 15.9 million putative peptides in BSt data set from the mismatched/unmatched reads, and Var2Pep only contained 2.7 million peptides for MS/MS database search. This process significantly reduced the number of possible decoys in the database to permit the database searches, which significantly contribute to increasing the identification rate of metaproteomic MS/MS data. Similar ideas have been applied to improve metaproteomics data identification: for example MetaLab (21) uses a reduced gut reference gene catalogue to improve MS/MS identification). We note Var2Pep may still be used when assembly graph is not available (i.e. only the contigs are available). In such a case, peptides encoded in short sequencing reads that share similarity with MS/MS supported peptides identified from the contigs can be included in the database searching for MS/MS analysis.
The main goal of our approaches is to optimize the use of matching metagenomic and/or metatranscriptomic data for metaproteomics data analysis. Our integrated pipeline results in two databases for MS/MS identification, one from Graph2Pro (Graph2Pro DB), and the other from Var2Pep (Var2Pep DB). Our pipeline provides one way of using these databases as shown in Fig. 1. Users can make use of these databases in different ways. In order to get high confident peptide identification, as well as variant peptides, it is highly recommended that users apply strict false discovery rate control, either in spectral level or peptide level. A reasonable concern on our approach of combining the results from two separated searches each at 1% FDR is that the actually FDR would be higher. We re-calculated the FDRs for combined results, based on the unique sets of identified spectra for all the data sets we tested. The results showed that the FDRs estimated using this approach were slightly higher, but still reasonably low (at about 1.5%).
It is possible to combine our approaches with those that are recently developed, including the utilization of spectral clustering to speed up the search as shown in (57). We used the MS-GFϩ search engine (48), which is one of the fastest MS/MS search engines. More recently developed tools for the fast database searching of metaproteomic MS/MS data such as ProteoStorm (26) can also be incorporated into the pipeline. Our significantly improved results indicate there is still a huge space for algorithmic improvement in the field of metaproteomics and proteomics in general.
The goal of metaproteomics is not only to identify proteins expressed in the microbial community, but also to estimate their abundances (i.e. their expression levels) under different conditions. Nevertheless, a protein can be quantified only if it can be identified by using the metaproteomic data. Therefore, the methods presented here that increase the coverage of protein identification will also help the subsequent steps for protein quantification. Furthermore, identification of peptide variants may contribute to accurate quantification of the proteins from which these peptides are derived. Tools such as MetaGOmics (58) can be applied to infer functional and taxonomic contents of the microbial communities based on identified proteins/peptides. Better identification of peptides and/or proteins from metaproteomic MS/MS data can improve the functional and pathway annotation of the corresponding microbial communities, as we have shown in (33).
We showed that metaproteomics (and matching metagenomics/metatranscriptomics) data can be used to identify hypothetical proteins that sometimes are very abundant in the microbial communities. These MS/MS supported hypothetical proteins provide interesting candidates for further studies of their functions using computational and experimental approaches. However, our analyses of the identified proteins also suggested that unless we significantly increase the depth of the metaproteomic experiments, what we identify are limited to those highly abundant (often well studied) proteins.
In conclusion, we developed a pipeline for metaproteomic MS/MS data analysis using matching metagenomic and/or metatranscriptomic sequencing data as the reference. Tests of our pipeline using publicly available meta-omics data sets showed that it is important to consider the assembly uncertainties (captured in assembly graph) and genomic variants to maximize the utilization of metagenomes and/or metatranscriptomes as the reference for metaproteomic data interpretation.

DATA AVAILABILITY
Our pipeline for peptide/protein identification from metaproteomic data using matching metagenomic and/or metatranscriptomic data as the reference can be downloaded from https://github.com/COL-IU/graph2pro-var. We also made available the results (including the reference databases and search results) from analyzing the two collections of microbiome data sets in a supplementary website at http://omics. soic.indiana.edu/MP/and Zenodo at http://doi.org/10.5281/ zenodo.2691363.