BacARscan: an in silico resource to discern diversity in antibiotic resistance genes

Abstract Antibiotic resistance has escalated as a significant problem of broad public health significance. Regular surveillance of antibiotic resistance genes (ARGs) in microbes and metagenomes from human, animal and environmental sources is vital to understanding ARGs’ epidemiology and foreseeing the emergence of new antibiotic resistance determinants. Whole-genome sequencing (WGS)-based identification of the microbial ARGs using antibiotic resistance databases and in silico prediction tools can significantly expedite the monitoring and characterization of ARGs in various niches. The major hindrance to the annotation of ARGs from WGS data is that most genome databases contain fragmented genes/genomes (due to incomplete assembly). Herein, we describe an insilicoBacterial Antibiotic Resistance scan (BacARscan) (http://proteininformatics.org/mkumar/bacarscan/) that can detect, predict and characterize ARGs in -omics datasets, including short sequencing, reads, and fragmented contigs. Benchmarking on an independent non-redundant dataset revealed that the performance of BacARscan was better than other existing methods, with nearly 92% Precision and 95% F-measure on a combined dataset of ARG and non-ARG proteins. One of the most notable improvements of BacARscan over other ARG annotation methods is its ability to work on genomes and short-reads sequence libraries with equal efficiency and without any requirement for assembly of short reads. Thus, BacARscan can help monitor the prevalence and diversity of ARGs in microbial populations and metagenomic samples from animal, human, and environmental settings. The authors intend to constantly update the current version of BacARscan as and when new ARGs are discovered. Executable versions, source codes, sequences used for development and usage instructions are available at (http://www.proteininformatics.org/mkumar/bacarscan/downloads.html) and GitHub repository (https://github.com/mkubiophysics/BacARscan).


Introduction
Antibiotics-the 'miracle drugs'-were modern medicine's revolutionary discovery. However, their frequent and indiscriminate use for clinical, veterinary, agricultural, poultry and other purposes has resulted in bacterial resistance for almost all the significant antibiotic classes [1,2]. Antibiotic resistance has become a global phenomenon spanning humans, animals and the environment. Hence, regular surveillance of antibiotic resistance genes (ARGs) in microbes and metagenomes from human, animal and environmental sources is vital to understand the epidemiology of ARGs and to foresee the emergence of new antibiotic resistance determinants.
A rapid decline in the cost of high-throughput whole-genome sequencing (WGS) technologies has yielded a plethora of information about ARGs in the microbial population [3][4][5]. Also, metagenomic studies have added immensely to our knowledge of the microbial ARG pool [6]. As the cost of genome sequencing is decreasing, it is expected that classical microbiological and molecular methods of antibiotic resistance determination might be replaced by WGS-based approaches to microbial ARGs identification methods [7][8][9]. Antibiotic resistance databases and in silico ARG prediction tools leverage the information generated through classical microbiological and molecular experimental methods to determine antibiotic resistance. This means the efficiency of an in silico tool of microbial ARGs identification methods depends on the quality of data generated by classical microbiological and molecular methods of antibiotic resistance determination. Currently, a significant hindrance in identifying and annotating ARGs from WGS data is the abundance of fragmented genes/ genomes in the genome databases (due to incomplete assembly). In the metagenomic studies of ARGs, the source of genetic material is either the environment, or site of infection or microbiome. The DNA is extracted from the sample and sequenced using the whole-genome shotgun sequencing technique. The fragments are then assembled into contigs. Sometimes a very high genomic diversity may result in a low level of genome coverage that may lead to a small number of contigs and many unassembled sequencing reads (singletons). Thus, an ideal bioinformatics tool, which can detect ARGs in WGS data, should also work effectively on short sequencing reads and fragmented contigs.
ARDB was the first in silico resource that provided a centralized repository to characterize ARGs. ARDB is no longer updated, and its data are incorporated in the CARD. Resfams is a database of hidden Markov models (HMMs) developed using the protein families associated with antibiotic resistance [16]. Meta-MARC is based on hierarchical HMMs, which can predict AMR in metagenomic data (either a short read or a longer assembled contig) into resistance class, group and mechanism. AMRFinderPlus identifies acquired AMR genes and resistance-associated point mutations in protein or assembled nucleotide sequences. CARD identifies and annotates ARGs using BLAST or Resistance Gene Identifier (RGI) modules. A similarity search is performed using BLAST against the CARD reference sequences in the CARD sequence, while the RGI option predicts ARG(s) based on homology and single nucleotide polymorphism (SNP) models [14,19]. Several studies have pointed out that BLAST works well in comparing sequences with a high degree of similarity (60% or higher) but does not identify a distant homolog [20,21]. This limits the capacity of BLAST-based methods to annotate novel ARGs. Because sequence alignment methods like BLAST work well in comparing sequences with a high degree of similarity (60% or higher) but do not identify a distant homologue [20,21]. Pawlowski et al. [22] reported five ARGs in a cave bacterial isolate, Paenibacillus sp. LC231, which does not have characterized homologues, and CARD could not identify them [22].
One major limitation of ARG databases is their bias towards human pathogens and model organisms, which is understandable because most research efforts are directed towards them [23]. In non-model or non-culturable organisms, remote homologues or highly divergent forms of ARGs of human pathogens and model organisms may be present. This effectively means either the absence or very less number of ARGs of lesser-studied microbes/pathogens. Hence identification of remote homologues or novel ARGs is difficult using the traditional sequence similarity methods [24]. This bias also makes it difficult to identify ARGs in less commonly studied microbes [25].
Also, several resources can identify/characterize resistance only against the prominent/frequently prescribed antibiotics [20]. Thus, we have described a new in silico tool for rapid monitoring, characterization and surveillance of all bacterial ARGs in the present study. Named as Bacterial Antibiotic Resistance scan (BacARscan), this tool has the edge over its predecessors as it can also discern ARGs in short sequencing reads and fragmented contigs. BacARscan can be easily integrated into a user-defined ARG annotation pipeline for the detection of ARG variants in microbial genomes. We also compared BacARscan with other ARG annotation tools: AMRFinderPlus, Meta-MARC, RGI-CARD and Resfams. The results indicated less false-positive (FP) prediction of ARG by BacARscan vis-a-vis other methods.

Overall study design
The complete pipeline of the tool is shown in Fig. 1 with a detailed description.

Dataset collection and cleaning
In the present study, the data on ARGs/proteins were retrieved from various antibiotic resistance databases such as ARG-ANNOT [26], CARD [14], CBMAR [12], INTEGRALL [27], RAC [28], Tetracycline þ MLS nomenclature, UCARE [29], Lahey Clinic, Resfams [16], ResFinder [13], HMP [30], LacED [31], MvirDB [32], Institut Pasteur, Patric [33] and published research papers. All ARGs were first grouped based on the mechanisms they employ to inactivate the effect of antibiotics profile. For example, hydrolysis is used to inactivate the b-lactams hence all b-lactamases were grouped; similarly, all antibiotic efflux proteins were grouped. Members of each group were then clustered using BLASTClust with an identity threshold of 90% and coverage of 95%. After removing protein clusters containing less than five sequences or fragmented sequences, we are finally left with 254 protein sequence clusters. Clusters with less than five sequences were removed because a small number of sequences might not provide all mutational information to a profile HMM.
Building protein profile ARGhmms from sequences of antibiotic resistance proteins All the protein sequences within the 254 clusters were subjected to multiple sequence alignment (MSA) using Muscle 3.8 [34] at default parameters. To check the accuracy of an alignment, we manually checked each alignment and removed sequences that were too divergent [35,36]. Each MSA was then converted into a profile HMMs using the hmmbuild module of HMMER (http:// hmmer.wustl.edu/) (version 3.1) at default parameters [37]. The library of protein profile ARG HMMs was referred to as pARGhmm.
Construction and functional annotation of the nucleotide version of ARGhmm (nARGhmm) In order to develop a tool that can be used directly on genomics/ metagenomics sequence data without translating into protein sequences, we created a nucleotide version of the profile HMMs library, named nARGhmm. nARGhmm was created using nucleotide sequences back translated from protein to gene sequences using the Backtranseq (Predicts the potential nucleic acid sequence from which a given peptide sequence originated) program EMBOSS. The procedure followed to build nARGhmm was the same as that followed for pARGhmm. For each antibiotic resistance cluster, we constructed two profile HMMs, one for proteins (pARGhmm) and another for genes (nARGhmm). Collectively, pARGhmm and nARGhmm were named as BacARscan. The schema used for building BacARscan from the set of curated antibiotic resistance proteins (ARPs) and gene sequences is shown in Fig. 1.

Dataset I: evaluation dataset
The efficiency of BacARscan was evaluated using a positive and negative dataset. The negative dataset was composed of protein clusters formed during clustering by BLASTclust and had less than five sequences. Since these proteins were not used to create the profile HMMs, they were designated as negative datasets. In order to confirm that the negative dataset is composed of non-ARG sequences, we did gene ontology (GO) term enrichment analysis of negative dataset. The results of enrichment analysis indicated that the most frequently occurring GO terms were not related to AMR, barring a minor fraction showing 'beta-lactamase activity' as molecular function and 'beta-lactam antibiotic catalytic process' as biological processes ( Supplementary Fig. S1). We Figure 1: The complete workflow depicting the methodology used for the development of profile HMMs-based prediction tool BacARscan.
In silico tool for AMR detection | 3 further confirmed this using RGI-CARD. Out of 65 816 proteins of negative dataset, 6007 ($10%) were characterized as ARG through RGI-CARD out of which 760 were categorized as perfect, 450 as strict and 4797 as loose. It was observed that a majority of ARGs (1574) were characterized as b-lactamases (perfect: 186, strict: 94 and loose: 1294). Details of the classification are shown in Supplementary Table S1. We feel this might have happened because b-lactamase is a highly diverse group of ARG, and the negative dataset has those b-lactamases that are less abundant in nature. Results of RGI-CARD and GO-term enrichment analysis indicated that the negative dataset is mostly (90%) composed of non-ARG sequences. Sequences used to build the pARGhmm, that is protein clusters that contained five or more sequences, were called positive datasets. The combined dataset of fulllength positive dataset sequences and both full-length and partial-negative dataset sequences were used for performance evaluation.

Dataset II: independent dataset
For benchmarking, we used two independent datasets consisting of (i) penicillin-binding proteins (PBPs) or DD-peptidase proteins, retrieved from the UniProtKB database on 15 April 22 using keyword search (penicillin-binding proteins). This dataset contains 60 reviewed non-fragmented PBPs whose existence was established at the protein level. (ii) Non-antibiotic resistant bacterial efflux (non-ARE) proteins were collected from our earlier work [38]. The final benchmark independent datasets comprised 60 PBPs and 389 non-ARE protein sequences.

Dataset III: annotation of ARGs in different strains of ESKAPE pathogens
To assess the ability of BacARscan to discern the diversity of ARGs among different ESKAPE pathogens (Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter baumannii, Pseudomonas aeruginosa and Enterobacter species), we annotated the ARGs present in five strains of an ESKAPE pathogen using BacARscan and compared its performance with other in silico resources, namely RGI-CARD and Resfams. For annotation, we used proteomes of those ESKAPE pathogens whose genomes were completely sequenced, with complete chromosomal and plasmid protein sequences available at National Center for Biotechnology Information (NCBI). Detailed information about ESKAPE pathogens used in this study is mentioned in Supplementary Table S2. During evaluation, the best-hit approach (lowest e-value and highest HMM score) was used for HMM-based methods (BacARscan and Resfams). But in the case of RGI-CARD, only perfect (e-value ¼ 0, sequence identity $99-100% and bit score above 500) and strict hits (e-value¼ $1e-80-1e-90, sequence identity $90-65% and bit score at least 500) were considered.
Dataset IV: annotation of ARGs using short sequence reads The performance of BacARscan was also evaluated on a dataset of short reads generated from the back-translation of protein sequences of positive (used for building profile HMMs: !5 cluster sequences) and negative (not used for building profile HMMs: <5 cluster sequences) evaluation dataset. The length of each read was 100 nucleotides (nt), and from each gene sequence, 20 short sequences were generated for evaluation.

Dataset V: validation dataset
Simulated short-read data of ARG.
This dataset was used to benchmark BacARscan vis-a-vis other ARG prediction and annotation methods. It was generated using the gene sequence of the 'protein homolog' model of the CARD (date: 29 July 2022). The 'protein homolog' model contains 4422 ARG sequences that do not include mutation as a determinant of resistance. Using InSilicoSeq (https://insilicoseq.readthedocs.io/) [39], short reads of 151-nt length at 20Â coverage were simulated, imitating the NovaSeq platform of Illumina. A total of 100 000 reads were randomly picked for benchmarking.
Simulated non-ARG short-read data.
We also created a non-ARG short-read dataset using the complete genome of a probiotic strain of Enterococcus faecium Strain T-110 (NCBI Genome Accession Number: CP006030). As per Natrajan and Parani, this strain, called T-110, does not contain ARGs against all major clinically relevant antibiotics [40]. We used the genome of T-110 to simulate a database of 2 million short reads of 151 bases using InSilicoSeq (https://insilicoseq.readthedocs.io/), imitating the NovaSeq platform of Illumina. The comparative evaluation was carried out among BacARscan, Meta-MARC and ResFinder using randomly selected 2M reads.

Performance measure parameters
During an evaluation, only predictions with e-values 1e-6 were considered significant because maximum precision, recall and F-score were observed at this value (details in comparative evaluation of prediction efficiency of pARGhmm and nARGhmm). To evaluate the efficiencies of both nARGhmm and pARGhmm, we classified the prediction results into four categories: true positives and true negatives (TP and TN, respectively), FP and false negatives FN, respectively. TP and TN represented the correct predictions, while FP and FN represented incorrect predictions. When an ARG sequence was rightfully associated/predicted to its actual class, the prediction was classified as TP, and when a non-ARG was predicted as non-ARG, it was classified as TN. When an ARG and a non-ARG were predicted as non-ARG and ARG, the prediction results were labelled FN and FP, respectively. Using this classification criterion, we calculated the precision, recall, F-score, falsepositive rate (FPR) and true-negative rate (TNR) or specificity to measure the prediction capability of BacARscan. i) Precision: It is the ratio of correctly predicted positive observations to the total number of predicted positive observations. High precision indicates a low FP prediction. Precision was calculated as: ii) Recall or sensitivity: It is the ratio of correctly predicted positive observations to the total number of positive samples submitted for prediction. High recall indicates a highly efficient predictor with very low false predictions. It was calculated as: iii) F-measure or F 1 or F-score: In order to balance the precision and recall values due to the unequal composition of positive and negative datasets, F-score was also calculated [41]. It is the harmonic mean of precision and recall, and F-score was calculated as: iv) FPR: The FPR is calculated as the ratio between the number of negative events wrongly categorized as positive (FPs) and the total number of actual negative events (regardless of classification).
v) TNR or Specificity: The TNR is the proportion of genuinely negative samples that give a negative result using the test in question. It is also known as specificity. It was calculated as:

Results and discussion
Performance assessment of BacARscan (pARGhmm and nARGhmm)

Evaluation of protein dataset
To assess the performance of pARGhmm, we used an evaluation dataset (positive and negative). For this, two different approaches were used to evaluate the performance, that is the top hit and multiple top hits (such as top 3rd, top 5th, top 7th hits, etc.) at different e-value and prediction score cut-offs. When the top hit (based on an e-value and prediction score) was used to evaluate the prediction capability, 228 TP and 26 FP predictions were found (Table 1). Furthermore, the performance of pARGhmm was also evaluated using multiple search results based on the lowest evalue and highest prediction score, and the consensus approach decided the outcome. With the majority approach, the maximum performance was achieved when the top five hits were used to evaluate the performance. The best precision and F-measure were 92.12% and 95.90%, respectively. Overall, the assessment results showed that the prediction performance increased till the decision was made based on the top five hits (based on e-value and prediction score); afterwards, the performance started decreasing. The consistency in search efficiency of each pARGhmm was evaluated using the leave-one-out cross-validation (LOOCV) approach. During LOOCV, one sequence was excluded from the process of model building. Then the performance of the model was evaluated against a dataset containing (i) excluded sequences, (ii) sequences used for building 253 HMMs and (iii) ARP sequences not used for building of pARGhmm. Using this approach, a perfect consistency ($100%) was observed in the performance of the majority of profile HMMs (Supplementary Fig. S2).

Evaluation of short sequencing reads
The performance of nARGhmm was also evaluated on a dataset consisting of short sequences of 100-nt length. The short sequence library was generated from used and non-used sequences for building the profile HMMs, which constitute the positive and negative evaluation datasets. Using the first hit for evaluation, the precision and F-measure values were 90.94 and 95.25, respectively (Table 1). When the performance of nARGhmm was evaluated using the majority approach, the best performance was achieved when the top five hits were used to quantify the performance.
Thus the comparative evaluation of prediction efficiency of nucleotide and protein modules of BacARscan (nARGhmm and pARGhmm) indicated a comparable performance in terms of both precision and F-measure ( Table 1). The results also indicated the comparable performance of BacARscan on fragmented and full-length datasets, indicating that BacARscan is an efficient tool for finding ARGs in both protein and gene datasets.

Assessment of BacARscan ARG prediction capability vis-a-vis other methods
We evaluated BacARscan with other existing methods using three datasets which include (i) full-length gene/protein homologous non-ARG sequences, (ii) nucleotide short reads of ARGs and (iii) nucleotide short reads of non-ARGs.

Benchmarking of BacARscan on homologous non-ARPs/ ARGs
To verify the capability of BacARscan to correctly differentiate ARPs/ARGs from homologous non-ARGs/proteins, we constructed an independent dataset consisting of PBPs and non-antibiotic efflux proteins [38] (non-ARE). Both PBPs and b-lactamases belong to the superfamily of serine penicillinrecognizing enzymes and have similar conserved protein folds [42,43]. It is pertinent to mention that PBPs and non-ARE proteins were not used to construct the HMMs. As done by the previous researchers [13][14][15][16], during benchmarking, we also The overall assessment of BacARscan using both versions of profile HMMs: Scanning of ARPs and ARGs using both versions of profile HMMs (pARGhmm and nARGhmm) against evaluation dataset (positive and negative protein sequences dataset) and short-read sequences independent dataset. The first column shows a couple of alternative top hits like top 1st, 3rd, 5th, 7th hits, etc. TP and FP are the calculation of accuracy of our prediction method. Precision known as exactness was calculated with the above predicted values of TP and FP rate, respectively. Recall and specificity were calculated and we found 100 and 0 in all the cases. F-measure known as harmonic mean was calculated with the help of all considered performance measures.
adopted the best-hit approach with an e-value cut-off of 1e-20. During analysis, out of a total of 60 PBP homologous gene/protein sequences of independent dataset, BacARscan, AMRFinderPlus, Meta-MARC, RGI-CARD, Resfams predicted 6, 12, 9, 15 and 4 as b-lactamases ARGs. In contrast, in non-ARE gene/protein sequences of independent dataset, BacARscan, AMRFinderPlus, Meta-MARC, RGI-CARD, Resfams predicted 23, 37, 26, 91 and 24, out of the 389 non-ARE genes/proteins as efflux-related ARGs ( Table 2). In terms of FPR and TNR, though Resfams outperformed all other methods when evaluated using PBP, the performance of BacARscan was also comparable.
On the other hand, on the non-ARE dataset, the performance of BacARscan was maximum, with 5.9% FPR and 94.08% TNR ( Table 2). Overall, results revealed that BacARscan could fairly discriminate between the ARGs and non-ARG consisting of homologous proteins.

BacARscan has a lower FPR on non-ARGs short-reads dataset than other algorithms
At the default e-value threshold of 1e-6, BacARscan, Meta-MARC and ResFinder predicted 3979, 22 331 and 1912 short reads, respectively, to belong to ARGs (Table 4). These reads belong to 19, 56 and 5 ARG genes of T-110, respectively, and 3, 18 and 3 ARGs, respectively, in T-110. When a more stringent e-value cut-off of 1e-20 was used, as done in earlier work [13,14,16], the number of reads predicted by BacARscan, Meta-MARC and ResFinder was reduced to 238, 9034 and 1648, respectively. When a highly stringent e-value of 1e-50 was used, the no reads were predicted by BacARscan and Meta-MARC as reads of ARG genes, but ResFinder still predicted 1500 reads as ARG. At all e-values, the FPR was <1% except ResFinder at an e-value of 1e-6. This means the rate of TN prediction was nearly 99%.

Annotation of ARGs in ESKAPE pathogens
The ESKAPE pathogens are globally recognized as a leading cause of nosocomial multidrug-resistant bacterial infections. The results of ARG annotations by BacARscan revealed a predominance of antibiotic efflux-related ARGs, followed by ARGs encoding b-lactamases or acetyltransferases in five strains of the ESKAPE pathogens (Fig. 2). Similar results were also reported in an earlier published study [44].
Our results also revealed that Resfams and BacARscan predicted a similar number of ARGs (with slight variations) for each antibiotic resistance mechanism (Fig. 2). We also observed that the number of ARGs predicted by HMM-based methods was much higher than the RGI-CARD (Table 5). Perfect and strict hits of RGI-CARD are categorized on the basis of their similarity to the matching CARD references. A perfect hit of RGI-CARD has all amino acids matching CARD reference, and a strict hit has a bitscore of !500. This means RGI-CARD annotated both perfect and strict hits using a very stringent e-value cut-off. We feel the difference in the number of ARGs predicted by HMM-based methods and RGI-CARD is due to the difference in e-value cut-off. When we compared the number of ARGs predicted by HMM-based methods at very stringent e-value cut-off (using 1e-80 or 1e-100), we found the number of ARGs was similar to the RGI-CARD (Supplementary Table S3). This indicated that at a very low evalue threshold, the performance of BacARscan and Resfams is comparable to the RGI-CARD. The results also indicated that a very low e-value should be used for proteome-wide annotation of ARGs.

Description of the BacARscan web server
The BacARscan tool is available as a web server at (http://protei ninformatics.org/mkumar/bacarscan), where the users can upload the gene/protein sequences for prediction and annotation. The output of BacARscan contains the complete annotation of the query sequence. BacARscan can process a maximum of 10 sequences in one go. For whole-genome/metagenomic/proteome scale annotation, a standalone version is required or available in the web server's download section or GitHub repository (https:// github.com/mkubiophysics/BacARscan).

Potential use of BacARscan
Most traditional antibiotic resistance determination methods are based on laboratory tests but the time and resources required to generate this data sometimes make it unsuitable for quick monitoring. The development of high-throughput sequencing methods has opened an alternative quick and cost-effective method of antibiotic resistance determination. BacARscan is an effort to provide an additional way to identify the ARGs in an -omics (proteomics/genomics and metagenomic) datasets. BacARscan can also be combined with traditional surveillance and thus can complement the traditional methods of ARG annotation. The current version of BacARscan supports prediction using only 254 ARG families, but in the future, we will extend it to the new ARGs or families. We hope that BacARscan will help predict ARGs and help in the progress of studies related to antimicrobial resistance (AMR).

Concluding remarks
In the present study, we have described a novel ARG annotation in silico resource named BacARscan, which can be used for rapid monitoring, surveillance and characterization of antibiotic resistance determinants in both genomics and proteomic datasets. Comparison with other in silico resources like AMRFinderPlus, Meta-MARC, Resfams and CARD revealed that BacARscan's ability to discern ARGs in -omics datasets was much more significant than its predecessors. One of the most notable improvements of BacARscan over other ARG annotation methods is its ability to work on genomes and short-reads sequence libraries with equal efficiency and without any requirement for assembly of short reads. Thus, BacARscan would help monitor the prevalence and diversity of ARGs in microbial populations and metagenomic samples from animal, human and environmental settings. Though the performance of BacARscan on different datasets looks promising, we understand that without a constant update, it would not remain helpful to the scientific community. The authors intend to constantly update the current version of BacARscan as and when new ARGs are discovered.  The number of ARGs predicted by BacARscan, Resfams and RGI-CARD. The e-value threshold: 1e-6 was used for BacARscan and Resfams, while in the case of RGI-CARD, perfect and strict cut-offs were used for the selection of the hits, respectively.

Supplementary data
Supplementary data is available at Biology Methods and Protocols online.

Data and software availability
The tools and all datasets used for evaluation are freely accessible without any restriction on the download page of the web server: http://proteininformatics.org/mkumar/bacarscan/downloads.html