Mapping Strengths and Weaknesses of Different Clustering Approaches to Deciphering Bacterial Chimerism

Abstract

Bacterial genomes are chimeras of DNA of different ancestries. Deconstructing chimeric genomes is central to understanding the evolutionary trajectories of their disparate components and thus the organisms as a whole in the light of their evolutionary contexts. Of specific interest is to delineate and quantify native (vertically inherited) and alien (horizontally acquired) components of bacterial genomes and also specify genomic fractions that represent different donor sources. An agglomerative clustering procedure that prioritizes grouping of proximal similar genomic segments has previously been invoked for this purpose in conjunction with a recursive segmentation procedure. Surprisingly, however, the relative strengths and weaknesses of different clustering approaches to deciphering bacterial chimerism have not yet been investigated, despite the need to robustly interpret tens of thousands of completely sequenced bacterial genomes and nearly complete genome assemblies available in the public databases. To bridge this knowledge gap and develop more robust approaches, we assessed different clustering methods, including segment order based (proximal) clustering, hierarchical clustering, affinity propagation clustering, and a novel network clustering approach on chimeric genomes modeled after bacterial genomes representing a broad spectrum of compositional complexity. Although segment order-based clustering and network clustering compared favorably with the other approaches in discriminating between native and alien DNA at genome optimized settings, network clustering did consistently better than other methods at parametric settings optimized on all test genomes together. Segment order-based clustering and hierarchical clustering outperformed other methods in alien DNA identification while preserving donor identity in the genomes. Our study highlights the strengths and weaknesses of different approaches and suggests how this can be leveraged to achieve a more robust deconstruction of bacterial chimerism.

Introduction

Frequent horizontal gene transfer (HGT) among prokaryotes renders their genomes chimeric, with each such genome a collection of vertically transmitted and horizontally acquired genes (Abby et al, 2012; Dagan et al, 2008; Ochman et al, 2000). Deconstructing prokaryotic chimerism provides insights into the plasticity of prokaryotic genomes, which is driven by the propensity of prokaryotes to acquire and assimilate foreign genomic elements from different lineages and thus gain new traits to adapt to changes in the environment.

One of the first steps in this direction is to identify evolutionarily distinct segments in genomes. Because these disparate segments represent different genomic contexts, that is, the contexts of their source genomes, methods that invoke compositional disparity, for example, biases in (oligo)nucleotide composition or codon usage, have been developed to delineate compositionally distinct segments in genomes.

These methods use a moving window scan, gene-by-gene analysis, or a recursive segmentation approach to delineate disparate segments in genomes (Azad and Lawrence 2012; Azad et al, 2007; Braun and Muller 1998; Karlin 1998; Karlin et al, 1998; Ravenhall et al, 2015; Waack et al, 2006; Zhang and Zhang 2004). In addition to detecting the boundaries of distinct segments, methods that relate segments originating from the same source, based on their shared composition, have also been developed (Azad and Li 2013; Boys and Henderson 2004; Gionis and Mannila 2003; Keith 2006; Nicolas et al, 2002).

A combination of recursive segmentation and agglomerative clustering was demonstrated to be the most efficient among all assessed methods in recovering segments originating from different sources in chimeric genomes (Azad and Li 2013). Using this method, a genome is first recursively segmented into compositionally distinct fragments (Fig. 1A). Distinct segment types, representing potentially different sources, are identified via an agglomerative clustering procedure to group compositionally similar segments (Fig. 1B) (Azad and Li 2013; Jani and Azad 2019; Jani et al, 2016).

FIG. 1.

(A) A genome is segmented into homogeneous fragments using a recursive segmentation procedure. Compositionally similar segments can be grouped into distinct clusters using (B) segment order-based clustering, which entails comparing segments in the order of their location along the genome, from 5′-end to the 3′-end. If two segments, whether contiguous or non-contiguous, in this comparison are found similar (shown as grey boxes), they are grouped together and this process is repeated recursively from the 5′-end to the 3′-end of the genome, until none of the segments or clusters of segments can be grouped any further, (C) hierarchical clustering, which entails a pairwise comparison of segments and identification of the pair of most similar segments among all pairs. These two most similar segments form a cluster, and along with the remaining single segment clusters, the pairwise comparison is performed recursively, allowing a merger of similar clusters in this process. The recursion is halted when any two clusters are deemed significantly dissimilar to be merged. Note that both segment order-based clustering and hierarchical clustering are performed within a statistical hypothesis testing framework (see the Materials and Methods section), and (D) both network clustering and affinity propagation clustering, which are graph-based clustering algorithms that treat genomic segments as nodes connected via edges. The strength of these edges in the network clustering is inherently based on the same divergence measure that was used in segment order-based clustering and hierarchical clustering. Note, however, that the resulting edge reduction and nodal clustering are conceptually distinct in these graph-based algorithms (see the Materials and Methods section).

This step is critical for partitioning the genome into core (native, vertically inherited) and accessory (horizontally acquired) components, with the accessory further partitioned into groups with each representing a distinct donor source. The current approach prioritizes the grouping of proximal segments or clusters of segments. This is based on the premise that the DNA segments from a donor source are likely to be localized to a region and therefore could be more robustly grouped via “proximal” clustering. Although this assumption may not always hold, this approach, referred to as segment order-based clustering henceforth, has been implemented with promising results (Azad and Li 2013; Jani and Azad 2019).

In the current implementation of segment order-based clustering, segments are assessed from 5′-end to 3′-end of a given genome sequence and compositionally similar segments are grouped recursively following this sequential order of the segments along the genome (Fig. 1B). One of the limitations of this approach is that once segments or clusters are merged, they are “frozen” that is, never changed even after the retrieval of more information to refine clustering in the latter rounds of clustering. Any erroneous mergers could have cascading effects of more erroneous groupings of segments or clusters as the clustering progresses.

This could lead to sub-optimal clusters or a clustering configuration that does not adequately represent the inherent compositional structure of the genome. In an extensive benchmarking using a variety of test datasets, this approach still yielded better results than other methods, including those based on optimization techniques (Azad and Li 2013; Jani and Azad 2019). Although it compared well with other methods designed to perform similar tasks, comparisons were not performed with other clustering approaches and therefore whether the segment order-based clustering is the best approach in practice—so far as deconstructing chimeric genomes is concerned—is yet to be established.

In this study, we aimed at bridging the knowledge gaps in the usage of clustering approaches in reconstructing evolutionarily disparate components in chimeric genomes. We compared the performance of the segment order-based clustering approach with hierarchical clustering, where the grouping of two most similar segments or clusters takes precedence in a recursive process, regardless of segment locations within a genome (Fig. 1C), and with network clustering, where a network of segments (nodes representing the segments and edges signifying compositional disparity between the segments) is partitioned into clusters or modules (Fig. 1D).

Note that unlike agglomerative clustering (segment order based or hierarchical) where clusters are built bottom-up beginning with single segment clusters, network clustering is a top-down approach that partitions the network into clusters or modules (Fig. 1D). Affinity propagation, a popular clustering method often compared with network clustering (Frey and Dueck 2007; Moschopoulos et al, 2011; Vlasblom and Wodak 2009), was also tested.

In the context of deciphering genome chimerism, these clustering approaches could have strengths and weaknesses. Segment order-based clustering is expected to recover well the proximal segments originating from a source; however, its effectiveness may be diminished by its sensitivity (SN) to the segment order, as mentioned earlier. Hierarchical clustering addresses the issue of segment order by giving precedence to the most similar segments regardless of their locations in the recursive segment grouping; however, the reconstruction by this approach may be confounded by the presence of segments from phylogenetically proximal donors in a genome of interest, not a rare instance considering that preferential donors to a recipient could more likely be phylogenetically proximal than distant.

These acquired segments from phylogenetically proximal donors could be compositionally similar and may confound reconstruction by hierarchical clustering. The proposed network clustering method was designed to address the limitations of both, segment order-based clustering and hierarchical clustering, namely, the segment order issue with segment order-based clustering, and the issue with bottom-up, segment-by-segment analysis that is intrinsic to both. Previous studies have shown the effectiveness of top-down approach that can simultaneously analyze genes or segments of large structures such as genomic islands, and the network clustering is based on this premise—first construct a segment network and then partition the network in a top-down manner.

This proposed approach, and affinity propagation that generates clusters through search for “exemplars” or representative examples for the clusters, appears promising in concepts; however, their effectiveness in deciphering genome chimerism needs to be demonstrated, and further benchmarked against the segment order-based and hierarchical clustering approaches.

Assessment of the clustering methods was performed using artificial chimeric genomes constructed within a generalized hidden Markov model framework as described in Azad and Lawrence (2005). In the absence of real genomes with evolutionary events known with certainty, artificial genomes serve as a valid test platform for benchmarking different clustering approaches. Artificial chimeric genomes comprised segments with a known evolutionary history and, therefore, have been used in the past to properly assess the performance of competing algorithms (Azad and Lawrence 2005).

As real genomes are often complex with varying levels of heterogeneity, we generated artificial chimeric genomes to mimic prokaryotic genomes with varying degrees of presumed horizontal transfer events. The comprehensive set of simulated genomes, thus, provided an ideal testing platform for assessing clustering methods at different parameter settings.

In what follows, we describe clustering methods evaluated in this study and their performance on a set of artificially chimeric genomes of varying complexity. Each clustering method was tested for its ability to reconstruct segments originating from different sources in artificial chimeric genomes. The effects of parameterization were also investigated, as it is important to understand the SN of a method at different parametric settings to genome heterogeneity given that the vertically inherited and horizontally acquired fractions may vary considerably across chimeric genomes.

Materials and Methods

Genome segmentation and first-pass clustering

Segmentation of genomes is based on the usage of a generalized information-entropic measure, Markovian Jensen-Shannon divergence (MJSD), to recursively partition a sequence into compositional homogeneous segments as described (Arvey et al, 2009; Thakur et al, 2007). MJSD of order m, D^m(p₁,p₂), quantifies the difference between probability distributions p₁ and p₂ estimated from respective DNA sequences and is defined as: $D^{m} (p_{1}, p_{2}) = H^{m} (π_{1} p_{1} + π_{2} p_{2}) - π_{1} H^{m} (p_{1}) - π_{2} H^{m} (p_{2}),$

where is the Shannon entropy function for Markov model of order m, x denotes the nucleotide succeeding oligonucleotide w of length m, P(x|w) is the probability of x given the preceding oligonucleotide w, and P(w) is the probability of oligonucleotide w. ≡ {A,T,C,G}. π_i is the weight associated with p_i, π₁ + π₂ = 1.

Given a genomic sequence S of length L to be segmented into subsequences S₁ and S₂ of lengths l₁ and l₂ respectively, p_i in D^m(p₁,p₂) represents the distributions P(w) and P(x|w) in S_i, which are estimated using the frequencies of oligonucleotides w and wx in S_i. When the weight π_i is equal to l_i/L, MJSD between S₁ and S₂ can be written as: $D^{m} (S_{1} S_{2}) = H^{m} (S) - (\frac{l_{1}}{L} H^{m} (S_{1}) + \frac{l_{2}}{L} H^{m} (S_{2})),$

where S = S₁⨁S₂, and ⨁ denotes concatenation. H^m(S) or H^m(S_i) is computed using probability distributions P(w) and P(x|w) in S or S_i as described earlier. The statistical significance of a value of this measure is assessed based on the probability distribution of this measure that was shown to approximate chi-square distribution $P (D^{m} \leq X) \approx χ_{υ}^{2} (2 L (l n 2) X)$ with υ degrees of freedom (Arvey et al, 2009). The p-value is, thus, $1 - P (D^{m} \leq X)$ .

Given a sequence S of length L, the split point in S that results in maximum value of D^m(S₁,S₂) between resulting subsequences S₁ and S₂ is identified. The position with maximum D^m value, denoted $D_{m a x}^{m}$ , is determined by computing D^m for each possible split beginning with nucleotide position10 and through L-10. The statistical significance of $D_{m a x}^{m}$ is inferred based on an analytic-numerical approximation of the probability distribution of $D_{m a x}^{m}$ that was shown to follow chi-square distribution function with fitting parameters (Arvey et al, 2009; Thakur et al, 2007): $P (D_{m a x}^{m} \leq X) \approx {\{χ_{υ}^{2} [2 L (ln 2) X β]\}}^{N_{e f f}},$

where $χ_{υ}^{2}$ is the chi-square distribution function with υ degrees of freedom, and β and N_eff are fitting parameters previously estimated using Monte Carlo simulations for model order up to 2 (Arvey et al, 2009; Thakur et al, 2007). The p-value for $D_{m a x}^{m}$ is, thus, estimated as $1 - P (D_{m a x}^{m} \leq X)$ .

If the p-value for $D_{m a x}^{m}$ is below a pre-determined significance threshold, sequence S is segmented at the position of maximum divergence between S₁ and S₂. This process is recursively repeated for each resulting segment, until none of the segments can be segmented further, that is, when p-value for $D_{m a x}^{m}$ is above the significance threshold for each segment, or the segment is too short to be reliably fragmented (<16 nucleotides).

Hyper-segmentation is allowed at a relaxed stringency to identify segment boundaries with precision (PR). This generates additional splits of otherwise homogeneous segments, and therefore to remedy this, grouping of similar adjacent segments is performed within the same statistical hypothesis testing framework (Azad and Li 2013). If the p-value for MJSD between two adjacent segments is lower than an established significance threshold for contiguous segment merger, the segments are deemed significantly different; otherwise, they are merged, that is, the boundary between them is purged.

This procedure is performed recursively until no two segments sharing a boundary are deemed similar (Azad and Li 2013). Note that oligonucleotide frequencies for a cluster with multiple segments were obtained by averaging the frequencies in these segments for each oligonucleotide, and the size was obtained as the mean of segment lengths.

Markov model of order 2 (m = 2) was used in all variants of the segmentation-clustering algorithm, as parameters for statistical hypothesis testing have already been determined for the 2nd order model and previous studies have shown promising results with 2nd order model algorithms (Arvey et al, 2009; Azad and Li 2013; Jani and Azad 2019; Jani et al, 2016; Thakur et al, 2007). Both recursive segmentation and first-pass clustering were performed at the significance threshold of 0.05. This threshold was established based on our assessment on artificial chimeric genomes, as described next.

Segment order-based clustering

Location-dependent clustering has been previously implemented in several MJSD-based genome segmentation-clustering programs (Azad and Li 2013; Jani and Azad 2019; Jani et al, 2016). This approach is based on the premise that the horizontally acquired sequences from a source are often localized in a chimeric genome and therefore sequential clustering, based on the order of segments, could reconstruct such structure as this prioritizes grouping of proximal similar segments, thus minimizing the effects of potential cross-clustering of apparently similar segments from different sources.

Given N segments, this algorithm begins with N clusters containing one segment each. Beginning with the first cluster (at the 5′-end of the genome sequence), it is compared with the next cluster in order (5′-end to 3′-end) along the sequence and their compositional difference is assessed using MJSD within the same statistical framework as used for the segmentation and first-pass clustering. If the p-value for this difference is less than a preset significance threshold, the clusters are deemed compositionally distinct; otherwise, they are merged. This is performed recursively until no two clusters can be merged further, that is, the p-value for MJSD between any two clusters is less than the significance threshold.

Hierarchical clustering

Hierarchical clustering was performed to group N clusters containing a single segment each. In a pairwise manner, the compositional difference between clusters in each possible pair was assessed using MJSD. The pair of clusters with minimum MJSD was identified; if the p-value for MJSD was less than a preset significance threshold, then the clusters were deemed compositionally distinct, otherwise these clusters were merged into a single cluster resulting in total N − 1 clusters.

This procedure was followed recursively until no two clusters can be merged further. Hierarchical clustering had previously been implemented to perform clustering of genes based on codon usage patterns to localize putative alien genes in prokaryotic genomes (Azad and Lawrence 2007; Azad and Lawrence 2012).

Markov clustering algorithm-based network clustering

A network of segments was constructed, with nodes representing segments and edges signifying compositional disparity between segments. This followed segmentation of a genome and then first-pass clustering at a significance level of 0.05, as described earlier. An edge connecting nodes was established for each node pair. To each edge was associated a weight based on p-value for MJSD between nodes (segments) connecting the edge. The MJSD p-values for all node (segment) pairs were first computed and were then renormalized to represent weights on edges connecting the nodes. The weight $π_{i j}$ for the edge connecting nodes i and j was obtained as: $π_{i j} = 1 - (\frac{P (D_{i j}^{m}) - P_{m i n} (D^{m})}{P_{m a x} (D^{m}) - P_{m i n} (D^{m})})$

where $P (D_{i j}^{m})$ is the p-value for the order m MJSD between segments i and j, and $P_{m a x} (D^{m})$ and $P_{m i n} (D^{m})$ are the maximum and minimum p-values, respectively, among all segment pair MJSD p-values in the network. A graph clustering algorithm, namely, Markov clustering algorithm (MCL), was applied to partition this network into modules or clusters of segments (van Dongen 2000; van Dongen and Abreu-Goodger 2012). The MCL clustering simulates random walks within graph structures under the influence of repeated expansion and inflation.

Expansion corresponds to the random walking within the graph, establishing transition probabilities between all nodes. Inflation exaggerates these probabilities by raising all probabilities to a specified value (inflation), followed by a scaling step to ensure that the resulting matrix remains stochastic, with the eventual goal of highlighting clusters (van Dongen 2000). The mcl version 14–137 suite of software was employed for this task (https://www.micans.org/mcl/index.html?sec_software).

Tab-delimited (.abc) file containing each node pair and associated edge weight was inputted to the program that performed clustering with the “mcl” command. Clustering granularity was controlled via the inflation (-I) parameter, with tested values ranging from 1.1 to 30 in 0.1 increments. As this value is increased, the number and the specificity of clusters also increase, that is, more granular cluster configuration appears with increasing inflation parameters (van Dongen and Abreu-Goodger 2012). The output was converted to the recommended native network format of mcl using the “mcxload” and “mcxdump” commands.

Affinity-propagation clustering

A matrix of similarity values based on the p-values for MJSD between segments was constructed for use with the affinity propagation clustering method (Bodenhofer et al, 2011; Frey and Dueck 2007). The similarity matrix contained the complement of p-values for MJSD between segments for all segment pairs, resulting from segmentation and first-pass clustering at the significance level of 0.05. The similarity matrix was inputted to the affinity propagation clustering program (“apcluster” package v1.4.8 through R v4.0.0) that generated clusters of similar segments (Bodenhofer et al, 2011). The performance of the apcluster function is controlled by two primary parameters, “p” and “q,” which control the likelihood of each data sample to become a cluster representative.

By default, apcluster initializes exemplar preferences, determined by the diagonal of the matrix, for all data points. When the p parameter is unspecified, exemplar preferences are set to the quantile assigned by q, with q = 0.5 representing the median and therefore leading to default behavior (Bodenhofer et al, 2011). Although unique preferences can be set for each data point by providing a vector to the p parameter, each segment was considered equally important and the p parameter thus remained unassigned for this benchmark. In all cases, the damping factor was set to 0.90, with the “maxits” and “convits” parameters set to 10,000 and 1000, respectively.

Construction of artificial chimeric genomes

Artificial genomes for assessing the performance of clustering methods were constructed using a generalized hidden Markov model framework, as previously described (Azad and Lawrence 2005). An artificial genome was modeled after a core genome comprised a set of core (native) genes in the genome, which were extracted based on a model selection criterion algorithm at a conservative setting. To model genic variation within a core genome, the core gene set was partitioned into gene classes with distinct mutational biases using a k-means clustering algorithm with relative entropy as the distance metric.

Gene models trained on distinct gene classes were incorporated within the framework of generalized hidden Markov model to generate an artificial genome modeled after the real core genome. Artificial chimeric genomes were constructed by simulating transfer of randomly sampled genes or clusters of genes from a pool of donor artificial genomes into recipient artificial genomes. Complexity of chimeric genomes resulted from both the number of donors and the proportion of horizontally acquired genes.

For this study, 35 artificial genomes modeled after the respective prokaryotic genomes were used (Supplementary Table S1, see tab “Genome Key”). The compositions of 11 artificial chimeric genomes constructed for assessment of clustering methods are also documented in Supplementary Table S1. The artificial chimeric genomes and clustering tools used in this study have been made available at https://github.com/djburks/SGM-Clustering-Programs.

Benchmarking of clustering methods

For each clustering method, the entire parameter space was explored to find the optimal parametric setting that yields the best performance of the method. For segment order based and hierarchical clustering, the significance threshold was varied from 0.999 to 10⁻¹⁵. For MCL clustering, inflation parameter threshold was varied from 1.1 to 30 in 0.1 increment. For affinity propagation clustering, the q-parameter threshold was varied from 0.001 to 0.999 in 0.001 increment.

The accuracy metrics, namely, SN, PR, and F₁-score (F₁), were computed at each threshold setting for each clustering method. Since the vertically inherited genes are the most numerous in a prokaryotic genome, the native genome component (backbone) is identified by the largest cluster, whereas the alien component is identified by segments in the remaining smaller clusters. We established the optimal threshold settings for both, identifying native component and identifying alien component, based on F₁, which is a harmonic mean of SN and PR, computed as follows: $S N = \frac{T P}{T P + F N}$ $P R = \frac{T P}{T P + F P}$

F_{1} = 2 \cdot \frac{S N \cdot P R}{S N + P R}

In terms of identifying the native component, true positives (TP), false negatives (FN), and false positives (FP) refer to the number of native nucleotides correctly identified as native (TP), number of native nucleotides incorrectly identified as alien (FN), and number of alien nucleotides incorrectly identified as native (FP), respectively, by a method. In terms of identifying the alien component, TP, FN, and FP refer to the number of alien nucleotides correctly identified as alien (TP), number of alien nucleotides incorrectly identified as native (FN), and number of native nucleotides incorrectly identified as alien (FP), respectively.

One aspect in deconstructing a chimeric genome is the identification of native and alien components of the genome. Without multiple genome comparison, the segmentation coupled with clustering infer these components in a single genome. The success hinges on reconstructing the backbone, which has to be recovered in a single cluster by the method, and thus the native and alien components are inferred based on the cluster size—the largest being native and the rest being alien. Note that native segments often group into two or more distinct clusters and attempts to merge them by relaxing the stringency may result in undesirable mergers.

We, therefore, assessed different clustering methods on their ability to merge native clusters and thus identify the native component robustly. At the algorithmic parameter setting where F₁ of a clustering method in identifying the native component was maximized, we further assessed the ability of the method in segregating segments from different sources (native and different donors).

For each source, the clusters that contain majority of their nucleotides from the source were identified; the nucleotides within these clusters were labeled positives and all the remaining nucleotides negatives and TP, FP, and FN (with respect to the source) were computed to assess SN, PR, and overall accuracy (F₁) of a method in identifying the source. We performed source-level assessment also at the algorithmic parameter setting, where F₁ of a clustering method in identifying the alien component was maximized.

Results

Comparative assessment of clustering methods in native (core) genome isolation

Rationale for native genome identification

Identifying the native or core component of a genome is an important goal in evolutionary genomics as it enables reconstructing the tree of life based on vertically inherited genetic information, specifically the microbial clades where such signals may be obfuscated by frequent horizontal gene exchange. Using compositional bias to identify the core genome may circumvent some limitations of comparative genomics approaches that rely on the sequenced genomes of close relatives to infer the core genome.

Depending on the genomes and criteria considered, the core genome may vary significantly, and therefore using additional information such as composition in conjunction with genome comparison may help in robust extraction of core genome. For example, a native gene that displays sporadic presence in close relatives due to gene loss may gain a support for inclusion in the native group based on its compositional similarity with other (well-supported) native genes. Since the evolutionary history of nucleotides is known with certainty in artificial chimeric genomes, we assessed different clustering methods for their ability to reconstruct the core components of artificial chimeric genomes.

Parametrization for segmentation and first-pass clustering

Markov model of order 2 (m = 2) was used for both segmentation and first-pass clustering, as described in the Materials and Methods section. Note that the core genome was assembled from the segments resident in the largest cluster after segmentation and clustering (see the Materials and Methods section). Both the recursive segmentation and the first-pass clustering to group similar adjacent segments outputted by segmentation were implemented at the significance level of 0.05. This threshold was established based on extensive experiments on artificial chimeric genomes. A significance threshold of 0.05 was observed to yield fragments that consistently match well with distinct segments in artificial chimeric genomes.

Performance of each clustering method at its genome-wise optimal setting

The clustering threshold at which the overall accuracy (F₁) in identifying the core component of a genome was maximized was determined for each artificial chimeric genome, for each method. The values of SN, PR, and F₁ in identifying core genomes at the optimal settings are provided in Table 1 and the values of SN, PR, and F₁ in identifying each source (recipient and different donors) at these settings are provided in Supplementary Table S2.

Table 1.

Sensitivity, Precision, and F₁ Score Values in Identifying Native Nucleotides in Test Genomes by Different Clustering Methods at the Genome-Specific Optimal Parameter Settings Where F₁ Scores Are Maximal

Test genome	Segment-order				Hierarchical				Network (MCL)				Affinity propagation
Test genome	Parameter (T3)	SN	PR	F1	Parameter (T3)	SN	PR	F1	Parameter (I)	SN	PR	F1	Parameter (q)	SN	PR	F1
Artificial chimeric genome 1	1.00E-08	0.999	0.996	0.997	1.00E-04	0.999	0.990	0.995	29.7	0.987	0.998	0.992	0.863	0.460	0.997	0.629
Artificial chimeric genome 2	0.019	0.997	0.962	0.980	0.034	0.996	0.966	0.981	27.4	0.987	0.971	0.979	0.712	0.507	0.973	0.667
Artificial chimeric genome 3	0.367	0.994	0.966	0.980	0.095	0.995	0.964	0.979	25.9	0.984	0.974	0.979	0.797	0.538	0.878	0.667
Artificial chimeric genome 4	0.449	0.974	0.959	0.966	0.012	0.994	0.951	0.972	27.2	0.994	0.951	0.972	0.887	0.406	0.962	0.571
Artificial chimeric genome 5	0.558	0.993	0.920	0.955	0.42	0.994	0.882	0.935	21.9	0.977	0.931	0.953	0.378	0.516	0.934	0.665
Artificial chimeric genome 6	0.62	0.987	0.922	0.953	0.863	0.987	0.916	0.950	21.2	0.970	0.934	0.952	0.769	0.420	0.916	0.576
Artificial chimeric genome 7	0.7	0.971	0.867	0.916	0.55	0.974	0.847	0.906	19.3	0.957	0.867	0.910	0.845	0.398	0.726	0.514
Artificial chimeric genome 8	0.313	0.988	0.981	0.985	0.018	0.995	0.971	0.983	23.6	0.987	0.982	0.984	0.83	0.458	0.980	0.624
Artificial chimeric genome 9	0.005	0.985	0.993	0.989	0.001	0.984	0.994	0.989	29	0.995	0.989	0.992	0.3	0.388	0.993	0.558
Artificial chimeric genome 10	0.029	0.985	0.997	0.991	1.00E-04	0.983	0.998	0.990	30	0.992	0.997	0.995	0.296	0.421	0.999	0.592
Artificial chimeric genome 11	1.00E-07	0.993	0.982	0.987	0.174	0.967	0.991	0.979	27.7	0.995	0.987	0.991	0.881	0.468	0.992	0.636

F1, F₁-score; MCL, Markov clustering algorithm; PR, precision; SN, sensitivity.

All clustering methods except affinity propagation clustering attained F₁ over 90% across all test genomes (Table 1). Segment order-based clustering and network clustering produced th highest overall accuracy (F₁ averaged over 11 test genomes: 0.972); average F₁ values for hierarchical clustering and affinity propagation clustering were 0.969 and 0.609, respectively. In genome-wise assessment, segment order-based clustering was found to have maximum F₁ for 6 of 11 genomes, followed by networking clustering (three genomes) and hierarchical clustering (two genomes).

As expected, we found that the optimal performance varies with genome composition. Although all clustering methods except affinity propagation clustering did well on genomes with substantial contribution from each of fewer donors (e.g., Artificial Chimeric Genome 1 in Supplementary Table S1), the performance declined with increasing complexity, particularly with increasing number of donors with fewer gene contributions by each (e.g., Artificial Chimeric Genome 7 in Supplementary Table S1).

Trend in overall performance of each clustering method across its entire threshold range and determination of threshold yielding maximal overall accuracy

In our aforementioned assessment, we observed that performance was optimized at different thresholds for different genomes by each method. Thus, a default threshold that can work well across genomes of different compositions cannot be established. However, it could still be possible to determine a threshold or a threshold range where the overall performance is acceptable across genomes of different compositions. We, therefore, computed the performance metrics averaged over all test genomes at each threshold (varied from 10⁻¹⁵ to 0.999 for segment order-based clustering and hierarchical clustering, from 1.1 to 30 (inflation) for network clustering, and 0.001 to 0.999 (q-value) for affinity propagation clustering).

The overall accuracy as assessed by F₁ measure is shown as a function of threshold for each method (Fig. 2). As expected, the best performance of both segment order-based clustering and hierarchical clustering was attained when the significance threshold was close to 0 (toward the rightmost end of F₁ plots in Fig. 2: 0.019 for the former and 0.002 for the latter; thresholds are represented here as a percentage of the maximum allowable threshold (complement of significance level for hierarchical and segment order clustering)).

FIG. 2.

Overall accuracy (F₁ averaged over all test genomes) in identifying native (core) nucleotides for each clustering method at increasing thresholds. The network clustering threshold (MCL inflation parameter) is represented as the percentage of the maximum allowable threshold of 30. The affinity propagation threshold is represented as a percentage of the maximum allowable threshold of 1. Agglomerative clustering (hierarchical and segment order) thresholds are represented here as a percentage of the maximum allowable threshold of 1–10⁻¹⁵ (complement of significance level). F₁, F₁ score; MCL, Markov clustering algorithm.

These are more relaxed stringencies for clustering, allowing the merger of multiple native clusters into a single native cluster. This is apparent in the SN plot (Supplementary Fig. S1), which shows that relaxing stringency (lower p-values) continually increases SN as the significance threshold approaches 0 (toward the rightmost end in Supplementary Fig. S1); however, the PR declines remarkably as alien clusters begin to coalesce with the native cluster as the stringency is successively relaxed (PR plot, Supplementary Fig. S2).

The overall performance of hierarchical clustering was on average lower than that of segment order-based clustering across the entire threshold range (Fig. 2), though this difference becomes much smaller at a significance threshold 0.001 or lower (Fig. 3). Notably, the F₁ varied in a stepwise manner with a threshold for both these clustering approaches, more so with hierarchical clustering (Fig. 2). F₁ varied in steps because of the fact that a significant increase in recall happened only when a threshold value was reached where small native cluster(s) merged with the largest cluster.

FIG. 3.

Overall accuracy (F₁ averaged over all test genomes) in identifying native (core) nucleotides for the hierarchical and segment order-based clustering methods at MJSD significance level range of 10⁻¹ to 10⁻¹⁵. MJSD, Markovian Jensen-Shannon divergence.

The stepwise pattern is more pronounced with hierarchical clustering where the largest (native) cluster remains invariant over longer threshold intervals until a smaller native cluster merges as the stringency is successively relaxed. The “punctuated” native cluster merger resulting in a ladder step pattern is less apparent with segment order-based clustering, where proximal native segments or clusters are continually merged as encountered from 5′ to 3′ end, recursively (Supplementary Fig. S1).

In contrast, network clustering displayed a different pattern of F₁ variation with threshold compared with agglomerative (segment order based or hierarchical) clustering (Fig. 2). In contrast to gradual stepwise variation in F₁ with threshold in agglomerative clustering, the F₁ variation with inflation parameter is shaped like a bucket with a rather abrupt transition from high to low or low to high F₁. The optimal performance (maximal F₁) was observed at inflation parameter values between 20 and 23 for artificial chimeric genomes.

Unlike agglomerative clustering where typically a single large native cluster is formed after successive rounds of cluster merger, analysis of cluster formation by network clustering revealed concurrent formation of at least two large native clusters that coalesce at the optimal threshold (typically between 20 and 23), yielding maximal F₁.

All methods except for affinity propagation attained high average F₁ (>0.9) in identifying native nucleotides, with network clustering with maximal average F₁ of 0.963 at threshold 23.3 slightly outperforming segment order-based clustering and hierarchical clustering that attained maximal average F₁ of 0.957 and 0.950 at significance thresholds 0.019 and 0.002, respectively (Fig. 2, Table 2). The maximal average F₁ of affinity propagation was achieved with the q parameter set to 0.831, but it was much lower than those of the other three methods at only 0.574.

Table 2.

Sensitivity, Precision, and F₁ Score Values in Identifying Native Nucleotides in Test Genomes by Different Clustering Methods, Generated at the Optimal Parametric Setting of Each Method Where the F₁ Score Averaged Over All Test Genomes Was Maximal

Test genome	Segment-order				Hierarchical				Network (MCL)				Affinity propagation
Test genome	Parameter (T3)	SN	PR	F1	Parameter (T3)	SN	PR	F1	Parameter (I)	SN	PR	F1	Parameter (q)	SN	PR	F1
Artificial chimeric genome 1	0.019	0.987	0.996	0.992	0.002	0.999	0.990	0.995	23.3	0.975	0.998	0.986	0.831	0.450	0.997	0.620
Artificial chimeric genome 2	0.019	0.997	0.962	0.980	0.002	0.997	0.958	0.978	23.3	0.983	0.975	0.979	0.831	0.488	0.972	0.650
Artificial chimeric genome 3	0.019	0.998	0.946	0.971	0.002	0.997	0.948	0.972	23.3	0.979	0.975	0.977	0.831	0.534	0.877	0.664
Artificial chimeric genome 4	0.019	0.983	0.923	0.952	0.002	0.995	0.927	0.960	23.3	0.980	0.960	0.970	0.831	0.360	0.841	0.504
Artificial chimeric genome 5	0.019	0.984	0.885	0.932	0.002	0.996	0.875	0.931	23.3	0.983	0.906	0.943	0.831	0.504	0.933	0.655
Artificial chimeric genome 6	0.019	0.994	0.852	0.918	0.002	0.995	0.852	0.918	23.3	0.982	0.911	0.945	0.831	0.317	0.962	0.477
Artificial chimeric genome 7	0.019	0.992	0.753	0.856	0.002	0.995	0.655	0.790	23.3	0.994	0.734	0.844	0.831	0.376	0.719	0.494
Artificial chimeric genome 8	0.019	0.995	0.945	0.970	0.002	0.996	0.931	0.962	23.3	0.987	0.982	0.984	0.831	0.434	0.981	0.602
Artificial chimeric genome 9	0.019	0.984	0.991	0.988	0.002	0.984	0.994	0.989	23.3	0.984	0.994	0.989	0.831	0.320	0.994	0.484
Artificial chimeric genome 10	0.019	0.985	0.995	0.990	0.002	0.983	0.998	0.990	23.3	0.984	0.998	0.991	0.831	0.408	0.999	0.579
Artificial chimeric genome 11	0.019	0.980	0.989	0.984	0.002	0.996	0.943	0.968	23.3	0.979	0.990	0.985	0.831	0.412	0.992	0.582

For agglomerative and network clustering, high performance in identifying native DNA was attained at a threshold range that allowed the merger of several native clusters into a single native cluster without incurring undesirable mergers (none or very few alien and native cluster mergers); these were inflation parameter range 20–23 for network clustering and significance level range 0.01–0.001 for agglomerative clustering. This was not the case for affinity propagation, where although the PR remained high across all tested thresholds (>0.9) (Supplementary Fig. S2), the SN was consistently low (< 0.54) and declined sharply after reaching a maximal average of 0.538, resulting in a similar trend observed with F₁ (Supplementary Fig. S1).

Genome-wise performance shows that the overall accuracies of the methods declined (Table 2, Supplementary Table S3) compared with the performance at genomic-specific optimal settings (Table 1, Supplementary Table S2), as expected. Network clustering produced maximum F₁ for 8 out of 11 genomes, outperforming other methods (Table 2). Segment order-based clustering, hierarchical clustering, and affinity propagation clustering produced maximum F₁ for 2, 1, and 0 genome(s), respectively (Table 2).

In Supplementary Table S3, we provide average SN, PR, and F₁ in identifying sources (recipient and different donors) at the same threshold setting that yielded maximal average F₁ in identifying native DNA for each method (Fig. 2). Network clustering whose maximal average F₁ in identifying native DNA was highest among all four methods identified sources more robustly than other methods (identified 71 sources with F₁ > 0.7, whereas the hierarchical, segment order-based, and affinity propagation clustering methods identified 69, 68, and 59 sources, respectively, with F₁ > 0.7).

Notably, F₁ for native source was greater for the former compared with the latter three clustering approaches. This demonstrates that a better discrimination of native from alien is attained via robust grouping of segments from different sources, particularly by robustly grouping the native segments into a single large native cluster in this particular instance.

Comparative assessment for alien genome identification

Rationale for alien genome identification

Quantifying alien component of prokaryotic genomes is central to understanding how prokaryotes evolve through HGT and acquire novel traits, for example, virulence, antibiotic resistance, or novel metabolic capabilities that help them adapt to new conditions or environments. Robust delineation and aggregation of alien segments is critical to achieving this goal. Therefore, similar to native DNA identification, we compared the four clustering methods to determine the approach and algorithm parameter setting that yield the maximal overall accuracy in alien DNA identification. Note that the approach and parameter setting that yielded maximal overall accuracy in native DNA identification may not necessarily represent the optimal framework for alien DNA identification as well.

Performance of each clustering method at its genome-wise optimal setting

First, at each threshold setting, the largest cluster was identified as native and the remaining clusters as alien (the native nucleotides being most numerous coalesce into a large cluster whereas the alien nucleotides coalesce into several small clusters each representing a distinct donor source, consistent with the composition of prokaryotic genomes previously reported). SN, PR, and F₁ for alien nucleotide detection were computed and the values for each artificial chimeric genome at the optimal setting (where F₁ was maximized) were tabulated (Table 3).

Table 3.

Sensitivity, Precision, and F₁ Score Values in Identifying Alien Nucleotides in Test Genomes by Different Clustering Methods at the Genome-Specific Optimal Parameter Settings Where F₁ Scores Are Maximal

Test genome	Segment-order				Hierarchical				Network (MCL)				Affinity propagation
Test genome	Parameter (T3)	SN	PR	F1	Parameter (T3)	SN	PR	F1	Parameter (I)	SN	PR	F1	Parameter (q)	SN	PR	F1
Artificial chimeric genome 1	1.00E-08	0.988	0.916	0.951	0.0001	0.984	0.821	0.895	29.7	0.799	0.959	0.872	0.863	0.088	0.975	0.161
Artificial chimeric genome 2	0.019	0.983	0.788	0.875	0.025	0.976	0.810	0.885	23.3	0.904	0.864	0.884	0.712	0.257	0.923	0.402
Artificial chimeric genome 3	0.367	0.973	0.864	0.915	0.095	0.979	0.856	0.913	25.9	0.934	0.897	0.915	0.551	0.328	0.954	0.489
Artificial chimeric genome 4	0.452	0.924	0.886	0.904	0.012	0.982	0.858	0.916	25	0.953	0.881	0.915	0.876	0.368	0.960	0.532
Artificial chimeric genome 5	0.722	0.979	0.809	0.886	0.275	0.983	0.706	0.821	21.4	0.939	0.842	0.888	0.378	0.461	0.919	0.614
Artificial chimeric genome 6	0.962	0.973	0.853	0.909	0.738	0.973	0.842	0.903	21.2	0.943	0.880	0.910	0.769	0.477	0.932	0.631
Artificial chimeric genome 7	0.997	0.970	0.863	0.914	0.937	0.973	0.838	0.900	19.3	0.957	0.865	0.908	0.968	0.562	0.980	0.714
Artificial chimeric genome 8	0.313	0.964	0.943	0.953	0.018	0.983	0.911	0.946	23.6	0.960	0.944	0.952	0.83	0.374	0.972	0.541
Artificial chimeric genome 9	0.005	0.958	0.978	0.968	0.001	0.954	0.982	0.968	29	0.986	0.967	0.977	0.003	0.355	0.992	0.523
Artificial chimeric genome 10	0.029	0.956	0.990	0.973	0.0001	0.952	0.993	0.972	30	0.978	0.991	0.984	0.296	0.368	0.999	0.538
Artificial chimeric genome 11	1.00E-07	0.977	0.946	0.961	0.043	0.906	0.973	0.938	27.7	0.984	0.961	0.973	0.881	0.378	0.988	0.547

Network clustering produced the highest F₁ in identifying alien nucleotides for 5 of 11 artificial chimeric genomes, displaying a better performance than the other clustering approaches; the average F₁ of 0.928 by segment order-based clustering was, however, ∼0.003 higher than network clustering with the highest alien nucleotide F₁ for 4 of 11 artificial chimeric genomes. This is primarily due to the performance discrepancy between network clustering and segment order-based clustering in detecting alien nucleotides in Artificial Chimeric Genome 1, the chimeric genome with the least amount of alien nucleotides (5%).

Hierarchical clustering with average F₁ of 0.914 generated the highest alien nucleotide F₁ for 2 of 11 artificial chimeric genomes; affinity propagation clustering with average F₁ of 0.517 was outperformed by the other methods on both aggregate and individual genome tests. The values of these metrics in identifying nucleotides from each source (recipient and different donors) within the same optimized settings are provided in Supplementary Table S4. The percentages (and relative percentages) of native and donor nucleotides in each cluster at these and other optimal settings are illustrated in Supplementary Figures S3–S6.

Trend in overall performance of each clustering method across its entire threshold range and determination of threshold yielding maximal overall accuracy

Similar to the analysis for native nucleotide identification, we also assessed the F₁ for alien DNA identification averaged over all test genomes at each threshold for each method (Fig. 4). SN, PR, and F₁ for each genome at the setting where average F₁ for alien DNA identification was maximized is provided in Table 4 for each method and the values of these metrics for each source in test genomes at the same setting are provided in Supplementary Table S5. At its optimal setting, network clustering generated F₁ (averaged over 11 genomes) of 0.903 outperforming other clustering methods (average F₁ for the segment order-based, hierarchical, and affinity propagation clustering methods were 0.881, 0.864, and 0.499, respectively; Table 4).

FIG. 4.

Overall accuracy (F₁ averaged over all test genomes) in identifying alien nucleotides for each clustering method at increasing thresholds. The network clustering threshold (MCL inflation parameter) is represented as the percentage of the maximum allowable threshold of 30. The affinity propagation threshold is represented as a percentage of the maximum allowable threshold of 1. Agglomerative clustering (hierarchical and segment order) thresholds are represented as a percentage of the maximum allowable threshold of 1–10⁻¹⁵ (complement of significance level).

Table 4.

Sensitivity, Precision, and F₁ Score Values in Identifying Alien Nucleotides in Test Genomes by Different Clustering Methods, Generated at the Optimal Parametric Setting of Each Method Where the F₁ Score Averaged Over All Test Genomes Was Maximal

Test genome	Segment-order				Hierarchical				Network (MCL)				Affinity propagation
Test genome	Parameter (T3)	SN	PR	F1	Parameter (T3)	SN	PR	F1	Parameter (I)	SN	PR	F1	Parameter (q)	SN	PR	F1
Artificial chimeric genome 1	0.019	0.934	0.793	0.858	0.002	0.821	0.983	0.895	23.7	0.960	0.782	0.862	0.753	0.976	0.078	0.145
Artificial chimeric genome 2	0.019	0.788	0.983	0.875	0.002	0.766	0.982	0.861	23.7	0.858	0.903	0.880	0.753	0.926	0.253	0.397
Artificial chimeric genome 3	0.019	0.779	0.988	0.871	0.002	0.787	0.986	0.875	23.7	0.902	0.925	0.913	0.753	0.947	0.322	0.480
Artificial chimeric genome 4	0.019	0.774	0.944	0.850	0.002	0.782	0.984	0.871	23.7	0.885	0.943	0.913	0.753	0.805	0.312	0.450
Artificial chimeric genome 5	0.019	0.716	0.953	0.818	0.002	0.684	0.987	0.808	23.7	0.761	0.958	0.848	0.753	0.916	0.454	0.607
Artificial chimeric genome 6	0.019	0.697	0.985	0.816	0.002	0.695	0.989	0.816	23.7	0.829	0.966	0.892	0.753	0.932	0.475	0.629
Artificial chimeric genome 7	0.019	0.700	0.990	0.820	0.002	0.516	0.991	0.679	23.7	0.632	0.992	0.772	0.753	0.891	0.551	0.681
Artificial chimeric genome 8	0.019	0.974	0.954	0.964	0.002	0.982	0.954	0.968	23.7	0.983	0.956	0.969	0.753	0.992	0.346	0.513
Artificial chimeric genome 9	0.019	0.984	0.956	0.970	0.002	0.993	0.952	0.972	23.7	0.993	0.954	0.973	0.753	0.999	0.365	0.535
Artificial chimeric genome 10	0.019	0.828	0.982	0.899	0.002	0.778	0.984	0.869	23.7	0.944	0.960	0.952	0.753	0.975	0.363	0.529
Artificial chimeric genome 11	0.019	0.967	0.939	0.953	0.002	0.815	0.984	0.892	23.7	0.970	0.948	0.959	0.753	0.989	0.359	0.527

Although these values are less than those from a genome-wise optimal setting (Table 2) as expected, they are still close and just within 3% for network clustering highlighting its ability to provide a stable parametric setting for application to genomes of variable composition in identifying the alien nucleotides. Notably, network clustering was found to have maximum F₁ for 9 of 11 genomes.

Assessment of the clustering methods in inferring alien DNA sources

Although identifying alien nucleotides could serve well the purposes of some studies, the goal of deciphering clusters arising from different donor sources may not be well served if the focus is on only native and alien discrimination. For example, a method may allow alien clusters from different sources to merge native clusters to achieve maximal F₁ in alien nucleotide identification. Therefore, in our second way in assessment, we evaluated the ability of the clustering methods in inferring alien DNA sources by identifying clusters that harbor segments from each donor source.

If DNA sequences from a source organism were assigned to multiple clusters, only those clusters with majority of resident nucleotides belonging to the source were deemed as representing the source organism. As artificial genomes are modeled after core genomes, we expect a method to generate as many clusters as sources contributing to an artificial chimeric genome. Therefore, we selected the largest cluster among clusters representing a source organism, and we then computed SN, PR, and F₁ in identifying each donor source for each artificial chimeric genome.

The values of these metrics, averaged over the donors, indicate the ability to identify alien DNA in a genome while preserving the identity of the donors in the process. These values at the optimal settings of clustering methods are provided for each genome in Table 5, and the respective values for each source in each genome are provided in Supplementary Table S6. Note that in the event that a donor was not identified as the majority representative of a single cluster, the F₁ for that particular donor was considered zero.

Table 5.

Sensitivity, Precision, and F₁ Score Values in Identifying Alien Nucleotides While Preserving Donor Identity in the Test Genomes by Different Clustering Methods, Generated at the Optimal Parametric Setting of Each Method Where the F₁ Score Averaged Over All Test Genomes Was Maximal

Test genome	Segment-order				Hierarchical				Network (MCL)				Affinity propagation
Test genome	Parameter (T3)	SN	PR	F1	Parameter (T3)	SN	PR	F1	Parameter (I)	SN	PR	F1	Parameter (q)	SN	PR	F1
Artificial chimeric genome 1	0.453	0.587	0.981	0.730	0.332	0.596	0.978	0.738	1.4	0.664	0.974	0.778	0.776	0.590	0.980	0.734
Artificial chimeric genome 2	0.948	0.588	0.935	0.686	0.747	0.568	0.959	0.683	22.6	0.540	0.897	0.646	0.874	0.487	0.928	0.612
Artificial chimeric genome 3	0.561	0.615	0.850	0.696	0.987	0.609	0.851	0.695	27.8	0.490	0.779	0.571	0.869	0.495	0.798	0.578
Artificial chimeric genome 4	0.937	0.526	0.809	0.616	0.914	0.559	0.864	0.646	28.1	0.500	0.806	0.587	0.946	0.482	0.891	0.611
Artificial chimeric genome 5	0.999	0.481	0.847	0.572	0.999	0.498	0.912	0.612	22.5	0.479	0.749	0.563	0.914	0.469	0.832	0.572
Artificial chimeric genome 6	0.998	0.423	0.828	0.502	0.999	0.458	0.851	0.536	20.7	0.403	0.668	0.480	0.942	0.406	0.786	0.508
Artificial chimeric genome 7	0.991	0.418	0.826	0.462	0.998	0.401	0.874	0.455	18.7	0.365	0.582	0.392	0.955	0.353	0.746	0.448
Artificial chimeric genome 8	0.313	0.719	0.924	0.794	0.216	0.649	0.937	0.752	1.7	0.607	0.949	0.723	0.851	0.522	0.960	0.673
Artificial chimeric genome 9	0.754	0.684	0.925	0.761	0.173	0.692	0.920	0.764	2.2	0.564	0.964	0.700	0.753	0.507	0.982	0.666
Artificial chimeric genome 10	0.373	0.711	0.866	0.763	0.902	0.710	0.824	0.746	4.2	0.579	0.940	0.700	0.889	0.519	0.959	0.668
Artificial chimeric genome 11	0.796	0.603	0.922	0.703	0.043	0.604	0.905	0.689	1.4	0.605	0.891	0.678	0.873	0.476	0.982	0.640

The cluster representing each donor is the largest cluster harboring segments primarily of that donor.

Segment order-based clustering produced the highest F₁ scores for a majority of genomes (6 of 11 genomes), whereas hierarchical, network, and affinity clustering methods produced highest F₁ on 4, 1, and 0 genome(s) respectively. The values of F₁ averaged over 11 genomes were similar for segment order-based clustering and hierarchical clustering (0.662 and 0.665, respectively), and they were higher than those for network clustering (0.619) and affinity propagation clustering (0.609). Further, we computed SN, PR, and F_1, averaged over all genomes, at different threshold settings for each method (Fig. 5, Supplementary Figs. S7 and S8).

FIG. 5.

Overall accuracy (F₁ averaged over all test genomes) in identifying alien nucleotides while preserving donor identity in the genomes for each clustering method at increasing thresholds. The network clustering threshold (MCL inflation parameter) is represented as the percentage of the maximum allowable threshold of 30. The affinity propagation threshold is represented as a percentage of the maximum allowable threshold of 1. Agglomerative clustering (hierarchical and segment order) thresholds are represented as a percentage of the maximum allowable threshold of 1–10⁻¹⁵ (complement of significance level). The cluster representing each donor is the largest cluster harboring segments primarily of that donor. For donors that did not make up the majority of nucleotides in any cluster, the values of metrics SN, PR, and F1 were deemed zero. PR, precision; SN, sensitivity.

F₁ of all four methods varied within 5% for over 90% of their threshold range (Fig. 5). The maximal average F₁ values were 0.623, 0.615, 0.576, and 0.563 for hierarchical clustering, segment order-based clustering, network clustering, and affinity propagation clustering, respectively. Thus, overall, agglomerative clustering methods outperformed other clustering approaches in alien DNA identification while preserving donor identity in the genomes.

Notably, the thresholds that yielded maximal F₁ in native and alien DNA detection differ noticeably for agglomerative clustering methods; in contrast, there was an overlap of the optimal threshold ranges for native and alien DNA detection by network clustering. Affinity propagation clustering, as with native genome isolation, performed worse than the other methods.

Discussion

Our analysis highlights the promises and challenges of different clustering approaches in unraveling bacterial chimerism. Although segment order-based clustering was introduced as a part of an integrative approach to segment and cluster bacterial DNA sequences, it is for the first time that hierarchical, network, and affinity propagation clustering techniques were assessed along-side segment order-based clustering in this context. Although hierarchical clustering was earlier employed for gene clustering to identify compositionally atypical genes (Azad and Lawrence, 2007), it was not tested on variable length genomic fragments outputted by a recursive segmentation program.

Diametrically opposite to agglomerative clustering is network clustering; it partitions the network to identify clusters rather than builds bottom up as with agglomerative clustering. Network clustering brought a new perspective in addressing this problem, as this is conceptually different from often invoked agglomerative clustering. The top-down approach that performs partitioning or segmentation of an organized dataset (e.g., a genome or a network) was earlier found to be more effective in localizing horizontally acquired regions in bacterial genomes in comparison to bottom-up methods such as agglomerative clustering (Arvey et al, 2009).

The horizontally acquired genomic segments were earlier identified by assessing the atypicality of each segment against the genome background (Arvey et al, 2009). Later studies showed that clustering after segmentation yields even more promising results; however, only the agglomerative clustering approach was tested (Azad and Li, 2013). The problem with the bottom-up (agglomerative) approach is that members of earlier formed clusters are not allowed to be reassigned at later rounds of clustering (except for the merger of two clusters into a new cluster) when additional information emerges that may rectify any past misassignments.

If clusters are not refined at successive steps, impurities in the clusters may only amplify as clustering proceeds. Both segment order-based and hierarchical clustering approaches are prone to be affected by this. Although reassignment is itself a non-trivial problem, another issue that afflicts segment order-based clustering is that clustering configuration varies by how the grouping is initiated (i.e., from 5′-end or 3′-end of a genome). This is not unexpected, as earlier rounds of clustering dictate the later rounds and thus the directionality impacts the clustering. Hierarchical clustering, however, proceeds independent of directionality.

Deciphering the underlying structure using agglomerative clustering has its pitfalls; however, this may still yield reasonably good outcomes. This was realized while probing genomic islands in bacterial genomes by previous studies that led to the development of a top-down, recursive partitioning approach that allowed assessment of entire information content to successively generate compositionally homogeneous segments (Arvey et al, 2009). Here, we explored whether a top-down rather than a bottom-up approach could again be invoked to more efficiently reconstruct the chimeric structure of bacterial genomes after segmentation.

We posited that this could be feasible within a network framework that will allow partitioning the network of segments into modules or clusters of similar segments. Implementation of networking clustering and comparative assessment revealed its promising aspects; at the same time, it also highlighted the complementarity of different approaches.

Agglomerative clustering performed well in discriminating between native and alien DNA at genome optimized thresholds, whereas the optimal threshold settings varied among genomes. Considering that optimizing a clustering program on a just sequenced, anonymous genome may not be feasible, we assessed the clustering methods for their performance averaged over the test genomes across the entire threshold range for each method. Our intent was to identify a threshold range for each method where its maximal or close to maximal overall performance could be attained. We observed the network clustering's overall performance maximizing within a certain inflation parameter range (averaging 25.7 ± 3.7 for native and 25.1 ± 3.7 for alien DNA identification, Figs. 2 and 4 and Tables 1, 2), reaching a maximal higher than by other methods for both native and alien identification.

This demonstrates that top-down partitioning post-segmentation, indeed, yields a better outcome. However, the accuracy is lower over a large range than agglomerative clustering before it begins to spike toward maximal that appeared more stable with native than alien DNA identification. Agglomerative clustering's performance is relatively less variable over the entire threshold range and reaches the maximal at low p-values as expected (right end of hierarchical clustering and segment-order clustering curves in Figs. 2 and 4).

These results highlight the strengths and weaknesses of different approaches, and future studies could focus on exploiting the complementary strengths of different clustering approaches to develop an integrative method for native and alien DNA identification. For example, as agglomerative clustering's performance is less variable and is optimal in low p-value range, one can first use agglomerative clustering to get a cue of cluster structure inherent to a genome of interest and then follow up with network clustering to “fine-tune” to optimal or nearly optimal cluster configuration.

This will help exploit the benefits of the network approach while obviating low accuracy risk as could occur within certain inflation parameter ranges with network clustering (Figs. 2 and 4). This is supported by our observation that with network clustering, the merger of large native clusters happened within the inflation parameter range of 16–23 across all artificial chimeric genomes yielding maximal F₁, and therefore fine tuning can be performed within this range.

Identifying native component of a genome is an important goal, and the alternative approach suggested here will complement the frequently used comparative genomics approach in establishing backbone genomes in different taxa (Chung et al, 2018; Daubin et al, 2002; Na et al, 2018; Segata and Huttenhower 2011). Our results showed that the overall accuracy (F₁) in identifying native DNA in an ensemble of artificial chimeric genomes could reach up to 0.972, with network clustering as the best overall performer.

Genome-wise optimal performers for native DNA identification varied among agglomerative and network clustering methods. Using both in concert, as suggested earlier, may enhance the reliability of clustering in inferring the core component of bacterial genomes. In addition to estimating native and alien nucleotides in bacterial genomes, estimating donors is another problem where the power of clustering can be leveraged upon. Clustering methods may be expected to group segments from each donor within a unique cluster but a donor genome may be chimeric and therefore multiple clusters may also be expected for a donor.

Here, as we used genomes modeled after the respective cores of the bacterial genomes, we expect methods to group segments from a donor into a single cluster. Note that the genomes were modeled on cores that were obtained by removing compositionally atypical genes; a core, thus, represents the mutational proclivity of the recipient genome with the composition shaped by the directional mutational pressure specific to that genome. Because of the superior performance of agglomerative clustering in identifying alien DNA while preserving donor identities in the major clusters of the donors, one can use this for more reliably estimating the donors.

It must also be emphasized that the performance of clustering methods exploiting compositional disparity in delineating native and alien DNAs declines with increasing complexity of chimeric genomes. For example, when the number of donors for the artificial chimeric Escherichia coli genome was increased to 34, with each contributing ∼1.5% resulting in the recipient E. coli genome comprising ∼50% alien DNA (Artificial Chimeric Genome 7 in Tables 1–5), the overall accuracy (F₁) declines in general.

This is expected as a lower number of donors with greater contribution by each will provide substantial proportion of DNA by each and thus a better statistics aiding alien DNA quantification by the clustering methods. Note that even when challenged with test genomes of a high level of complexity, the best performing clustering methods (agglomerative/network) could still attain overall accuracy in ∼68–92% range for native or alien DNA identification (Tables 1–4), indicating their potential utility in interpreting even complex bacterial genomes.

In summary, agglomerative and network clustering approaches possess complementary strengths that can be exploited to understand the different facets of genome chimerism, which may not be possible with any single method. Future efforts could focus on integrating different approaches toward the goal of still better interpretation of bacterial genomes.

Footnotes

Author Disclosure Statement

The authors declare they have no conflicting financial interests.

Funding Information

The authors did not receive funds from any funding organization for this work.

Supplementary Material

Abbreviations Used

References

Abby

, Tannier

, Gouy

, et al. Lateral gene transfer as a support for the tree of life. Proc Natl Acad Sci U S A, 2012; 109:4962–4967. doi: 10.1073/pnas.1116871109.

Arvey

, Azad

, Raval

, et al. Detection of genomic islands via segmental genome heterogeneity. Nucleic Acids Res. 2009; 37:5255–5266. doi: 10.1093/nar/gkp576.

Azad

, Lawrence

. Use of artificial genomes in assessing methods for atypical gene detection. PLoS Comput Biol, 2005; 1:e56. doi: 10.1371/journal.pcbi.0010056.

Azad

, Lawrence

. Detecting laterally transferred genes: Use of entropic clustering methods and genome position. Nucleic Acids Res, 2007; 35:4629–4639. doi: 10.1093/nar/gkm204.

Azad

, Lawrence

. Detecting laterally transferred genes. Methods Mol Biol, 2012; 855:281–308. doi: 10.1007/978-1-61779-582-4_10.

Azad

, Lawrence

, Thakur

, Ramaswamy

Segmentation of genomic DNA sequences. In: Advanced Computational Methods for Biocomputing and Bioimaging. (Pham TD, Yan H, Crane DI. eds.) Nova Science Publishers: New York; 2007.

Azad

, Li

Interpreting genomic data via entropic dissection. Nucleic Acids Res, 2013; 41:1–21.

Bodenhofer

, Kothmeier

, Hochreiter

APCluster: An R package for affinity propagation clustering. Bioinformatics, 2011; 27:2463–2464. doi: 10.1093/bioinformatics/btr406.

Boys

, Henderson

. A Bayesian approach to DNA sequence segmentation. Biometrics, 2004; 60:573–581; discussion 581–588. doi: 10.1111/j.0006-341X.2004.00206.x.

10.

Braun

, Muller

H-G

. Statistical methods for DNA sequence segmentation. Stat Sci, 1998; 13:142–162.

11.

Chung

, Munro

, Tettelin

, et al. Using core genome alignments to assign bacterial species. mSystems, 2018; 3:e00236-18. doi: 10.1128/mSystems.00236-18.

12.

Dagan

, Artzy-Randrup

, Martin

Modular networks and cumulative impact of lateral transfer in prokaryote genome evolution. Proc Natl Acad Sci, 2008; 105:10039–10044. doi: 10.1073/pnas.0800679105.

13.

Daubin

, Gouy

, Perrière

A phylogenomic approach to bacterial phylogeny: Evidence of a core of genes sharing a common history. Genome Res, 2002; 12:1080–1090. doi: 10.1101/gr.187002.

14.

Frey

, Dueck

Clustering by passing messages between data points. Science (80-), 2007; 315:972–976. doi: 10.1126/science.1136800.

15.

Gionis

, Mannila

Finding recurent sources in sequences. In: Proceedings of the Seventh Annual International Conference on Computational Molecular Biology—RECOMB'03. ACM Press: New York, New York, USA; 2003; pp. 123–130. doi: 10.1145/640075.640091.

16.

Jani

, Azad

. IslandCafe: Compositional anomaly and feature enrichment assessment for delineation of genomic islands. G3 (Bethesda) 2019; 9:3273–3285. doi: 10.1534/g3.119.400562.

17.

Jani

, Mathee

, Azad

. Identification of novel genomic islands in liverpool epidemic strain of pseudomonas aeruginosa using segmentation and clustering. Front Microbiol, 2016; 7:1210. doi: 10.3389/fmicb.2016.01210.

18.

Karlin

Global dinucleotide signatures and analysis of genomic heterogeneity. Curr Opin Microbiol, 1998; 1:598–610. doi: 10.1016/s1369-5274(98)80095-7.

19.

Karlin

, Mrázek

, Campbell

. Codon usages in different gene classes of the Escherichia coli genome. Mol Microbiol, 1998; 29:1341–1355. doi: 10.1046/j.1365-2958.1998.01008.x.

20.

Keith

JM.

Segmenting eukaryotic genomes with the Generalized Gibbs Sampler. J Comput Biol, 2006; 13:1369–1383. doi: 10.1089/cmb.2006.13.1369.

21.

Moschopoulos

, Pavlopoulos

, Iacucci

, et al. Which clustering algorithm is better for predicting protein complexes?. BMC Res Notes, 2011; 4:549. doi: 10.1186/1756-0500-4-549.

22.

S-I

, Kim

, Yoon

S-H

, et al. UBCG: Up-to-date bacterial core gene set and pipeline for phylogenomic tree reconstruction. J Microbiol, 2018; 56:280–285. doi: 10.1007/s12275-018-8014-6.

23.

Nicolas

, Bize

, Muri

, et al. Mining Bacillus subtilis chromosome heterogeneities using hidden Markov models. Nucleic Acids Res, 2002; 30:1418–1426. doi: 10.1093/nar/30.6.1418.

24.

Ochman

, Lawrence

, Groisman

. Lateral gene transfer and the nature of bacterial innovation. Nature 2000; 405:299–304. doi: 10.1038/35012500.

25.

Ravenhall

, Škunca

, Lassalle

, et al. Inferring horizontal gene transfer. PLoS Comput Biol, 2015; 11:e1004095. doi: 10.1371/journal.pcbi.1004095.

26.

Segata

, Huttenhower

Toward an efficient method of identifying core genes for evolutionary and functional microbial phylogenies. Gibas C (ed.). PLoS One, 2011; 6:e24704. doi: 10.1371/journal.pone.0024704.

27.

Thakur

, Azad

, Ramaswamy

Markov models of genome segmentation. Phys Rev E Stat Nonlinear Soft Matter Phys, 2007; 75. doi: 10.1103/PhysRevE.75.011915.

28.

van Dongen

(2008). Graph Clustering Via a Discrete Uncoupling Process. SIAM Journal on Matrix Analysis and Applications,, 30(1), 121–141. doi: 10.1137/040608635.

29.

van Dongen

, Abreu-Goodger

Using MCL to extract clusters from networks. Methods Mol Biol, 2012; 804:281–295. doi: 10.1007/978-1-61779-361-5_15.

30.

Vlasblom

, Wodak

. Markov clustering versus affinity propagation for the partitioning of protein interaction graphs. BMC Bioinformatics, 2009; 10:99. doi: 10.1186/1471-2105-10-99.

31.

Waack

, Keller

, Asper

, et al. Score-based prediction of genomic islands in prokaryotic genomes using hidden Markov models. BMC Bioinformatics, 2006; 7:142. doi: 10.1186/1471-2105-7-142.

32.

Zhang

, Zhang

C-T

. A systematic method to identify genomic islands and its applications in analyzing the genomes of Corynebacterium glutamicum and Vibrio vulnificus CMCP6 chromosome I. Bioinformatics, 2004; 20:612–622. doi: 10.1093/bioinformatics/btg453.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

5.85 MB

0.00 MB

0.19 MB

0.21 MB

0.17 MB

0.01 MB

0.05 MB

0.06 MB

0.05 MB