Ayu: a machine intelligence tool for identification of extracellular proteins in the marine secretome

Domain: earth_science
Published: 2025-03-21
Editorial outcome: minor_revision

Abstract
Microbes are the engines driving the elemental cycles. In order to interact with their environment and the community, microbes secrete proteins into the environment (known collectively as the secretome), where they remain active for prolonged periods of time. Despite the environmental relevance of microbes, our knowledge of the marine secretome remains limited due to a lack of effective in silico methods for the study of secreted proteins. An alternative approach to characterise the secretome is to combine modern machine learning tools with the evolutionary adaptation changes of the proteome to the marine environment. In this study, we identify and describe adaptations of marine extracellular proteins, which vary between phyla, resulting in differences in ATP costs, amino acid composition and nitrogen and sulphur content. We develop ‘Ayu’, a machine prediction tool that does not employ homology-based predictors and achieves better and quicker performance than current state-of-the-art software. When applied to oceanic samples (Tara Oceans dataset), our method was able to recover more than double the proteins compared to the most widely used method to identify secreted proteins. The application of this tool to open ocean samples allows better characterisation of the composition of the marine secretome.

Introduction
The marine environment is the stage for critical geochemical processes key for maintaining planetary habitability, such as the production of atmospheric oxygen, the remineralization of organic carbon and the cycling of nitrogen, phosphorus and sulphur. Microbes have a central role in this process, as their genomes code for the enzymes that catalyse these chemical reactions. It is not surprising, then, that a large research effort has been poured into understanding the capabilities of marine microbe communities and their metabolic and genomic capabilities, resulting in global sampling expeditions like GEOTRACES, Tara Oceans and Malaspina.

However, in our quest for understanding what microbes are capable of, the question of where seems to be key. Studies in the laboratory with pathogenic bacteria indicate that up to 30% of the bacterial genome encodes proteins that are released to the extracellular milieu. This subset of proteins (the secretome) is how bacteria interact with their environment, and as such are involved in a myriad of processes such as provision of nutrients through recognition, degradation and uptake of large extracellular molecules; communication and competition between bacterial cells; bacterial adhesion and biofilm production. Exoproteome composition is also a key mechanism for bacterial adaptation. Studies of the exoproteome of Ruegeria and Synechococcus strains revealed that the exoproteome composition of each strain reflected their adaptations to their ecological niche and their growth conditions. These discoveries have been corroborated in metaproteomics studies on both epipelagic and bathypelagic seawater samples. Measurements of the extracellular enzymatic activity in the ocean indicate these reactions are largely catalysed by dissolved (cell-free) enzymes, with this ratio of cell-free to cell-attached enzymatic activity increasing with depth. Furthermore, these experiments have also shown that the effects of the exoproteome can extend beyond the cell that secreted them, as exoproteins present a half-life of up to 20 days away from the cell, suggesting that in order to understand the nutrient utilisation capabilities of a marine community, the history of the water mass might be more important than the present community composition.

Despite the relevance of the secretome, its study is limited by the lack of appropriate methodology. Most marine prokaryotes cannot be cultured in lab conditions, and although shotgun metaproteomic approaches have been used to study the exoproteome, their throughput is low compared to proteomic assays based on bacterial cultures and most of the material recovered belongs to virions. Even proteomics assays in controlled environments present difficulties, as without careful methodology and quantification, it is difficult to determine which proteins found in the exoproteome are secreted or merely a product of cell lysis. A reasonable approach would be to exploit the vast amounts of metagenome and metatranscriptome datasets available, but we are faced then with the challenge of predicting subcellular localization from the amino acid sequence.

Although the popularity of artificial intelligence and machine learning has led to the development of many tools for this purpose, a recent review by Hui et al. compared several state-of-the-art available tools for Gram-negative bacteria and argued that “More enthusiasms have been put in new algorithms rather than the biological side” [sic]. Many of these tools are too narrow in their scope, focusing only on one secretion system or one bacterial strain, which limits their use in metagenomic data. Furthermore, many of these tools are only available to the scientific community as web servers, which are not suited to the high volumes of data required for omics datasets. Even pSORTb 3.0, the gold standard for subcellular localization predictions, relies on homology searches against its curated profile and protein databases to obtain its predictions, which severely limit the throughput of proteins it is able to process. It is no wonder, then, that most publications studying the secretome in omics datasets almost exclusively use SignalP (a mature and robust software to detect Sec/Tat signals) to identify secreted proteins, even though translocation through the general secretion pathway does not guarantee secretion to the extracellular milieu, and many of the proteins with signal peptides stay in the periplasm or attached to the outer membrane.

Yet, the peculiarities of the marine environment present an opportunity to improve protein localization predictions. It is known that the amino acid composition (AAC) of a protein is in part adapted to the physicochemical properties of its location. Seawater exhibits many distinguishing features, chief of which is its average salt concentration of 3.5%. Salt has a denaturing effect in non-adapted proteins, mainly attributed to the disturbance of the water layer surrounding the protein, lowering solubility and promoting protein aggregation. Cytoplasmic proteins are protected from the effects of salts, as marine bacteria are salt-out strategists which maintain a relatively salt-free cytoplasm. However, this is not the case for proteins which operate in the periplasm, which is not osmoregulated, or in the extracellular milieu. Consequently, both extracellular and periplasmic proteins must be adapted to this salt content. Previous studies support this point: an analysis of proteomes found a correlation between isoelectric point (pI) and salinity, a survey of the proteins of the halophile gammaproteobacteria Chromohalobacter salexigens revealed that only the periplasmic and secreted proteins presented adaptations to salt, and a study comparing close phylogenetic neighbours of freshwater-saltwater pairs found differences in the isoelectric points and AAC of their encoded proteins, with differences being more pronounced in secreted than in cytoplasmic proteins.

Hence, in this study we first characterise the specific adaptations of extracellular proteins to the marine environment. Then, with this information, we developed a machine learning tool (‘Ayu’), designed to exploit the signal left by these adaptations to predict secreted proteins in large marine metagenomic datasets, comparing its performance to state-of-the-art tools for subcellular location prediction. And finally, we applied this tool to environmental samples for the Tara Oceans expedition to uncover how much and which proteins composed the actual marine secretome.

Results
Differences in AAC between subcellular localizations and habitat

A biplot of the weighted log-ratio analysis of grouped AACs for marine proteins can be found in Supplementary Fig. S1. Although this collection of logratios does not explain a majority of the variance (42.5% within the first two components), the plot shows a gradient of changes from cytoplasmic to extracellular proteins, with periplasmic proteins situated between them. The loadings from the plot can be used to determine which amino acids are driving these differences: i.e., extracellular proteins are relatively enriched in polar (Ser, Thr, Asn, Glu), negatively charged (Asp, Glu) and aromatic (Phe, Tyr, Trp) amino acids. In contrast, cytoplasmic proteins contain more positively charged (Arg, Lys, His) and small hydrophobic (Val, Ile, Leu, Met, Cys) amino acids.

However, these changes could just reflect the different physicochemical properties of the cytoplasm compared to the other subcellular locations, and not to any specific effect of the extracellular environment. Therefore, we compared the AAC of non-marine proteomes to our marine representatives. As our control group, we chose the ESKAPEE pathogens (Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter baumanii, Pseudomonas aeruginosa, Enterobacter spp. and Escherichia coli), as they are an extensively studied group of bacteria, and none of their members are commonly found in marine environments. Figure 1A shows radar plots comparing the AAC of our marine dataset and the proteome of the seven ESKAPEE pathogens. While both groups follow the systemic differences previously described for the marine proteomes, further differences between both groups of proteomes can be observed depending on the subcellular location (Supplementary Fig. S2). While cytoplasmic proteins from both groups present a similar composition, periplasmic and extracellular marine proteins show an increase in negatively charged and aromatic amino acids, and a decrease in some positively charged (Arg, Lys) and hydrophobic (Ala, Leu, Val) amino acids, compared to ESKAPEE proteins found in the same subcellular location.

To test if these differences were statistically significant, a Dirichlet regression analysis was performed, adjusting for subcellular localization and cell wall structure (that is, if the proteins were coded by a gram positive or gram negative bacterium). The complete results of the analysis can be found in Supplementary Text S1. The effect of habitat was modelled as an interaction between habitat and subcellular location, with the objective of assessing how each cellular compartment was affected separately. This interaction was found to be statistically significant by comparing the full model to a reduced model with the interaction removed (difference in deviance = 1794, p-value < 2.2e −16). Overall, the model indicates that while subcellular location and cell wall structure influence AAC, these differences are intensified for marine proteins relative to ESKAPEE proteins, especially for periplasmic and extracellular proteins. Amino acids E, G, I, K, P, T, V are especially affected, with their effects increasing at least 20% (p-value < 1e−5 for all comparisons). Additionally, the marine environment affected periplasmic and extracellular proteins differently, as the AAC of extracellular proteins was influenced more strongly towards polar, aromatic and charged amino acids. When considering both the changes in AAC explained by subcellular location and the influence of the marine environment, we find that in 12 amino acids (L, V, F, Y, S, T, N, Q, D, E, K, R) the distribution varies following their order of exposition to the environment (Extracellular > Periplasmic > Cytoplasm), compared to 7 in the non-marine ESKAPEE proteins (L, Y, G, S, T, N, R).

Overall, these results prove that the marine environment has a specific effect on the proteins exposed to it. However, the question of which effectors might force these constraints still stands. Salinity is the obvious culprit, as it both has a well-known effect in protein function and it is known that bacterial cytoplasms are relatively salt-free compared to the periplasm and the extracellular milieu. In fact, the pattern of amino acid substitution observed in marine extracellular proteins (Decrease in positively charged and aliphatic amino acids, increase of negatively charged and small amino acids) is remarkably similar to that reported for salt adaptation in halophilic proteins. The increase in aromatic residues is the exception to the rule, but there are studies reporting a link between hydrophobic interactions between aromatic residues and enzyme halotolerance, possibly by forming weak polar interactions with other residues. Another possible explanation is that the increase in aromatic residues is a product of an adaptation to multiple stresses. In that respect, there are multiple reports where the addition of aromatic residues to a marine protein increased its tolerance to high and low temperatures.

Differences in AAC between taxonomic groups

As the Dirichlet regression model indicated a significant taxonomy effect in AAC, we inspected the proteomes of a few taxa with enough proteins (n > 30) in each location to provide a fair estimate. Supplementary Fig. S3 shows a Log Ratio Analysis biplot that exposes this interaction between taxonomy and AAC, with differences between phyla caused by the proportion of positively charged (His, Arg, Lys) and small apolar amino acids (Ala, Gly). Although the amount of proteins available only allow for comparisons at the level of order, we can dig deeper by performing analyses on groups of interest. Supplementary Fig. S4 shows the amino acid distribution in the orders Synechococcales, Flavobacteriales, Bacilliales and Micrococcales. Comparing the four sets of extracellular proteins, it can be observed that while all groups follow the trends described previously (i.e., increase in negatively charged and polar residues compared to positively charged and hydrophobic residues as we increase exposure to the extracellular milieu) each order follows its own amino acid distribution. For example, Synechococcales adapt to the extracellular medium by increasing the ratio of Glutamic acid at the expense of polar amino acids, while extracellular proteins in Flavobacteriales and Rhodospirillales (Supplementary Fig. S5) are enriched in polar amino acids.

To explain these differences in AAC across phyla, we might turn to the differences in lifestyle between the shown taxa. Previous studies have speculated that extracellular proteins produced by bacteria will be on average cheaper than their cytosolic counterparts, as these proteins cannot be recycled, and because cooperation is more likely to evolve when the costs of said cooperation are outweighed by the advantages provided by their kin in the community. Here, we have found that the cost of extracellular proteins varies significantly between phyla. We propose that these differences originate from the different trophic strategies of the bacteria that produce them. It is not unreasonable that Alteromonas, a copiotrophic bacteria with chemotaxis and motility to detect locations with high nutrient availability, will get a better return for their investment in the production of extracellular enzymes than an oligotrophic organism such as Prochlorococcus, which cannot guarantee that its extracellular proteins will provide a profit. Extracellular enzymes are an important part of bacterial adaptation to their ecological niche, as their genetic properties reveal: extracellular protein repertoire is highly divergent compared to intracellular proteins, and they tend to be lost and gained at a high rate in bacterial genomes. As different trophic strategies will be under different selection pressures, it would stand to reason that extracellular proteins from oligotrophic organisms are selected for cost.

Differences in protein properties derived from AAC

As AAC is the basis for all characteristics of a protein, it would be reasonable to expect that the differences between subcellular locations be reflected in other protein properties. The classical example of this would be the isoelectric point (pI), which is a combination of the dissociation constant (pKa) values of the constituent amino acids. In fact, it is already well known that bacterial proteomes present a bimodal distribution, with both peaks corresponding to cytosolic and integral membrane proteins. In the marine dataset, we found an evident shift towards the acidic end of the scale as the location gets closer to the extracellular milieu (Fig. 1B). These results are consistent with those found in proteomes of closely-related bacterial lineages that inhabit either freshwater or marine waters.

Another property that might help distinguish between secreted and non-secreted proteins is their cost, as extracellular proteins tend to be composed of amino acids that are less expensive to produce. We calculated the average ATP cost, average nitrogen and sulphur content for the marine dataset and found mixed results. On the one hand, there is a small (effect size 0.0121) but significant (p-value < 0.00483, kruskal–wallis test) reduction in ATP cost for both extracellular and periplasmic relative to cytoplasmic proteins. However, the magnitude of the effect is highly dependent on phylum, confirming previous reports. In the marine dataset, only Synechococcales, Flavobacteriales and Rhodospirillales show a large reduction in ATP cost (Supplementary Fig. S6, effect size > 0.2 for all three groups). This provides further support to the observation mentioned above - that differences in AAC between taxa are related to trophic strategy. With regards to nitrogen and sulphur content, there is a clear decrease of these elements in extracellular proteins compared to their intracellular counterparts (p-value < 8.29E−193, effect size = 0.1 for sulphur and 0.178 for nitrogen) (Supplementary Fig. S7). However, there is a decrease in average carbon content in extracellular proteins compared to the other groups, probably as a side effect of amphiphilic amino acids being replaced by the smaller amino acids alanine and glycine (Supplementary Fig. S7).

It has been reported that in the gut microbiota, extracellular proteins are longer than their intracellular counterparts. We confirm that this observation stands in our marine dataset (Supplementary Fig. S7). Interestingly, this increase in length is not accompanied by an increase in molecular weight, as the amino acids more prevalent in extracellular proteins tend to have smaller side chains. A possible explanation for this phenomenon is that larger proteins present less diffusion length, therefore increasing the time that the protein stays close to the bacteria. Garcia-Garcera et al. discovered that bacteria inhabiting poorly structured habitats tend to produce protein with lower diffusion lengths, as there is no community to share the burden of producing extracellular effectors. In this paper, diffusion length was calculated with the Einstein-Stokes equation, in which diffusion length is a function of temperature and molecular weight. We found no significant differences in molecular weight, but perhaps more sophisticated methods to calculate protein diffusion rates, such as HYDROPRO or HullRad might uncover differences in this property between extracellular and cytoplasmic proteins Another, more parsimonious explanation is that secreted proteins are larger since they need to incorporate sorting signals in their sequence.

Sequence order effect

A major drawback of composition-based protein descriptors is that they ignore the distribution of amino acids in the protein sequence. As it has been demonstrated that including this information in a prediction model can significantly increase its performance for cellular location prediction tasks, it is of interest to assess their possible contributions to subcellular localization prediction. The most basic descriptor of sequence order would be di/tripeptide composition, as it incorporates neighbourhood information. The top 40 dipeptides from our marine dataset that vary the most between subcellular locations are shown in Supplementary Fig. S8. As a general rule, most of these dipeptides are combinations of the single amino acids with the most variation, but it includes some unusual pairings that can help distinguish between cellular locations. An example of this are the dipeptide pairs NP and PN, which are most abundant in extracellular and periplasmic proteins.

In order to assess sequence order over larger distances, a commonly used approach is to calculate autocorrelation in the sequence (that is, how similar the sequence is to itself by given intervals), after translating each amino acid to a numerical value provided by a propensity scale, such as those provided by AAIndex. Although these descriptors are less intuitively interpretable and capture less information than 3D-based protein analyses, they still capture information not detected by amino acid composition. For example, autocorrelation plots have been found to reflect the tertiary structure of a protein and have been used to predict secondary structure composition. Furthermore, they are much faster to calculate and only require the amino acid sequence. As the effect of the chosen propensity scale is marginal at best, we decided to only use two types of autocorrelation descriptors: partial Quasi Sequence Order (pQSO), an adaptation of QSO which uses the physicochemical distances between pairs of amino acids calculated by Schneider and Wrede, and Pseudo amino acid composition (pPAAC), an adaptation of Chou’s PAAC metric which combines hydrophobicity, hydrophilicity and residue mass propensity scales. A non-parametric MANOVA analysis of the ILR-transformed features reveal that although each variable has a small effect size (relative effects mean = 0.462, std = 0.025 for pQSO, relative effects mean = 0.503, std = 0.026 for pPAAC) and a lot of variability, the combination of at least 20 autocorrelation measures is enough to reliably distinguish between subcellular locations (p-value < 1.5E−20 for both QSO and PseAAC). As the correlation between these two sets of features was small (< 0.3, Pearson correlation), we included both sets of predictors in the final model.

Machine learning model design and validation

With a validated set of protein descriptors, we tested if this information could be used to improve on current subcellular location prediction methods. xgBoost was our algorithm of choice since it presents many advantages: can be used with non-parametric data, supports multi classification, is reasonably resistant to overfitting if the right hyperparameters are used and it has been used extensively for protein classification problems. As our analysis had revealed a gradient of adaptations with the order extracellular > periplasmic > cytoplasmic, we suspected predictions would improve by framing the problem as an ordinal classification. Therefore, two classification strategies were implemented: a multiclass classifier, treating each subcellular location as an independent class, and an ordinal classifier, which is aware of the intrinsic ordering between the classes.

The results of the five-fold cross-validated testing of both Ayu implementations (ordinal and multiclass), pSORTb 3.0 and BUSCA are collected on Table 1. In general, when comparing MCC and Kappa scores, all Ayu implementations (MCC > 0.89, Kappa > 0.89 for all implementations) significantly overperform compared to pSORTb3 (MCC = 0.64, Kappa = 0.64) and BUSCA (MCC = 0.43, Kappa = 0.42). Examining the precision-recall metrics for each group (Fig. 2A), it is clear that while pSORTb3 obtains a reasonable precision in all three subcellular localizations (0.97 for cytoplasmic predictions, 0.9 for periplasmic, 0.93 for extracellular), its recall score is fairly low for periplasmic (0.58) and extracellular (0.46) proteins; meaning pSORTb3 mis-classifies more than half of the non-cytoplasmic proteins as false negatives. BUSCA only achieves good scores for cytoplasmic proteins (precision = 0.93, recall = 0.9), presenting low precision and recall scores for extracellular proteins (precision = 0.44, recall = 0.61). The reason for a worse than random precision score is probably due to the fact that BUSCA only uses signal peptide information to predict the extracellular localization, therefore including periplasmic proteins that also contain a signal peptide.

Both versions of Ayu (multiclass and ordinal) present an improvement over the other classifiers, although with differences in recall and precision. The ordinal version presents better precision scores than both multiclass versions (0.97/0.92/0.93 Ordinal, 0.97/0.91/0.89 Multiclass for Cytoplasmic/Periplasmic/Extracellular), but the multiclass version achieves better recall (0.99/0.82/0.7 Ordinal, 0.99/0.85/0.81 Multiclass for Cytoplasmic/Periplasmic/Extracellular). Likewise, the application of the SMOTE algorithm to ameliorate the imbalance between protein classes results in a small improvement in the multiclass implementation of Ayu (MCC score improvement of 0.02). As the use of SMOTE only affects training time does not impact prediction time, the SMOTE version of the multiclass implementation was kept for the final version of Ayu, as we consider the tradeoff between recall and precision to be better in the multiclass version.

As xgBoost belongs to the algorithm family of boosted trees, we are able to obtain feature importance scores, which contain information about which feature descriptors are more useful to discriminate between classes. Figure 2B shows the permuted feature importance results for the Ayu SMOTE multiclass version. “Gain” represents the contribution of each feature to the classification, while “Weight” indicates how many times the feature appears in a tree across the ensemble of trees. Therefore, the graph shows that while Signal peptide information is by far the most important feature to discriminate between the three cellular locations, Dipeptide composition is the feature that appears most in the trees, probably due to the fact that it includes more features.

However, this does not mean that all features are equally important for the classification to all locations. An example of this can be found on Supplementary Fig. S9, which plots the permuted feature importance gains for the individual binary classifiers in Ayu Ordinal. For the cytoplasmic vs non-cytoplasmic predictor, the most important protein descriptors are the presence of a signal peptide, the ratio of transmembrane to non-transmembrane regions in the protein and the dipeptide composition. On the other hand, the classifier tasked with discriminating between periplasmic and extracellular proteins finds that amino acid and dipeptide composition are more important than the presence of a signal peptide, corroborating the magnitude of the differences between proteins in different subcellular locations.

Application to real world (marine) dataset: Tara Oceans

In order to test the performance of Ayu in a real world metagenomic dataset, we applied our prediction tool on 6 Tara Oceans metagenomic and metatranscriptomic datasets, comprising three sampling sites of the prokaryotic fraction at surface and mesopelagic depths of different ocean basins (Supplementary Data S1). Out of the 46,775,154 total proteins found in the combined dataset, 73% of the sequences belong to bacterial genes, 8% to viral genes and 3% to archaeal genes, with the rest having no taxonomic classification. The results of the subcellular location classification for bacterial proteins with Ayu are summarised in Fig. 3A. Around 15.7% of the proteins were classified as transmembrane proteins by manual classification (see Methods). Out of the remaining proteins, 65.2% are classified as cytoplasmic, while 12.5% of proteins are classified as non-cytoplasmic (5.5% extracellular, 7.0% periplasmic). 5.6% of the proteins are not classified into any category. The proportion of non-integral membrane proteins secreted into the periplasm or to the extracellular milieu (12.5%) reported here falls in the range of the size of the secreted proteome reported in previous studies (10–30%). Interestingly, even though Ayu was trained only with bacterial proteins, the predictions for viral proteins are accurate, with DNA replication and auxiliary metabolic genes classified as cytoplasmic and virion structural proteins classified as periplasmic or extracellular. As virion proteins are in contact with the same extracellular milieu as bacterial extracellular proteins, this is further proof that our tool is using adaptations to the extracellular conditions to classify proteins.

As Ayu uses signal peptide information as one of its features, we are able to ascertain how many proteins predicted to each cellular location include a signal peptide that is detected by SignalP. For bacterial proteins, only 79% of periplasmic and 54.7% of extracellular proteins contain a signal peptide (Fig. 3A), reflecting the importance of GSP-independent secretion systems in the marine secretome, and the relevance of Ayu. To test the prevalence of these cases, we clustered proteins at 70% identity with 95% coverage (see Methods), presuming that proteins that cluster together should contain the same secretion signals. Our results show that out of 53,902 proteins grouped in clusters with at least 1 protein with a signal peptide, only 43,361 (80%) coded for one. These results suggest that Ayu is able to complement signal peptide prediction to recover more intra-cluster extracellular protein diversity.

The aforementioned clustering process also produced several protein clusters of predicted extracellular proteins without a signal peptide. These clusters sum up to 39,871 proteins, representing almost half of the total extracellular proteins detected in this dataset. Although only 53% of the proteins detected in this way can be annotated, it is still possible to find proteins that further prove the validity of the prediction method. Interestingly, while the extracellular proteins with a signal peptide are dominated by catabolic enzymes (peptidases, glycosyl hydrolases, sulfatases, phosphatases) and outer membrane-associated proteins (assembly factor BamB/BamE, TonB-dependent receptors), the bulk of the signal peptide-free protein fraction is composed of flagellins (FliC/FljB/LafA), flagellar basal body rod proteins (FlgC/FlaE/FlgF) and flagellar hook-associated proteins. These proteins are secreted via the T3SS pathway, which translocates proteins directly from the cytoplasm. Proteins with secretion signals for other pathways can also be found, including proteins with RTX repeats (involved in secretion via T1SS) and domains related to secretion pathways II, III, V, VI and VII. Surprisingly, some proteins contain sorting domains related to T9SS, a widespread secretion pathway in Bacteroidetes that is dependent on the GSP for translocation into the periplasm. The absence of signal peptide in these proteins suggests that either these proteins reach the periplasm via another method or the signal peptide used is not detected by SignalP 6. A group of extracellular proteins of particular interest are those related with interaction with other members of the microbial community. Examples of this group are an homolog of the protein Reb, involved in interactions between bacteria and eukaryotes; and several proteins encoding the Nif11 domain, found in the leader peptide of various bacterial microcins. These results provide further evidence for the predictions of Ayu to not be mere false positives.

Finally, we studied and compared the metagenomic to the metatranscriptomic dataset from the same Tara Oceans samples to test for differential patterns based on gene content or expression. Overall, there was a relatively high expression of genes identified as coding for secreted proteins, confirming the relevance of the secretome in the environment. We found that the ratio of metatranscriptome RPKG (Reads Per Kilobase of sequence per Gigabase of dataset) to metagenome RPKG is an order of magnitude larger in mesopelagic samples compared to surface samples (log2 fold metaT/metaG in surface = −0.93, log fold metaT/metaG in mesopelagic = 0.84) (Supplementary Fig. S10), which is consistent with previous studies reporting protein activity and transcription to increase with depth. As a general rule, metagenome samples have a wider spread of genes in the gene pool (8674 genes with > 1 RPKG in surface samples, compared to 7292 in mesopelagic), but expression of the genes present is higher in mesopelagic samples. We have found differences between taxonomic clades (Supplementary Fig. S11). For example, genes assigned to the orders Alteromonadales, Vibrionales and Flavobacteriales show the previously reported increase in mesopelagic transcription (effect size > 0.4 for all three orders), while those genes assigned to orders Pelagibacterales, Synechococcales and Bacillales show a lower increase (effect size = 0.17 for Synechococcales). These results are consistent with previous studies reporting Gammaproteobacteria and Flavobacteriales as the main producers of extracellular proteins in the mesopelagic. Finally, we found no clear distinction in gene functions between surface and mesopelagic samples (Fig. 3B), which is also consistent with previous studies reporting high functional redundancy of microbial functions with depth.

Methods
Data collection, annotation and curation

A dataset was compiled to both study the adaptations of bacterial proteins to the marine environment and train a machine learning tool. An overall flowchart of the data collection process can be found in Supplementary Fig. S12. First, an extensive bibliographic search was conducted to select a collection of bacteria that met all of the following characteristics: a) the bacteria had been isolated from the marine environment, or its genome sequence assembled from a marine sample; b) the bacteria was a mild halophile, as defined in refs. 38, 81, or had been only described in marine environment. This process resulted in a collection of 105 Gram-positive and 929 Gram-negative bacteria (Supplementary Data S2). Protein sequences were then downloaded from UniProtKB or the NCBI Genomes database.

The gold standard for protein datasets is the use of proteins with experimentally proven locations. Unfortunately, out of the 17,047 proteins coded by our collection of marine bacteria, only 38 of them have their location experimentally proven, and as we are interested in studying the adaptation to the marine environment, we classified the proteins following these criteria. First, proteins with a reviewed cellular location annotation in UniProtKB were divided into three groups (cytoplasmic, periplasmic and secreted). Moonlighting proteins (sequences reported to be present in more than one subcellular location) or proteins with ambiguous annotation were removed from the dataset.

According to previous reports, proteins that share function and have at least ≥25% global sequence identity tend to share subcellular location. Unreviewed marine proteins were aligned to a subset of prokaryotic proteins with manually curated subcellular locations downloaded from UniProtKB (“reference dataset”, 1,131,685 proteins) using ggsearch. For each marine protein with a match to a reference sequence that passed the 25% identity threshold, domains and features of both sequences were detected using InterProScan, and the cellular location from the reference sequence was propagated to the marine one only if the content and synteny of domains and features was the same for both proteins. Finally, proteins were clustered to 30% sequence identity and coverage alignment of at least 80% using CD-HIT. From each cluster, the protein with the best match to a reference sequence was selected to be part of the training dataset. The resulting dataset contains 17,047 proteins (1934 Gram-positive, 15,113 Gram-negative), divided into 14,410 cytoplasmic, 1873 periplasmic, 764 extracellular proteins.

Differences in signal peptide content is a defining characteristic of the different locations, as by definition cytoplasmic proteins will never code for one, while proteins the other two categories might. Therefore, we performed a circularity test for the dataset as defined by Riezler and Hagmann. Briefly, two Generalised Additive Models (GAMs) were fit to the dataset, one including signal peptide information as an explanatory variable (GAM_wSP) and one without (GAM_woSP). The tests show that the dataset does not meet any of the two criterions established: although the scaled deviance is remarkably higher in GAM_wSP than in GAM_woSP (D2 = 0.524 vs D2 = 0.352), it is not close to the value of near 1 that we would expect if the signal peptide information was able to predict the entire dataset. Supplementary Fig. S13 shows that features in both GAMs are still contributing to the classification of the model, as evidenced by the fact that they present a non-zero feature shape. Therefore, the second criterion (nullification of the contributions of other features) can be also ruled out.

From this compiled protein dataset, 70% of the proteins (stratified by subcellular location) were reserved for data exploration, feature extraction and model training, while the remaining 30% were only used for further evaluation of the model (15% for validation during training, 15% for testing).

In order to test if these amino acid adaptations were only found in the marine environment, the ESKAPEE group (Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter baumanii, Pseudomonas aeruginosa, Enterobacter spp. and Escherichia coli) were used as a control group. This group was chosen as they are an extensively studied group of bacteria, due to their clinical importance, and none of its members are typically found in marine open waters. Proteins from these organisms were downloaded from PSORTdb 4.0, including both experimentally-proven and predicted subcellular locations.

Statistics and compositional analysis

Standard statistical analyses (Kruskal–Wallis test, post hoc Dunn’s test, eta squared effect size calculations) were performed using the R package statix. Non-parametric MANOVA was performed using the R package npmv. As many of the feature datasets employed in this paper are compositional, that is, data that carry only relative to a total, special care has been taken to use compositionally-appropriate methods when necessary. When interpretability of the results was not a concern, the scikit-bio package was used to treat compositions first with multiplicative replacement to remove all zero values, then an isometric log ratio (ILR) transformation was applied to the composition, in order to transform the data into linearly independent, subcompositionally coherent features that can be analysed with standard statistical methods. The exception to this is the compositional data obtained from TMBed. As we were only interested in the ratio between transmembrane and extramembrane components, the composition was modelled as an amalgam logratio with the following formula (1):

SLR = log ((pH + pB) / (pS + po + pi)) (1)

In which pH, pB, pS, po, pi represent the proportions for transmembrane alpha helices, beta barrels, signal peptide, cytoplasmic and external regions respectively. As ILR-transformed data is significantly more difficult to interpret, we used non-transformative methods for data exploration. The weighted logratio variance plot was calculated using the easyCODA package, while regression analysis was performed using the DirichletReg R package. Full regression results can be found in Supplementary Text S1.

For the circularity check, the methodology presented in ref. 89 was applied. Generalized Additive Models (GAMs), feature shape plots and scaled deviance calculation were performed using the mgcv R package.

Protein feature extraction

A list of the protein descriptors used in this paper can be found in Supplementary Data S3. Each protein descriptor was chosen by their capacity to discriminate between proteins located in different cellular compartments (see Results). Each protein is represented by a 466 length feature vector. All features were extracted using in-house scripts, following the methodology stated in their respective papers.

Many of the protein features that account for sequence order include AAC in their feature list, introducing redundancy in the feature set. With the objective of including multiple sequence order features and removing redundant information from the feature dataset, slightly modified formulas for quasi sequence order (QSO) and pseudo amino acid composition (PseAAC) were employed. Descriptions of these modified formulas can be found in Supplementary Text S2. Transmembrane regions, signal peptide information and isoelectric points for each protein were also calculated and included as features, using the programs TMBed, SignalP 6.0 and IPC 2.0 respectively.

Model training, optimization and validation

The features extracted in previous steps were used to train xgBoost models using the python package xgBoost. Both a multiclass and an ordinal classifier were trained, with the latter being implemented following the scheme described in ref. 101. Briefly, N-1 binary classifiers are trained to predict N categories, and the probability for prediction of the middle categories is defined by the ensemble of the N-1 binary classifiers. For this specific application, two classifiers were trained, one to discriminate between cytoplasmic and periplasmic+extracellular proteins and another to distinguish between periplasmic and extracellular proteins (Supplementary Fig. S14).

Hyperparameters for each model were optimised with a five-fold cross validation grid search, as implemented in scikit-learn. Early stopping of the fitting step was implemented to reduce the risk of overfitting. The train-test partition and the splits for cross-validation were performed using graphpart, an homology-aware partitioner, using the recommended 30% identity threshold. As stated previously, the uniprot training dataset is severely imbalanced for two of the three classes (1:1:10 ratio extracellular:periplasmic:cytoplasmic). As imbalanced ratios between classes might cause difficulties during the fitting process, an oversampling SMOTE algorithm, as implemented in the python package imbalanced learn, was applied to the training split before fitting the multiclass classifier. Feature importance analysis was also performed for all models to assess the performance of each feature group. Finally, fitted models were evaluated against the validation partition, using the metrics precision (Pr), recall (Rc), precision-recall curves, Matthews correlation coefficient (MCC), using Gorodkin’s k-category definition and Cohen’s Kappa (Kappa).

Comparison against other cellular localization classifiers

The performance of our tool was compared against other subcellular localization predictors (pSORTb 3.0 and BUSCA). The marine testing dataset was predicted using pSORTb 3.0 and BUSCA, as they are the two non-ensemble predictors that have the best reported performance to date and are readily available. As BUSCA is only available as a web server with a limited throughput, an in-house script was used to run predictions. Performance metrics mentioned in the previous section were calculated for both methods per class using scikit-learn. Precision-recall curves were only calculated for pSORTb, as BUSCA only provides probability scores for the predicted category.

Sublocalization prediction in Tara Oceans samples

To test the performance of the machine learning model in a real marine dataset, the proteins contained in Ocean Microbial Reference Catalogue v2 (OM-RGC.v2) were downloaded and classified into different cellular locations using the machine learning models built in previous steps (See “Model training, optimization and validation”). The dataset was functionally annotated against the CDD database profile database using hmmscan. A profile hit was kept if the e-value < 1e−5 and the alignment covered at least 75% of the profile. Proteins from OM-RGC.v2 were clustered into protein clusters to 70% identity and 95% coverage using mmseqs2. Membrane proteins were identified by calculating transmembrane regions for all proteins with TMBed, classifying a protein as transmembrane if it contained at least 5% of its protein length in a beta barrel or a transmembrane helix. Proteins with only one transmembrane helix on the N-Terminal were included as part of the membrane location only if they contained a domain with a GO term with located the protein to a membrane, as signal peptides can be misclassified as transmembrane helices. Finally, recruitment values in metagenome and metatranscriptome samples provided in the OM-RGC.v2 dataset were used to compare gene abundance and gene transcription rates between protein subcellular locations and depth.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "General comment:\nThis manuscript describes a novel tool for predicting secreted proteins in marine bacteria. It is a very welcome addition to the\nprokaryotic protein sorting prediction toolbox, where the reference method, PSORTb 3.0, is more than a decade old.\nHowever, I have serious doubts about the separation between training and test data and hence about the reported predictive\nperformance.\nIt is stated that 30% of the dataset were only used for evaluation of the model, and that \"proteins were grouped into clusters\nwith 50% sequence identity and at a coverage alignment of at least 80%, using the MMseqs2 software suite\".\nThere are two problems with this. First, 50% identity and 80% coverage constitutes a rather high threshold; if this is\ncompared to other recent subcellular location prediction methods (mainly in eukaryotes), 30% or 40% identity would be more\nappropriate. Second, MMseqs2 simply cannot be used in this way. Only the representative sequence of each cluster (as\ndefined by MMseqs2) have a maximum identity of 50% to other representative sequences; this restriction does not apply to\nother members of the clusters. In other words, there may very well be sequences in the test set that are considerably more\nthan 50% identical to sequences in the training set.\nA proper separation between training and test can be achieved by either (i) reducing the dataset by only keeping the\nrepresentative sequence of each MMseqs2 cluster, or (ii) using a homology partitioning scheme such as GraphPart\n(https://doi.org/10.1093/nargab/lqad088).\nSpecific comments, major:\n- l.88-90: \"can only be run via a docker container, which severely limits the throughput of proteins it is able to process\": No, a\ndocker container does not in itself severely limit the throughput. What is limiting the throughput of PSORTb 3.0 is its\ndependence on BLAST.\n- l.160-164: The sentence beginning with \"When combined with...\" is difficult to understand. What is the \"control group\"\nhere?\n- l.232-234: \"a small ... but significantly lower ... reduction\": did you mean \"a small but significant reduction\"? If not: what is\nthe reduction lower than?\n- l.248-249: \"this increase in size is not accompanied by an increase in molecular weight\": I don't understand this. There are\nmore amino acids, but not a higher molecular weight? Is the increase in number of amino acids compensated by a higher\nproportion of small side chains?\n- l.355-358: According to your own definition in the Abstract, the secretome comprises proteins secreted into the\nenvironment. In other words, periplasmic proteins should not be counted in. But then, it is not true that your estimate of the\nproportion of the secretome (8%) falls within the range reported in previous studies (10-30%).\n- l.369-372: \"There are two possible explanations...\" There is a third possible explanation: that the prediction by Ayu is a\nfalse positive.\n- l.403-411, the results depicted in Figure S9: I fail to see the relevance of this to subcellular location - the metric is higher in\nmesopelagic than in surface samples, but the trend is the same for all three locations?\n- l.502-504: \"As the number of...\": some details about \"domain annotation and synteny\" would be nice.\n- l.582 (and other places): You report Matthews correlation coefficient (MCC), but this measure is only defined in the two-\ncategory case, and you have three. Did you use the k-category definition described by Gorodkin\n(https://doi.org/10.1016/j.compbiolchem.2004.09.006)? In that case, cite that paper. If not, describe how you calculated one\nnumber for three categories.\n- l.605-606: \"Membrane-attached proteins were identified by annotation\": which annotation?\n- The selection of the dataset, as described in Figure S11: A sequence is only classified as \"Membrane\" if it has >3 predicted\ntransmembrane domains. Doesn't that leave a number of transmembrane proteins with 1, 2 or 3 TM helices in the dataset?\nSpecific comments, minor:\n- Abstract l.30: \"The application of this powerful novel to open ocean samples\": a word is missing, such as \"tool\" or \"method\"\n- l.71: \"shotgun metaproteomic have been used to study\" it seems there is a word missing, maybe \"techniques\"\n- l.84: \"focusing only on one prediction system\": should probably have been \"secretion system\"\n- l.98: \"adapter\" should be \"adapted\"\n- l.99: \"Salt is\" should be \"Salt has\"\n- l.104: \"osmorregulated\" should be \"osmoregulated\"\n- l.111: \"isoelectric\" missing \"point\"\n- l.150: \"cell wall type\": please be more precise - do you mean the Gram-positive / Gram-negative distinction?\n- l.159; \"was distinct for periplasmic than for extracellular proteins\": this does not read well, please rephrase\n- l.178-179: \"tolerance to temperature\": high or low temperature?\n- l.185: abbreviation \"LRA\" has not been defined\n- l.186: \"phylums\" should be \"phyla\"\n- l.189-190: \"Figure S4 shows...\": No, these are the orders shown in Figure S5. S4 shows some other orders.\n- l.236-237: You mention three orders, but Figure S6 shows families. And I see no Rhodospirillales in Figure S6.\n- l.238: \"above mentioned\" should be \"mentioned above\"\n- l.240: \"compounds\" should be \"elements\".\n- l.249: \"phenomena\" should be \"phenomenon\"\n- l.256-257: \"but perhaps more...\": the end of this sentence is incomplete.\n- l.265 and 3 other places: I would advise against using the term \"sublocation\" for bacteria, but instead simply say \"location\".\nIn my book, \"sublocation\" is a sub-category of a location, such as found in eukaryotes, where \"nucleus\" may be subdivided\ninto \"nuclear matrix\", \"nucleolus\", etc.\n- l.446: a word missing in this sentence, probably \"which\".\n- l.487: \"UniRefKB\" should be \"UniProtKB\".\n- l.506: \"excluded all other\": missing \"from\"\n- l.543-544: \"Full regression results can be found in Full regression results can be found in\": repetition.\n- l.554: \"ACC\" should probably have been \"AAC\".\n- Figure S1: the different shapes of the points are not reflected in the legend.\n- Figure S3 caption: \"accordin\" should be \"according\"\n- Supplementary Text 2: \"described by Chou in [ref]\": reference missing.\nGitHub page:\n- \"Once all protein features have been processed, protein sublocations can be predicted with the command:\" is followed by\nan empty line.", "reviewer_2": "The paper presents an analysis of the composition of extracellular proteins from marine bacteria and a new predictor to\ndiscriminate them starting from protein sequences of marine prokaryotes.\nThe research topic is relevant and important and the presentation is generally clear.\nMy main concerns regard the development and validation of the prediction method:\n- Training dataset. The annotation of the dataset is only partially derived from experiments. It is unclear how many proteins\nare experimentally annotated and how many are assigned to a location through the predictive workflow reported in figure S1\nand the following manual annotation. Rules of manual annotation should be further detailed and it should be specified how the\nassigment is performed when contrasting annotations are present for proteins in the same cluster (if this is the case). It is\nunclear whether a representative per cluster is selected and, in case, on which basis. Also, I wonder why the number of\nproteins from Gram+ and Gram- bacteria are quite similar despite the fact that selected Gram- organisms are largely higher.\n- Splitting among training and validation sets. The splitting criteria are not clear. I can assume it has been randomly\nperformed and in this case the complete independence between the sets is not granted, since the starting dataset contains\nproteins up to 50% identical (if a representative per cluster is selected). To ensure an unbiased evaluation of the method\nperformance and generalization, proteins in the validation set should share a much lower similarity with proteins in the\ntraining set. Possibly, the results should be reported only on proteins that have been experimentally annotated, to avoid\ncircularity between the initial annotation and the prediction (see next point).\n- Input features. One of the input features (SignalP6 prediction) is also used to annotate the training/validation sets. This\ncircularity increases the risk to unbiased evaluation of the method, if non-experimentally annotated proteins are used in the\ntesting phase.\n- Benchmark. Table 1 reports the scoring indexes of the method. It is unclear whether they are computed in cross-validation\nor they refer to the evaluation on the 30% of the dataset \"used for further evaluation of the model\". Since hyperparameters\nare set during the cross-valdation procedure, cross-validation results risk to be overfitted. The table reports MCC values that\nare routinely computed for binary classification. However in this case the discriminated classes are three and it then unclear\nwhat the score refers to. Some of the considered methods (i.e., Busca) does not separate secreted proteins between\nperiplasmic and extracellular proteins and could therefore be omitted.\nSome minor point:\n- the variance of the compositional data could be indicated to better evaluate the significance of their differences\n- the statistics discussed in section \"Differences in AAC between subcellular localizations and habitat\" should be provided\n- supplementary texts should be reformatted and clarified."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nGeneral comment:\n1.1 This manuscript describes a novel tool for predicting secreted proteins in marine\nbacteria. It is a very welcome addition to the prokaryotic protein sorting prediction\ntoolbox, where the reference method, PSORTb 3.0, is more than a decade old.\nHowever, I have serious doubts about the separation between training and test data\nand hence about the reported predictive performance.\nIt is stated that 30% of the dataset were only used for evaluation of the model, and\nthat \"proteins were grouped into clusters with 50% sequence identity and at a\ncoverage alignment of at least 80%, using the MMseqs2 software suite\".\nThere are two problems with this. First, 50% identity and 80% coverage constitutes a\nrather high threshold; if this is compared to other recent subcellular location\nprediction methods (mainly in eukaryotes), 30% or 40% identity would be more\nappropriate. Second, MMseqs2 simply cannot be used in this way. Only the\nrepresentative sequence of each cluster (as defined by MMseqs2) have a maximum\nidentity of 50% to other representative sequences; this restriction does not apply to\nother members of the clusters. In other words, there may very well be sequences in\nthe test set that are considerably more than 50% identical to sequences in the training\nset.\nA proper separation between training and test can be achieved by either (i) reducing\nthe dataset by only keeping the representative sequence of each MMseqs2 cluster, or\n(ii) using a homology partitioning scheme such as GraphPart\n(https://doi.org/10.1093/nargab/lqad088).\nWe thank the reviewer for their suggestions. We have implemented the reviewer’s remarks\nas part of our revision to the training dataset: First, clustering is now performed with CD-HIT\ninstead of MMSeqs, since the training dataset is now smaller, clustering identity is now 30%,\nfollowing the work of other predictor methods; and only one representative sequence is\nchosen from each cluster. Second, we have implemented homology-aware partitioning via\ngraphpart as part of the machine learning pipeline.\nSpecific comments, major:\n1.2 - l.88-90: \"can only be run via a docker container, which severely limits the\nthroughput of proteins it is able to process\": No, a docker container does not in itself\nseverely limit the throughput. What is limiting the throughput of PSORTb 3.0 is its\ndependence on BLAST.\nAmended in the text, it It now reads: “Even pSORTb 3.0, the gold standard for subcellular\nlocalization predictions, relies on homology searches against its curated profile and protein\ndatabases to obtain its predictions, which severely limits the throughput of proteins it is able\nto process.”\n1.3 - l.160-164: The sentence beginning with \"When combined with...\" is difficult to\nunderstand. What is the \"control group\" here?\nThe control group in this case refers to the non-marine ESKAPEE proteins. We have\nreformatted the paragraph so it is more intelligible. The paragraph It now reads: “ When\nconsidering both combined the changes in AAC explained by subcellular location and the\ninfluence of the marine environment, we find that in 12 amino acids (L, V, F, Y, S, T, N, Q,\nD, E, K, R) the distribution varies following their order of exposition to the environment\n(Extracellular > Periplasmic > Cytoplasm), compared to 7 in the non-marine ESKAPEE\ngroup (L, Y, G, S, T, N, R).”\n1.4 - l.232-234: \"a small ... but significantly lower ... reduction\": did you mean \"a small\nbut significant reduction\"? If not: what is the reduction lower than?\nWe meant the former. The text It now reads: “On the one hand, there is a small (effect size\n0.0161) but significant (p-value < 0.00266, kruskal-wallis test) reduction in ATP cost…”\n1.5 - l.248-249: \"this increase in size is not accompanied by an increase in molecular\nweight\": I don't understand this. There are more amino acids, but not a higher\nmolecular weight? Is the increase in number of amino acids compensated by a higher\nproportion of small side chains?\nExactly. The text has been modified to make this statement more clear, it It now reads:\n“Interestingly, this increase in size is not accompanied by an increase in molecular weight,\nas the amino acids more prevalent in extracellular proteins tend to contain smaller side\nchains.”\n1.6 - l.355-358: According to your own definition in the Abstract, the secretome\ncomprises proteins secreted into the environment. In other words, periplasmic\nproteins should not be counted in. But then, it is not true that your estimate of the\nproportion of the secretome (8%) falls within the range reported in previous studies\n(10-30%).\nThe definition of secretome varies significantly between authors and fields of study. We\nchose to use the definition that has been used for the study of extracellular enzymes in the\nmarine environment, as we believe it better represents the interests of other researchers in\nthe field and the data we used for the environmental/field application part of this study.\nHowever, the reviewer is correct in pointing out the incongruence in definitions between\nabstract and results. As there are no available experimental studies that determine the size\nof the secretome as we define it, we have decided to alter the text in the results to reflect this\nfact. Lines 375-376 now read: “The proportion of non-integral membrane proteins secreted\ninto the periplasm or to the extracellular milieu (15.5%) reported here falls in the range of the\nsize of the secreted proteome size reported in previous studies (10-30%)”.\n1.7 - l.369-372: \"There are two possible explanations...\" There is a third possible\nexplanation: that the prediction by Ayu is a false positive.\nIt would be a possibility; however, due to the fact that many of these proteins code for\ndomains related to extracellular activity(either because of their function or because they are\na sorting signal), we are inclined to believe that they are genuine. Nonetheless, we have\nremoved this statement from the paper, as it does not add anything relevant to the results.\n1.8 - l.403-411, the results depicted in Figure S9: I fail to see the relevance of this to\nsubcellular location - the metric is higher in mesopelagic than in surface samples, but\nthe trend is the same for all three locations?\nThe relevance for this plot is that the trend is not between cellular locations but between\ndepths (i.e., there is more production of extracellular proteins at lower depths, as it has been\nreported previously). We have adapted the text so that the statement and its relevance is\nclearer.\n1.9 - l.502-504: \"As the number of...\": some details about \"domain annotation and\nsynteny\" would be nice.\nWe have completely overhauled the methodology to produce the training set. In this version,\nclassification into different versions is not based on specific domains, but on homology to\nproteins with manually curated subcellular location and a complete match of domain /\nfeatures content and synteny to the aforementioned proteins.\n1.10 - l.582 (and other places): You report Matthews correlation coefficient (MCC), but\nthis measure is only defined in the two-category case, and you have three. Did you\nuse the k-category definition described by Gorodkin\n(https://doi.org/10.1016/j.compbiolchem.2004.09.006)? In that case, cite that paper. If\nnot, describe how you calculated one number for three categories.\nWe are indeed using the k-category definition. A citation for the paper has now been added\nto the manuscript.\n1.11 - l.605-606: \"Membrane-attached proteins were identified by annotation\": which\nannotation?\nIn the previous version of the manuscript, we detected Gram-positive cell wall-attached\nproteins by the presence of a LPXTG domain at the C-T. However, this step was not\nperformed as part of the analysis, since non-transmembrane proteins attached to the\nmembrane are now predicted by Ayu as being located in their respective cellular\ncompartment. We have removed all mentions of this step.\n1.12 - The selection of the dataset, as described in Figure S11: A sequence is only\nclassified as \"Membrane\" if it has >3 predicted transmembrane domains. Doesn't that\nleave a number of transmembrane proteins with 1, 2 or 3 TM helices in the dataset?\nWe are wary of including proteins that have only one Transmembrane helix, especially if it is\nfound in the C-T (because they tend to be confused with a signal peptide). We have\nchanged our classification methodology to include these proteins if they code for a domain\nthat locates them to the membrane.\nSpecific comments, minor:\n1.13 - Abstract l.30: \"The application of this powerful novel to open ocean samples\": a\nword is missing, such as \"tool\" or \"method\"\nAmended in the text. It now reads: “The application of this powerful novel tool to open ocean\nsamples…”\n1.14 - l.71: \"shotgun metaproteomic have been used to study\" it seems there is a word\nmissing, maybe \"techniques\"\nAmended in the text. It now reads: “although shotgun metaproteomic approaches have been\nused…”\n1.15 - l.84: \"focusing only on one prediction system\": should probably have been\n\"secretion system\"\nAmended in the text accordingly.\n1.16 - l.98: \"adapter\" should be \"adapted\"\nAmended in the text accordingly.\n1.17 - l.99: \"Salt is\" should be \"Salt has\"\nAmended in the text accordingly.\n1.18 - l.104: \"osmorregulated\" should be \"osmoregulated\"\nAmended in the text accordingly.\n1.19 - l.111: \"isoelectric\" missing \"point\"\nAmended in the text accordingly.\n1.20 - l.150: \"cell wall type\": please be more precise - do you mean the Gram-positive /\nGram-negative distinction?\nWe indeed mean the distinction between gram-positive / gram-negative. The sentence has\nbeen edited and it now reads: “adjusting for subcellular localization and cell wall structure\n(that is, if the proteins were coded by a gram positive or gram negative bacteria).”\n1.21 - l.159; \"was distinct for periplasmic than for extracellular proteins\": this does not\nread well, please rephrase\nThe sentence has been edited and it now reads: “Additionally, the marine environment\naffected periplasmic and extracellular proteins differently, as the AAC of extracellular\nproteins was influenced more strongly towards being stronger in polar, aromatic and charged\namino acids.”\n1.22 - l.178-179: \"tolerance to temperature\": high or low temperature?\nBoth high and low temperatures. Edited in the text, it now reads: “there are multiple reports\nwhere the addition of aromatic residues to a marine protein increased its tolerance to high\nand low temperatures”.\n1.23 - l.185: abbreviation \"LRA\" has not been defined\nAmended in the text accordingly.\n1.24 - l.186: \"phylums\" should be \"phyla\"\nAmended in the text accordingly.\n1.25 - l.189-190: \"Figure S4 shows...\": No, these are the orders shown in Figure S5. S4\nshows some other orders.\nFigure numbers have been adapted to match the description in the text.\n1.26 - l.236-237: You mention three orders, but Figure S6 shows families. And I see no\nRhodospirillales in Figure S6.\nThe mismatch between figure and text was due to the fact the previous version of Figure S6\nwas submitted for revision. Figure S6 has been replaced with the intended version.\n1.27 - l.238: \"above mentioned\" should be \"mentioned above\"\nAmended in the text accordingly.\n1.28 - l.240: \"compounds\" should be \"elements\".\nAmended in the text accordingly.\n1.29 - l.249: \"phenomena\" should be \"phenomenon\"\nAmended in the text accordingly.\n1.30 - l.256-257: \"but perhaps more...\": the end of this sentence is incomplete.\nAmended in the text. It now reads: “but perhaps more sophisticated methods, such as\nHYDROPRO51 or HullRad52 might be able to detect differences in diffusion rates between\nintra and extracellular proteins.”\n1.31 - l.265 and 3 other places: I would advise against using the term \"sublocation\" for\nbacteria, but instead simply say \"location\". In my book, \"sublocation\" is a\nsub-category of a location, such as found in eukaryotes, where \"nucleus\" may be\nsubdivided into \"nuclear matrix\", \"nucleolus\", etc.\nFollowing the advice of the reviewer we have substituted all mentions of “sublocation” to\n“cellular location”.\n1.32 - l.446: a word missing in this sentence, probably \"which\".\nAmended in the text accordingly.\n1.33 - l.487: \"UniRefKB\" should be \"UniProtKB\".\nAmended in the text accordingly.\n1.34 - l.506: \"excluded all other\": missing \"from\"\nAmended in the text accordingly.\n1.35 - l.543-544: \"Full regression results can be found in Full regression results can be\nfound in\": repetition.\nAmended in the text accordingly.\n1.36 - l.554: \"ACC\" should probably have been \"AAC\".\nAmended in the text accordingly.\n1.37 - Figure S1: the different shapes of the points are not reflected in the legend.\nAmended in Figure S1.\n1.38 - Figure S3 caption: \"accordin\" should be \"according\"\nAmended in Figure S3.\n1.39 - Supplementary Text 2: \"described by Chou in [ref]\": reference missing.\nAmended in the text.\nGitHub page:\n1.40 - \"Once all protein features have been processed, protein sublocations can be\npredicted with the command:\" is followed by an empty line.\nAmended in the GitHub page.\nReviewer #2 (Remarks to the Author):\n2.1 The paper presents an analysis of the composition of extracellular proteins from\nmarine bacteria and a new predictor to discriminate them starting from protein\nsequences of marine prokaryotes.\nThe research topic is relevant and important and the presentation is generally clear.\nMy main concerns regard the development and validation of the prediction method:\nTraining dataset. The annotation of the dataset is only partially derived from\nexperiments. It is unclear how many proteins are experimentally annotated and how\nmany are assigned to a location through the predictive workflow reported in figure S1\nand the following manual annotation. Rules of manual annotation should be further\ndetailed and it should be specified how the assigment is performed when contrasting\nannotations are present for proteins in the same cluster (if this is the case). It is\nunclear whether a representative per cluster is selected and, in case, on which basis.\nAlso, I wonder why the number of proteins from Gram+ and Gram- bacteria are quite\nsimilar despite the fact that selected Gram- organisms are largely higher.\nThe reviewer raises an excellent point about the need for experimentally proven datasets.\nAlthough we agree that using a validation set with experimentally proven proteins would be\ndesirable, we are limited to the study of marine prokaryotes in order to assess the effect of\nsalt adaptation, and there are only around 40 proteins with an experimentally proven cellular\nlocation in our collection of organisms. Yet, we have now improved our manual annotation\nprocess so that it is more specific and more strict at all steps of the process (See “Data\ncollection, annotation and curation” in Materials and Methods of the current revision). With\nregards to the question about the number of proteins in Gram + and Gram - bacteria, they\nwere quite similar because we did not cluster Gram + and Gram - proteins together. This has\nnow changed in the new version, and the new protein numbers reflect that (2346 Gram +,\n15682 Gram -).\n2.2 - Splitting among training and validation sets. The splitting criteria are not clear. I\ncan assume it has been randomly performed and in this case the complete\nindependence between the sets is not granted, since the starting dataset contains\nproteins up to 50% identical (if a representative per cluster is selected). To ensure an\nunbiased evaluation of the method performance and generalization, proteins in the\nvalidation set should share a much lower similarity with proteins in the training set.\nPossibly, the results should be reported only on proteins that have been\nexperimentally annotated, to avoid circularity between the initial annotation and the\nprediction (see next point).\nThe reviewer raises two excellent points with regards to the homology between training /\nvalidation sets and the value of a test set composed of experimentally annotated proteins.\nThe first point was also raised by reviewer #1 and is answered in detail in question 1.1. As\nfor the second point, although we agree that using a validation set with experimentally\nproven proteins would be desirable, we are limited to the study of marine prokaryotes in\norder to assess the effect of salt adaptation, and there are only around 40 proteins with an\nexperimentally proven cellular location in our collection of organisms.\n2.3 - Input features. One of the input features (SignalP6 prediction) is also used to\nannotate the training/validation sets. This circularity increases the risk to unbiased\nevaluation of the method, if non-experimentally annotated proteins are used in the\ntesting phase.\nWe thank the reviewer for their suggestion. Although the new method for annotating the\ntraining/validation sets does not employ SignalP6 prediction, we have applied a test for\ncircularity to confirm that this was not an issue. The outcome of this test was successful, and\nthe results and description of the test can be found in the Materials & Methods section.\n2.4 - Benchmark. Table 1 reports the scoring indexes of the method. It is unclear\nwhether they are computed in cross-validation or they refer to the evaluation on the\n30% of the dataset \"used for further evaluation of the model\". Since hyperparameters\nare set during the cross-valdation procedure, cross-validation results risk to be\noverfitted. The table reports MCC values that are routinely computed for binary\nclassification. However in this case the discriminated classes are three and it then\nunclear what the score refers to. Some of the considered methods (i.e., Busca) does\nnot separate secreted proteins between periplasmic and extracellular proteins and\ncould therefore be omitted.\nThe scoring indexes are calculated during testing with the validation dataset, only using the\nhyperparameter tuning and training of the models. The MCC used in this paper is Gorodkin’s\nk-category definition, which can be used for multi classification. The text has been updated\nto reflect this fact, as stated in comment 1.10. Regarding the inclusion of BUSCA, although it\nis true that it does not distinguish between periplasmic and extracellular proteins for\nGram-positive bacteria, we consider that it is important to provide as many comparisons to\nother state-of-the art classifiers as possible, and it is one of the only classifiers that is readily\navailable.\n2.5 - the variance of the compositional data could be indicated to better evaluate the\nsignificance of their differences.\nThe statistics are provided in Supplementary Text S1. A sentence has been added to the text\nto better direct readers to that section of the paper.\n2.6 - the statistics discussed in section \"Differences in AAC between subcellular\nlocalizations and habitat\" should be provided.\nSee question 2.5.\n2.7 - supplementary texts should be reformatted and clarified.\nThe Supplementary texts have been reformatted for clarity."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I am pleased to see that the authors have responded appropriately to almost all my comments from the first review. However,\nthe proofreading could have been better. I still have a number of minor comments.\nThe process of reviewing the revised manuscript was made unnecessarily complicated by the fact that the authors had not\nmarked the changes in the manuscript. Furthermore, there are no line numbers in the revised manuscript, making it difficult\nto refer to specific places.\nSpecific comments:\nThe one comment I had that is not addressed in the revision is 1.40: '\"Once all protein features have been processed, protein sublocations can be predicted with the command:\" is followed by an empty line.'\nThe authors wrote that this was amended in the GitHub page, but it is not. This should arguably be the most important\ncommand in the README file, and the field is empty.\nTable 1: Caption missing!\np.4: \"ESKAPEE pathogens\" should be explained the first time it is used.\np.4: \"by a gram positive or gram negative bacteria\": should be \"bacterium\" in singular.\np.4: \"p-value 2.2e-16\": equal to or less than?\np.4: \"Amino acids E, G, I, K, P, T, Vare especially\": space missing between \"V\" and \"are\".\np.6: \"Interestingly, this increase in size is not...\": I would suggest \"length\" instead of \"size\".\np.6: \"but perhaps more sophisticated methods, such as HYDROPRO51 or HullRad52.\": This sentence is incomplete.\np.7: BUSCA should be cited the first time it is mentioned.\np.7: \"significantly over perform compared to\": \"overperform\" should be one word.\np.8: \"periplasmic (0.58) extracellular (0.46)\": an \"and\" is missing.\np.8: \"he use of SMOTE\": something is missing, perhaps it should have been \"Since the use of SMOTE\"\np.9: SignalP6 should be cited the first time it is mentioned.\np.10: \"there was a relatively high expression of these genes\": It is not entirely clear what \"these\" refers to here.\np.11: \"out of the X proteins coded by our collection\": number missing!\np.12: \"domains and features both sequences\": an \"of\" is missing.\np.12: \"as defined in [ref].\": citation missing.\np.12 : \"Figure S12 shows that\": should be \"Figure S13\"\np.12: \"the fact that their feature shapes have a non-zero feature shape\": should it have been \"the fact that they have a non-\nzero feature shape\"?\np.12: \"15% for testing during training, 15% for validation\": I would suggest swapping the words \"testing\" and \"validation\" to\nconform with standard usage of these words in the field, see e.g.\nhttps://en.wikipedia.org/wiki/Training,_validation,_and_test_data_sets\np.13: \"using the R package statix89\": period missing.\np.13: TMBed should be cited the first time it is mentioned.\np.13: \"Generalized Applied Models (GAMs)\": \"Applied\" should be \"Additive\".\np.14: \"performed using graphpart 102, an homology-aware partitioner\": using which threshold?\nFigure S7 caption: \"sulphut\" should be \"sulphur\"\nFigure S14: swap \"testing\" and \"validation\", see my last comment to p.12.", "reviewer_2": "Authors satisfacorily replied to my concerns."}, "author_response": "RESPONSE TO REVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nI am pleased to see that the authors have responded appropriately to almost all my\ncomments from the first review. However, the proofreading could have been better. I\nstill have a number of minor comments.\nThe process of reviewing the revised manuscript was made unnecessarily\ncomplicated by the fact that the authors had not marked the changes in the\nmanuscript. Furthermore, there are no line numbers in the revised manuscript,\nmaking it difficult to refer to specific places.\nSpecific comments:\nThe one comment I had that is not addressed in the revision is 1.40: '\"Once all protein features have been processed, protein sublocations can be predicted with the\ncommand:\" is followed by an empty line.'\nThe authors wrote that this was amended in the GitHub page, but it is not. This should\narguably be the most important command in the README file, and the field is empty.\nWe thank the reviewer for their kind comments and we apologize for any inconveniences\ncaused by the manuscript (unmarked edits, no line numbers). The README file in the\nGitHub page has been extensively improved, going further from the changes mentioned in\n1.40.\n1.1: Table 1: Caption missing!\nAdded in table.\n1.2: p.4: \"ESKAPEE pathogens\" should be explained the first time it is used.\nAmended. Line 155 now reads: \"As our control group, we chose the ESKAPEE pathogens\n(Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter\nbaumanii, Pseudomonas aeruginosa, Enterobacter spp. and Escherichia coli), as they are\nan extensively studied group of bacteria37, and none of their members are commonly found\nin marine environments.\"\n1.3: p.4: \"by a gram positive or gram negative bacteria\": should be \"bacterium\" in\nsingular.\nAmended in text.\n1.4: p.4: \"p-value 2.2e-16\": equal to or less than?\nLess than 2.2e-16. Amended in text.\n1.5: p.4: \"Amino acids E, G, I, K, P, T, Vare especially\": space missing between \"V\"\nand \"are\".\nAmended in text.\n1.6: .6: \"Interestingly, this increase in size is not...\": I would suggest \"length\" instead\nof \"size\".\nAmended in text.\n1.7: p.6: \"but perhaps more sophisticated methods, such as HYDROPRO51 or\nHullRad52.\": This sentence is incomplete.\nAmended in text. Line 278 now reads: \"We found no significant differences in molecular\nweight, but perhaps more sophisticated methods to calculate protein diffusion rates, such as\nHYDROPRO52 or HullRad53 might uncover differences in this property between extracellular\nand cytoplasmic proteins.\"\n1.8: p.7: BUSCA should be cited the first time it is mentioned.\nAmended in text.\n1.9: p.7: \"significantly over perform compared to\": \"overperform\" should be one word.\nAmended in text.\n1.10: p.8: \"periplasmic (0.58) extracellular (0.46)\": an \"and\" is missing.\nAmended in text.\n1.11: p.8: \"he use of SMOTE\": something is missing, perhaps it should have been\n\"Since the use of SMOTE\"\nAmended in text. Line 357 now reads: \"As the use of SMOTE only affects training time does\nnot impact prediction running time, the SMOTE version of the multiclass implementation\nwas kept for the final version of Ayu.\"\n1.12: p.9: SignalP6 should be cited the first time it is mentioned.\nAmended in text.\n1.13: p.10: \"there was a relatively high expression of these genes\": It is not entirely\nclear what \"these\" refers to here.\n\"These\" refers to extracellular proteins. Line 444 now reads \"Overall, there was a relatively\nhigh expression of these genes identified as coding for secreted proteins...\" in order to\nimprove readability.\n1.14: p.11: \"out of the X proteins coded by our collection\": number missing!\nAmended in text.\n1.15: p.12: \"domains and features both sequences\": an \"of\" is missing.\nAmended in text.\n1.16: p.12: \"as defined in [ref].\": citation missing.\nAmended in text.\n1.17: p.12 : \"Figure S12 shows that\": should be \"Figure S13\"\nAmended in text.\n1.18: p.12: \"the fact that their feature shapes have a non-zero feature shape\": should it\nhave been \"the fact that they have a non-zero feature shape\"?\nIndeed. Amended in text.\n1.19: p.12: \"15% for testing during training, 15% for validation\": I would suggest\nswapping the words \"testing\" and \"validation\" to conform with standard usage of\nthese words in the field, see e.g.\nhttps://en.wikipedia.org/wiki/Training,_validation,_and_test_data_sets\nAmended in text.\n1.20: p.13: \"using the R package statix89\": period missing.\nAmended in text.\n1.21: p.13: TMBed should be cited the first time it is mentioned.\nAmended in text.\n1.22: p.13: \"Generalized Applied Models (GAMs)\": \"Applied\" should be \"Additive\".\nAmended in text.\n1.23: p.14: \"performed using graphpart 102, an homology-aware partitioner\": using\nwhich threshold?\nUsing the recommended 30% identity threshold. Line 635 now reads \" The train-test partition\nand the splits for cross-validation were performed using graphpart103, an homology-aware\npartitioner, using the recommended 30% identity threshold.\"\n1.24: Figure S7 caption: \"sulphut\" should be \"sulphur\"\nAmended in figure.\n1.25: Figure S14: swap \"testing\" and \"validation\", see my last comment to p.12.\nAmended in figure.\nReviewer #2 (Remarks to the Author):\nAuthors satisfacorily replied to my concerns."}]

Decision letter