Macrel: antimicrobial peptide screening in genomes and metagenomes

Macrel classifiers

Features

Local features, those dependent on the order of the peptide sequence, were inspired by the composition-transition-distribution (CTD) framework Dubchak et al. (1995, 1999). Physicochemical properties of a peptide at its N terminal are informative for the prediction of its antimicrobial activity (Bahar & Ren, 2013; Bhadra et al., 2018). Therefore, we defined features based on the normalized position of the first amino acid in a group of interest.

Global features, which are independent of amino acids primary sequence, were chosen to capture well-described AMP characteristics, such as the typical AMPs composition of approximately 50% hydrophobic residues, usual positive charge and folding into amphiphilic ordered structures (Zhang & Gallo, 2016). The mechanism of antimicrobial activity also was summarized in Macrel’s features by global descriptors of stability, amphiphilicity and predisposition of a peptide to bind to membranes.

Therefore, Macrel combines 6 local and 16 global features (see Table S1), grouped as:

1. A new local feature group (3 local features), defined as the relative position of the first occurrence of residues in three groups of amino acids defined by their free energy of transition in a peptide from a random coil in aqueous environment to an organized helical structure in a lipid phase FET—see Fig. 2. The groups are: (1, lowest FET): ILVWAMGT, (2, intermediate): FYSQCN, (3, highest): PHKEDR (Von Heijne & Blomberg, 1979).

2. Solvent Accessibility (3 local features), obtained by the distribution at first occurrence of residues organized in groups by solvent accessibility as described by Bhadra et al. (2018), using the groups: (1, buried): ALFCGIVW, (2, exposed): RKQEND, and (3, intermediate): MSPTHY.

3. Amino acid composition (9 global features) as the fraction of amino acids in groups defined by their size (area/volume), polarity, charge and R-groups: acidic, basic, polar, non-polar, aliphatic, aromatic, charged, small, tiny (Jhong et al., 2019; Nagarajan et al., 2019).

4. Charge and solubility (2 global features): peptide charge (Ebenhan et al., 2014; Chung et al., 2020) and isoelectric point (Fan et al., 2016; Wenzel et al., 2014; Chung et al., 2020).

5. Indexes for multiple purposes (3 global features): instability, aliphaticity, propensity to bind to membranes (Boman (Jhong et al., 2019; Chung et al., 2020; Boman, 2003)).

6. Hydrophobicity (2 global features): hydrophobicity (KyteDoolittle scale) and hydrophobic moment at 100° to capture the helix momentum (Ebenhan et al., 2014; Dathe et al., 1997).

Benchmarking

Implementation and availability

Macrel is implemented in Python 3 and R (R Core Team, 2018). Peptides (Osorio, Rondon-Villarreal & Torres, 2015) is used for computing features, and the classification is performed with scikit-learn (Pedregosa et al., 2011). For ease of installation, we made available a bioconda package (Grüning et al., 2018). The source code for Macrel is archived at DOI 10.5281/zenodo.3608055 (with the specific version tested in this manuscript being available as DOI 10.5281/zenodo.3712125).

The complete set of scripts used to benchmark Macrel is available at https://github.com/BigDataBiology/macrel2020benchmark and the newly simulated generated dataset of different sequencing depths is available at Zenodo (DOI 10.5281/zenodo.3529860).

Macrel: (Meta)genomic AMPs classification and REtrievaL

Here, we present Macrel (for (Meta)genomic AMPs Classification and REtrievaL, see Fig. 1), a simple, yet accurate, pipeline that processes either genomes or metagenomes/metatranscriptomes and predicts AMP sequences. We test Macrel with standard benchmarks for AMP prediction as well as both simulated and real sequencing data to show that, even in the presence of large numbers of (potentially artifactual) input smORFs, Macrel still outputs only a small number of high-quality candidates.

Macrel can process metagenomes (in the form of short reads), (meta)genomic contigs, or peptides. If short reads are given as input, Macrel will preprocess and assemble them into larger contigs. Automated gene prediction then extracts smORFs from these contigs which are classified into AMPs or rejected from further processing (see Fig. 1 and “Methods”). Putative AMPs are further classified into hemolytic or non-hemolytic. Unlike other pipelines (Jhong et al., 2019), Macrel can not only quantify known sequences, but also discover novel AMPs.

Macrel is also available as a webserver at https://big-data-biology.org/software/macrel, which accepts both peptides and contig sequences, and retrieves AMPs coded by their own genes.

Novel set of protein descriptors for AMP identification

Two binary classifiers are used in Macrel: one predicts AMP activity and another the hemolytic activity (which is invoked only for putative AMPs). These are feature-based classifiers and use a set of 22 variables that capture the amphipathic nature of AMPs and their propensity to form transmembrane helices (see Table S1).

Peptide sequences can be characterized using local or global features: local features depend on the order of the amino acids, while global ones do not. Local features have been shown to be more informative when predicting AMP activity and its targets, while global features are more informative when predicting the potency of a given AMP (Bhadra et al., 2018; Fjell et al., 2009; Boone et al., 2018). Thus, Macrel combines both, including 6 local and 16 global features (see “Methods” and Table S1):

• Free energy transition (FET) (3 local features). This is a novel feature group, which was designed to capture the fact that AMPs usually fold from random coils in the polar phase to well-organized structures in lipid membranes (Nagarajan et al., 2019). Each amino acid is assigned to one of three groups of increasing free-energy change (Von Heijne & Blomberg, 1979). The three features consist of the position of the first amino acid in each group, normalized to the length of the sequence (see Fig. 2C). Earlier works had shown that the N-terminal is particularly informative for determining AMP activity (Bahar & Ren, 2013; Bhadra et al., 2018). We adopted the fractional position encoding from the more general CTD framework (Dubchak et al., 1995, 1999).

• Solvent accessibility (3 local features). Computed in the same way as the FET features, with amino acids groups representing levels of solvent accessibility.

• Amino acid composition (9 global features). As AMPs usually have biased amino acids composition (Nagarajan et al., 2019; Jhong et al., 2019), we used the fraction of amino acids falling into nine partially overlapping classes defined by charge, size, polarity, and hydrophobicity (see “Methods” and Table S1).

• Charge (1 global feature). AMPs typically contain approximately 50% hydrophobic residues (Zhang & Gallo, 2016; Malmsten, 2014; Pasupuleti, Schmidtchen & Malmsten, 2012), and their net charges are crucial to promote the peptide-induced membrane disruption (Malmsten, 2014; Pasupuleti, Schmidtchen & Malmsten, 2012; Ringstad, Schmidtchen & Malmsten, 2006).

• Membrane binding and solubility in different media (6 global features). These capture predisposition of peptides bind to membranes, and their solubility (Ebenhan et al., 2014; Dathe et al., 1997; Jhong et al., 2019).

All 22 descriptors used in Macrel are important for classification (see Fig. S1A). The fraction of acidic residues, charge, and isoelectric point were the most important variables in the hemolytic peptides classifier. Those variables tend to capture the electrostatic interaction between peptides and membranes, a key step in hemolysis. For AMP prediction, charge and the distribution parameters using FET and solvent accessibility are the most important variables. This is consistent with reports that cationic peptides (e.g., lysine-rich) show increased AMP activity (Bhadra et al., 2018; Jhong et al., 2019; Nagarajan et al., 2019).

Compared to other tools, Macrel achieves the highest specificity, albeit at lower sensitivity

To evaluate the feature set and the classifier used in the context of the pipeline as a whole, we benchmark both the classifier implemented in Macrel, built with the training set adapted from Bhadra et al. (2018) (see “Macrel Prediction Models”), which consists of 1 AMP for each 50 negative examples, and a second AMP classifier (denoted Macrel^X), which was built using the same features and methods, but using the training set from Xiao et al. (2013), which contains 770 AMPs and 2,405 non-AMPs (approximately 1:3 ratio).

Benchmark results show that the AMP classifier trained with a more balanced dataset performs better than most of alternatives considered on this balanced benchmark, with AmPEP (Bhadra et al., 2018) achieving the best results (see Table 1).

In terms of overall accuracy on this benchmark, the AMP classifier implemented in Macrel is comparable to the best methods, with different trade-offs. In particular, Macrel achieves the highest precision and specificity at the cost of lower sensitivity. Although we do not possess good estimates of the proportion of AMPs in the smORFs predicted from real genomes (or metagenomes), we expect it to be much closer to 1:50 than to 1:3. Therefore, we chose to use the higher precision classifier in Macrel for AMP prediction from real data to minimize the number of false positives in the overall pipeline.

Antimicrobial peptides, as they are likely to interact with cell membranes, can cause hemolysis, which can impact its potential uses, particularly in clinical settings (Zhang & Gallo, 2016; Ruiz et al., 2014; Oddo & Hansen, 2017). Therefore, for convenience, Macrel includes a classifier for hemolytic activity. This model is comparable to the state-of-the-art (see Table 2).

High specificity is maintained when controlling for homology

Although we used out-of-bag estimates (see “Methods”) to control for exact overlap between training and testing sets in the previous section, we still included similar sequences in training and testing, leading to an overestimate of generalization potential. To control for this effect, Macrel and three methods (those where the ability to retrain the model was provided by the original authors) were tested using a stricter, homology-aware, scheme where training and testing datasets do not contain homologous sequences between them (80% or higher amino acid identity, see “Methods”).

As expected, the measured performance was lower in this setting, but Macrel still achieved perfect specificity. Furthermore, this specificity was robust to changes in the exact proportion of AMPs:non-AMPs used in the training set, past a threshold (see Table S4; Fig. 3D). Considering the overall performance of iAMP-2L model, future versions of Macrel could incorporate a combination of features from Macrel and iAMP-2L.

Using blastp as a classification method was no better than random, confirming that homology-based methods are not appropriate for this problem beyond very close homologs.

Macrel recovers high-quality AMP candidates from genomes and metagenomes

To evaluate Macrel on real data, we ran it on 484 reference genomes that had previously shown to be abundant in the human gut (Coelho et al., 2019). This resulted in 171,645 (redundant) smORFs. However, only 8,202 (after redundancy removal) of these were classified as potential AMPs. Spurio (Höps, Jeffryes & Bateman, 2018) classified 853 of these (circa 10%) as likely spurious predictions.

Homology searches confirmed 13 AMP candidates as homologs from those in DRAMP database. Among them, a Laterosporulin (a bacteriocin from Brevibacillus), a BHT-B protein from Streptococcus, a Gonococcal growth inhibitor II from Staphylococcus, and other homologs of antimicrobial proteins. Seven of these confirmed AMPs were also present in the dataset used during model training.

To test Macrel on short reads, we simulated metagenomes composed of these same 484 reference genomes, at three different sequencing depths (40, 60, and 80 million reads) using abundance profiles estimated from six different real samples (Coelho et al., 2019) (for a total of 18 simulated metagenomes). The number of predicted smORFs increased with sequencing depth, with about 20k smORFs being predicted in the case of 80 million simulated reads (see Fig. 4A). Despite this large number of smORF candidates, only a small portion of them (0.17–0.64%) were classified as putative AMPs.

In total, we recovered 1,376 sequences for a total of 547 non-redundant AMPs predicted from the simulated metagenomes. Of these, only 44.5% are present in the underlying reference genomes. However, after eliminating singletons (sequences predicted in a single metagenome), this fraction rose to 80.4%. Thus, we recommend singleton elimination as a procedure to reduce false positives. Although fewer than half of pre-filtered AMP predictions were present in the reference genomes, only 12% of all predictions were marked as spurious by Spurio (see “Methods”). We manually investigated the origin of these spurious predictions and found that most of the spurious peptides are gene fragments from longer genes due to fragmentary assemblies or even artifacts of the simulated sequencing/assemblies. Interestingly, even the mispredictions were confirmed as AMPs by using the web servers of the methods tested in benchmark. In fact, circa 90% of all AMP candidates (including spurious predictions) were co-predicted by at least one other method than Macrel, and 61% were co-predicted by at least other four methods.

Having established that the rate of false positives can be kept low after singleton elimination, we investigated the recall of Macrel, namely whether it was able to recover the AMPs that were present in the underlying genomes. Post hoc, we estimated that almost all (97%) were in genomes with a coverage of at least 4.25 (while only 9% of the non-recovered AMPs had this, or a higher, coverage, see Fig. 4E). Nonetheless, in some exceptional cases, even very high coverage was not sufficient to recover a sequence.

Macrel predicts putative AMPs in real human gut metagenomes

To evaluate Macrel on real data, we used 182 previously published human gut metagenomes (Heintz-Buschart et al., 2016). Of these, 177 (97%) contain putative AMPs, resulting in a total of 3,934 non-redundant sequences (see Table S3). The fraction of smORFs classified as AMPs per metagenome ranged 0.1–1.65%, a range similar to that observed in simulated metagenomes (see “Macrel Recovers High-Quality AMP Candidates from Genomes and Metagenomes”).

After eliminating singletons, 1,373 non-redundant AMP candidates remained, which we further tested with alternative methods. In total, 92.8% of the AMPs predicted with Macrel were also classified as such by at least one other classifier, and 65.5% of the times, half or more of the tested state-of-art methods agreed with Macrel results (see Table S5). iAMPpred and CAMPR3-RF showed the highest agreement and co-predict 74.4% and 65.7% of the AMPs predicted by Macrel, respectively.

Ten percent of all predicted AMPs (414 peptides, or 10.5%) were flagged as likely spurious (see “Methods”). The fraction of non-singleton AMPs predicted as spurious was slightly lower (8%, a non-significant difference). Our final dataset, after discarding both singletons and smORFs identified as spurious (see “Methods” and Table S3), consists of 1,263 non-redundant AMPs.

As 36 metatranscriptomes produced from the same biological samples are also available, we quantified the expression of the 1,263 AMP candidates. Over 53.8% of the predicted AMPs had detectable transcripts (see Fig. S2). For 72% of these, transcripts were detected in more than one metagenome.

Taken together, we concluded that Macrel could find a set of high-quality AMPs candidates, which extensively agrees with other state-of-art methods, many of which are being actively transcribed.

Macrel requires only moderate computational resources

The tests reported here were carried out on machines corresponding to standard consumer hardware (Amazon WebServices, t2.large, which feature 8 GB of RAM and 2 cores) to show that Macrel is a pipeline with modest computational requirements. The execution time, although naturally dependent on the input size, was not greater than 25.5 h (recall that the largest simulated metagenomes contained 80 million reads, see Fig. 4F). The assembly steps consumed 75–80% of the execution time, while read trimming and gene prediction occupied another considerable part (10–15%).