Defining objective clusters for rabies virus sequences using affinity propagation clustering

Rabies is caused by lyssaviruses, and is one of the oldest known zoonoses. In recent years, more than 21,000 nucleotide sequences of rabies viruses (RABV), from the prototype species rabies lyssavirus, have been deposited in public databases. Subsequent phylogenetic analyses in combination with metadata suggest geographic distributions of RABV. However, these analyses somewhat experience technical difficulties in defining verifiable criteria for cluster allocations in phylogenetic trees inviting for a more rational approach. Therefore, we applied a relatively new mathematical clustering algorythm named ‘affinity propagation clustering’ (AP) to propose a standardized sub-species classification utilizing full-genome RABV sequences. Because AP has the advantage that it is computationally fast and works for any meaningful measure of similarity between data samples, it has previously been applied successfully in bioinformatics, for analysis of microarray and gene expression data, however, cluster analysis of sequences is still in its infancy. Existing (516) and original (46) full genome RABV sequences were used to demonstrate the application of AP for RABV clustering. On a global scale, AP proposed four clusters, i.e. New World cluster, Arctic/Arctic-like, Cosmopolitan, and Asian as previously assigned by phylogenetic studies. By combining AP with established phylogenetic analyses, it is possible to resolve phylogenetic relationships between verifiably determined clusters and sequences. This workflow will be useful in confirming cluster distributions in a uniform transparent manner, not only for RABV, but also for other comparative sequence analyses.


Introduction
Virus taxonomy differs from other types of biological classification because the International Committee on Taxonomy of Viruses (ICTV) not only regulates a Code of Nomenclature, but also considers and approves the creation of novel virus taxa (currently orders, families, subfamilies, genera and species). Thanks to long-lasting efforts of the ICTV [1], the classification of viruses has become clearer and more transparent [2]. Lyssaviruses, negative sense RNA viruses, represent one of 18 currently known virus genera within the family Rhabdoviridae of the order Mononegavirales. Based on the species concept in virus taxonomy whereby demarcation criteria are established to discriminate between different virus species within a genus [3], the Rhabdovirus study group defines a new species of lyssavirus among other things by more than 18-20% nucleotide divergence within the N-gene compared to the existing lyssavirus species [4]. The Lyssavirus genus comprises 14 recognized and two putative lyssavirus species [2,5,6], of which rabies lyssavirus represents the prototypical lyssavirus species.
Numerous viral variants of rabies virus (RABV) cause tens of thousands of human deaths annually on a global scale [7]. Nevertheless, there are no further diversification criteria below the species level [2], and for lyssaviruses not even a standard definition for genetic grouping (e.g. lineage, clade, variant, strain, cluster) has been established. In this study, the term 'cluster' is utilized throughout the manuscript to define sub-species demarcation.
In order to genetically characterize and sub-classify RABV isolates, multiple studies were conducted, resulting in approximately 21,000 datasets of partial and full genome RABV sequences obtained from the databases of the NCBI. Different studies analyzed the extent of relationship among selected samples using phylogenetic analyses to verify the results. Historically, the majority of phylogenetic analyses have been conducted at regional levels, i.e. for Europe [8,9], African regions [10,11], Asia [12] and The Americas [13,14]. Phylogenetic analyses used mostly N and G gene sequences (Table 1). This is due to a number of factors including the submission of diagnostic RT-PCR amplicon sequences [15]. In fact, the phylogenetic analysis of partial genome sequences requires less computational power and is more cost-effective. Subsequently, even broad-scale phylogenetic analyses are based on partial genes (Table 2) [16]. This approach was further supported by a study indicating that all lyssavirus genes are equivalent for phylogenetic analysis [17]. However, in reality the sequences submitted to NCBI have minimal or no sequence overlaps, resulting in datasets which cannot be compared; furthermore, phylogenetic trees often had low statistical support, as demonstrated previously for Arctic RABVs [18]. Recent analyses have suggested that full genomes provide better statistical support and are a comprehensive instrument for addressing the evolution, spread and genome-wide heterogeneity of RABV [19]. Increasing sequencing capacities including high-throughput sequencing enable generation of full RABV genome sequences for phylogenetic and evolutionary studies ( Table 2). Complete genome analyses provide the ultimate opportunity to detect specific nucleotide substitution patterns and identification of specific motifs across all genes [20]. In addition, full genomes are useful for comparison of selection pressures on different genes and are therefore helpful to analyze cross-species transmission events [20][21][22][23][24][25], and endemic transmission amongst others [26]. The most comprehensive analysis in this respect investigated the evolutionary history of RABV by contributing 170 newly generated full genome sequences [21].
In all previous phylogenetic studies, cluster allocation was either based on host species, region of origin or statistical (bootstrap) support. However, the allocation of RABVs into clusters can be highly subjective, because the thresholds of statistical support vary depending on the respective dataset, often the sequence length. In combination with the individual designation of clusters at a regional and local level, comparison and combination of published phylogenetic studies on RABV is often difficult. In addition, the increasing number of available RABV sequences represents a challenge for conventional computation of phylogenetic inferences. Similar problems have been described for other viruses and alternative solutions developed [27][28][29]. However, to reveal unbiased criteria for cluster definitions we preferred a workflow based on a novel non-hierarchical mathematical clustering method: affinity propagation (AP) clustering in combination with standardized phylogenetic analyses. AP is a tool that was developed for clustering similarity measures between all pairs of input samples based on the concept of "message passing" between data points [30]. The method does not require a vector space structure and so called 'exemplars', samples that are most representative for a cluster, are chosen among the observed data samples and not computed as hypothetical averages of cluster samples. These characteristics make AP clustering particularly suitable for applications in bioinformatics [31]. Therefore, this approach has already been successfully applied for various tasks in bioinformatics, e.g. for microarray and gene expression data [30,32,33]. Here, an extended panel of existing and newly obtained full genome RABV sequences was used to demonstrate the application of AP for RABV clustering and the results compared with previous studies.

RABV full genome sequences
For AP clustering, a total of 516 RABV full-genome sequences were obtained from the NCBI sequence database using the tools implemented in Geneious ( [34], version 10.0.9, http://www. geneious.com). All datasets were manually checked for completeness and the respective missing metadata, e.g. year of isolation, geographical origin, was manually completed where possible from literature (S1 Table). In addition, a further panel of 46 RABV isolates from previously underrepresented geographical areas such as Near East, Europe, Southern America and some African regions was sequenced using high-throughput sequencing (Table 1). RABVs were obtained from the virus archive of the WHO Collaborating Centre for Rabies Surveillance and Research at the Friedrich-Loeffler-Institut (FLI), Greifswald, Germany, or from the WHO Collaborating Centre for the Characterization of Lyssaviruses, Animal and Plant Health Agency (APHA) Weybridge, United Kingdom. All datasets were subsequently stored in a Geneious (version 10.0.9) database for further use. Due to the high passaging and modifications that are required to produce vaccine strains, these data are not included.

NGS sequencing
To obtain full genome sequences of 46 isolates, high-throughput sequencing at the two reference laboratories was applied as follows: FLI NGS methods were conducted as described previously [35]. After total RNA extraction from cell culture supernatant DNase (Qiagen, Hilden, Germany) treatment was performed as recommended by the supplier. Briefly, total RNA was extracted from cell culture supernatant using RNeasy Mini Kit (Qiagen) along with on-column DNase (Qiagen) treatment following the supplier recommendations. Subsequently, Agencourt RNAclean XP beads (Beckman Coulter, Fullerton, USA) were used to concentrate and clean the RNA. A maximum of 750 ng RNA was used as input for cDNA synthesis using the cDNA synthesis system kit (Roche, Mannheim, Germany) along with random hexamer primers (Roche). Sequencing libraries were generated using the SPRI-TE instrument with SPRIworks II cartridges (Beckman Coulter) and appropriate adapters as described elsewhere [36]. An Illumina MiSeq platform using a MiSeq reagent kit, version 3 (Illumina, San Diego, USA) was used for sequencing according to the manufacturer's instructions in 2x300 bp paired end mode. The raw reads were quality trimmed and mapped along an appropriate RABV reference using the 454 Sequencing Systems Software suite (version 3.0, Roche). Mapped reads were selected for de-novo assembly (454 Sequencing Systems Software suite) in order to generate full genome RABV sequences.
APHA NGS methods were conducted as described previously [37,38]. Briefly, RNA was extracted directly from clinical brain material using TRizol, then host genomic DNA and ribosomal RNA was depleted using DNAse (Qiagen) and Terminator 5'-Phosphate-dependent Exonuclease (Epicentre Biotechnologies, Madison, USA) respectively. Preparation of ds cDNA and sequencing libraries is described elsewhere [38]. Sequencing was carried out as above for FLI using 2x150bp paired end mode. An iterative mapping process was implemented as described previously [38] to generate RABV full-genomes.

Comparative sequence analysis
The sequences were annotated in Geneious ( [34], (version 10.0.9, http://www.geneious.com) and submitted to the European Nucleotide Archive (ENA) under study accession number PRJEB22369 (see Table 1 for isolate specific accession numbers). Altogether, 562 full genome RABV sequences were aligned codon-based using ClustalW [39] implemented in MEGA 6 [40]. Subsequently, MEGA 6 was used to conduct phylogenetic analyses, including model test and neighbor joining (NJ) phylogenetic tree calculations [41]. With regard to the model test, the Tamura 3-Parameter (T92) evolution model was applied to all datasets [42]. Both, gamma distribution (five categories) and invariant sites were considered. Statistical tree topology support was derived from 500 bootstrap replications.

AP clustering
Pairwise distances as calculated in the phylogenetic analysis were merged into a distance matrix and imported to the statistical software R [43]. For further analyses the package "apcluster" was used essentially as described [31]. By default, the AP clustering algorithm determines one sequence among the set of input elements for each potential cluster, which is most representative for this cluster. In AP terminology, these sequences are termed "cluster exemplars" [31]. Since this method was initially developed for analysis of similarity matrices, the distance matrix from the sequence alignment had to be converted by inverting the values. In addition, all values embedded in the matrix were squared to improve robustness and discriminatory power of the analysis. Subsequently, the AP algorithm computes the minimum (pmin) and maximum (pmax) of the input preference (p), which is defined as the tendency of each individual sequence to become a "cluster exemplar". To define the optimal input preference (p), the number of clusters for the complete preference range (pmin-pmax) was calculated stepwise [31]. Optimal input parameter for intraspecific analyses, i.e. the optimal number of clusters, was defined as the largest range of input parameters for which a constant number of clusters is calculated. This range is termed "plateau" throughout the manuscript. Methodologically, the beginning of the lower bound of the "two" cluster plateau cannot be defined and therefore the length of this plateau was not considered further. According to the now defined optimum input parameter, AP calculated the respective number of clusters and allocated any input sequence to only one of these.

Mapping of AP clusters
Results of AP were summarized on country level and exported as CSV file into ArcGis Desktop (Esri, version 10.3.1, Redlands, California, USA). Data was visualized in pie charts per country.

Combination of AP with phylogeny
Prior to phylogenetic analysis, each sequence dataset was completed with the assigned AP cluster using the group allocation function in MEGA 6. Within the phylogenetic tree, a color code for each assigned AP cluster was applied. For better visualization, the tree was condensed to the AP cluster level where applicable, i.e. each condensed branch in the dendrogram exclusively contained sequences unambiguously allocated to one AP cluster.

Phylogenetic analyses of RABV full genomes
All available RABV full genome sequences were used to assess the global evolutionary relationships of RABV. In total, 562 full genome sequences were analyzed, including 46 sequences newly obtained as part of this study. The RABV cohort originated from 76 different countries reflecting a broad host range mainly comprising bats, dogs, raccoons, and foxes (S1 Table). The individual sequence lengths ranged from 11,796 to 11,931 nucleotides in length, representing a lack of terminal sequences which can be problematic to obtain, particularly if amplicon sequence is utilized. Altogether, RABV sequences had a pairwise identity of 86.3% (range: 78.65% -99.97%) with an average G+C content of 45.3%. The low G+C content (under 50%), indicates an uneven distribution of nucleotides therefore not all evolutionary models are suitable for analysis. Therefore, the T92 model [42] was selected, also indicated by the model test in MEGA 6 [40]. Interestingly, the results described below were equally reproducible with other, more complex and time consuming, evolutionary, models e.g. Generalized Time Reversable (GTR).
Phylogenetic analyses based on neighbor joining (S1 Fig) confirmed previously described cluster distributions (Table 2), whereby two principal RABV clusters were observed. One cluster exclusively contains RABV sequences from The Americas, whereas the other cluster has a worldwide distribution and comprises RABV sequences of terrestrial mammals. In previous studies, the latter was further split into various sub-clusters based on either host or regional allocation, e.g. three [22], five [23], six [44], or up to eleven [45], illustrating methodological differences in the approach of sub-clustering. So far, cluster allocations have only partly been based on metadata including host species and geographic origin ( Table 2). Both metadata represent challenges, as host species do not necessarily identify the reservoir host but often the individual spill-over or 'dead end' host, whereas the spatial origin of samples is often not detailed enough to enable precise analysis. Additionally, cluster allocation can be by genetic distance within clusters [27] or by statistical (bootstrap) support of resulting nodes in phylogenetic trees [46]. However, delineation of clusters is not defined. In these analyses, 82% (463 out of 564) of all nodes showed bootstrap values higher than 70% (S1 Fig), widely recognized as a threshold for reliability [28], demonstrating that bootstrap support of nodes does not result in the delineation of meaningful clusters, if large numbers of full genome sequences are analyzed [47]. Using genetic distances for cluster allocation alone also does not resolve this issue, as no standards for lyssavirus species are defined for this purpose, and an arbitrarily defined threshold would also be subjective.

Affinity propagation clustering (AP)
To address the issue of undefined demarcation within the lyssavirus species, we used AP as a mathematical method for clustering RABV full genome sequences. Besides its computational efficiency [32], which substantially reduces the turnaround time, the main advantage is that it overcomes the described subjective criteria for cluster allocations with the help of mathematical algorithms. The results are non-hierarchically ordered clusters, which can have unequal cluster sizes [30,32]. By application of AP to pairwise genetic distances from an alignment of all 562 full RABV genomes, the most stable distribution after iteration over all possible input parameters was determined as four clusters (Fig 1A).
These were termed according to previous studies as Arctic/Arctic-like, Cosmopolitan, Asian (together forming the terrestrial cluster) and New World ( Table 2). Only one other phylogenetic analysis also identified four main clusters on the basis of full genome sequences [22]. However, other studies identified more than four main clusters for RABV on a global scale [20,21,24,25]. The four clusters identified in our study show a defined spatial distribution (Fig 2), confirming previous analyses that demonstrated the influence of geographic origin rather than the host [48].
In fact, allocation of clusters to known rabies reservoir hosts at this cluster level is problematic. For example, our data indicated that in the Arctic/Arctic like cluster the reservoir host arctic fox is only represented by 59.2% of RABV sequences. In the Cosmopolitan and Asian clusters 44% and 62% of the sequences originate from dog RABVs, respectively, while sequences in the New World cluster were derived to 51.7% from bat RABVs. The remaining sequences include further reservoir species as well as spillover hosts. But sub-dividing these main clusters, may lead to more host-specific clusters.

Combination of AP and classical phylogeny
As AP assigns sequences to robust and reproducable clusters based on transparent input parameter pre-analysis, we tested the hypothesis that a combination of AP with established classical phylogenetic analyses can overcome the inherent limitations of the latter methods alone. Therefore, a comparison of AP and classical phylogeny was undertaken to assess the overall extent of agreement between both methods. The obtained AP cluster distribution was transferred to the calculated phylogenetic dendrogram on an individual sequence basis, and visualized in a compressed tree (Fig 3).
As a result, the phylogenetic tree divides into two main branches. One branch comprises all full-genome sequences allocated to the New World AP cluster ( Table 2). The other comprises of sequences from Old World rabies cases, the further branching of sequences in the condensed phylogenetic tree is not congruent with the result obtained by AP except for the Arctic/ Arctic-like AP cluster. Sequences allocated to the Asian and Cosmopolitan AP clusters are separated in the dendrogram, whereas sequences of the Asian AP cluster comprise two separate sub-branches. The Cosmopolitan AP cluster sequences are even more diverse (three branches) (Fig 3). These differences are likely a result of the non-hierarchical clustering method of AP in contrast to phylogenetic analysis [32].
The division of RABV sequences into the two closely related branches, Asian and the Cosmopolitan AP cluster (highlighted in Fig 3), seems unclear. These sequences were all RABV variants from the Indian subcontinent, i.e. Nepal, India and Sri Lanka, which should form one joint phylogenetic cluster as suggested previously [21]. In contrast, AP allocates the Nepalese RABV isolates to the Cosmopolitan AP cluster and the Sri Lankan and Indian RABV isolates to the Asian AP cluster. To resolve this discrepancy the predefined AP cluster exemplars (New World AP cluster: JQ685974, Cosmopolitan AP cluster: KR906775, Asian AP cluster: KX148265, and Arctic AP cluster: LT598537) were required in the analyses. The resulting similarity matrix indicates that the similarities of the entire Indian subcontinent sequences to both,  set (562, A) and a reduced number of sequences (392, B). Optimal input preference for intraspecific analyses, i.e. the optimal number of clusters, was defined as the largest plateau (here four AP clusters), with the exception of the two cluster plateau (shaded gray) as methodologically, the beginning of the lower bound of the two cluster plateau cannot be defined for certain. In A the increasing length of the five cluster plateaus suggests the existence of an additional AP cluster which is not yet supported by sufficient data. the Cosmopolitan (84.86% ± 0.23%) and the Asian (84.98% ± 0.17%) AP cluster exemplars are almost equal (Table 3). Moreover, the distances of the Indian subcontinent sequences to those AP clusters are relatively large, but below the distance to both remaining cluster exemplars (New World: max 17.05%, Arctic: max 15.81%). Therefore, the analyses suggest that the degree of similarity/distance is a weak argument for allocation of these sequences into either the Cosmopolitan or Asian AP cluster. Even geographical allocation cannot resolve this problem ( Fig  2). According to these data, AP supports four main clusters. However, whether the two branches really form a separate cluster as suggested recently [21] can only be answered by including more full-genome sequences from those regions in future analyses. An advantage of AP clustering is that the dynamic evolution of further verifiable clusters can already be deduced from the input parameter iteration.

Dynamics of AP clustering
In AP, the determination of optimal numbers of clusters is a result of the previously defined iteration over all possible input parameters. To analyze the dynamics of AP in response to the number and diversity of sequences, different datasets were analyzed comprising (i) 392 fullgenome sequences ( Fig 1B) and (ii) a combined dataset including additional 170 isolates from Troupin et al. [21] (Fig 1A, S1 Table). For both datasets, four AP clusters were supported ( Fig  1A and 1B). However, when the larger dataset was used, the plateau supported an increase to five AP clusters (Fig 1A). Interestingly, the extra putative cluster represents the Africa-2 cluster previously defined by phylogenetic analysis [9,21,44] and indicated in Fig 3.   Fig 2. Global distribution of all 562 RABV full-genome sequences according to the results of AP clustering. The width of a pie chart is representing the total number of sequences from a specific country. Forty-six newly generated sequences from previously underrepresented areas in the Near East, Europe, Southern America and Africa were included in this study. The allocation to the AP clusters, i.e. New World cluster (blue), Arctic/Arctic-like (grey), Cosmopolitan (red), and Asian (green) is indicated. The nomenclature of AP clusters was based on previously assigned names. Samples from the previously described Indian subcontinent are highlighted with a red circle (Cosmopolitan sequences), and a green circle (Asian sequences).
https://doi.org/10.1371/journal.pntd.0006182.g002 A similar result is observed when partial sequences were analyzed. As mentioned above, the unclear situation for RABV sequences from the Indian subcontinent (Fig 3), i.e. Nepal, India and Sri Lanka, can be resolved by analyzing more sequences from this region. As further fullgenome sequences from this area are currently not available for these analyses, an extended panel of G-gene sequences obtained from full genome sequences (S1 Table) and additional isolates from Nepal (N = 11) and Sri Lanka (N = 50) (S2 Table) was evaluated. The analyses further supported the existence of an additional AP clusters (Fig 4). While this additional AP cluster clearly comprises sequences from the Indian subcontinent, the adjacent sixth plateau again was represented by the African-2 phylogenetic cluster. Thus, RABV may be globally separated into five or even more AP clusters if further full genome sequences are available for analysis in the future. The extended analysis of sequences highlights both the dynamics of the AP system and the robustness of cluster allocation. Of course, this is equally true for phylogenetic analysis where the length of sequence and increased geographic coverage of rabies sequences improve the analysis.   Table) and additional sequences from Nepal (N = 9) and Sri Lanka (N = 49) (S2 Table). G-gene analysis supported the existence of a fifth AP cluster as well as an additionally increased adjacent plateau. As the length of two cluster plateau cannot be defined, it is shaded in gray. https://doi.org/10.1371/journal.pntd.0006182.g004 Affinity propagation clustering for rabies

Application of AP for lyssavirus species differentiation
The lyssaviruses are a rapidly growing genus in the Rhabdoviridae family and we were interested to analyze whether AP is capable of differentiating the lyssavirus species. To this end, a set of 21 representative lyssavirus species full genome sequences were analyzed. The results confirmed the current species delineation and thus could be considered as a species demarcation by ICTV (Fig 5). Interestingly, for Lagos bat virus (LBV) AP did not segregate the sequences further, although previous studies had suggested that LBV lineage D may represent a separate lyssavirus species [49]. This analysis shows that AP clustering can provide an independent method to resolve questions regarding classification of lyssaviruses into species. However, further full genome sequences from all recognized and putative lyssavirus species would greatly facilitate these analyses.

Conclusions
The application of AP clustering, phylogenetic analyses, and the combination of both approaches revealed concordant results for RABV sub-species demarcation. The suggested approach offers advantages to phylogenetics in respect to transparency of grouping of RABV isolates and speed. As for the latter, for example the time to calculate AP cluster was about tenfold decreased compared to Neighbour Joining when applying it to the lyssavirus species sequence dataset (Fig 5). Also, the AP workflow simplifies downstream applications, such as mapping of clusters using GIS tools. In addition, we propose a combination of AP with established phylogenetic analyses to resolve phylogenetic relationships between verifiably determined clusters and sequences. This workflow could help to substantiate a transparent cluster distribution, not only for RABV but also for other comparative sequence analyses including virus species delineation as exemplarily shown for lyssaviruses (Fig 5). To this end, the cluster allocation based on AP might be implemented in phylogenetic software packages or sequence analysis pipelines (S2 Fig), and could help to substantiate transparent phylogenetic clustering as suggested previously [47]. In any case, the robustness of sequence analysis is increased by enlarging the dataset preferentially with full genome sequences from previously underrepresented areas accompanied by detailed metadata.
Supporting information S1 Table. Details of RABV full genome sequences obtained from the National Center for Biotechnology Information (NCBI) database. Datasets used to demonstrate the dynamics of AP in Fig 1A and 1B are highlighted in red. (XLSX) S2 Table. Details of partial G sequences of RABV from Nepal and Sri Lanka obtained from the National Center for Biotechnology Information (NCBI) database used to resolve the unclear situation for RABV sequences from the Indian subcontinent (Figs 3 and 4).