Characterizing Scales of Genetic Recombination and Antibiotic Resistance in Pathogenic Bacteria Using Topological Data Analysis

Pathogenic bacteria present a large disease burden on human health. Control of these pathogens is hampered by rampant lateral gene transfer, whereby pathogenic strains may acquire genes conferring resistance to common antibiotics. Here we introduce tools from topological data analysis to characterize the frequency and scale of lateral gene transfer in bacteria, focusing on a set of pathogens of significant public health relevance. As a case study, we examine the spread of antibiotic resistance in Staphylococcus aureus. Finally, we consider the possible role of the human microbiome as a reservoir for antibiotic resistance genes.


Introduction
Pathogenic bacteria can lead to severe infection and mortality and presents an enormous burden on human populations and public health systems. One of the achievements of twentieth century medicine was the development of a wide range of antibiotic drugs to control and contain the spread of pathogenic bacteria, leading to vastly increased life expectancies and global economic development. However, rapidly rising levels of multidrug antibiotic resistance in several common pathogens, including Escherichia coli, Klebsiella pneumoniae, Staphylococcus aureus, and Neisseria gonorrhoea, is recognized as a pressing global issue with near-term consequences [15,23,26]. The threat of a post-antibiotic 21st century is serious, and new methods to characterize and monitor the spread of resistance are urgently needed.
Antibiotic resistance can be acquired through point mutation or through horizontal transfer of resistance genes. Horizontal exchange occurs when a donor bacteria transmits foreign DNA into a genetically distinct bacteria strain. Three mechanisms of horizontal transfer are identified, depending on the route by which foreign DNA is acquired [16]. Foreign DNA can be acquired via uptake from an external environment (transformation), via viral-mediated processes (transduction), or via direct cell-to-cell contact between bacterial strains (conjugation).
Resistance genes can be transferred between strains of the same species, or can be acquired from different species in the same environment. While the former is generally more common, an example of the latter is the phage-mediated acquisition of Shiga toxin in E. coli in Germany in 2011 [18]. Elements of the bacterial genome that show evidence of foreign origin are called genomic islands, and are of particular concern when associated with phenotypic effects such as virulence or antibiotic resistance.
The presence of horizontal gene transfer precludes accurate phylogenetic characterization, because different segments of the genome will have different evolutionary histories. Bacterial species definitions and taxonomic classifications are made on the basis of 16S ribosomal RNA, a highly conserved genomic region between bacteria and archaea species [25]. However, the region generally accounts for less than 1% of the complete genome, implying that the vast majority of evolutionary relationships are not accounted for in the taxonomy [5]. Because of the important role played by lateral gene transfer, new ways of characterizing evolutionary and phenotypic relationships between microorganisms are needed.
Topological data analysis (TDA), and persistent homology in particular, has been shown to be an effective tool for capturing horizontal evolutionary processes at the population level by measuring deviations from treelike additivity. Initial work in this direction characterized recombination in viral evolution, particularly in influenza, where genomic reassortment can lead to the emergence of viral pandemics [3]. Further work established foundations for statistical inference in population genetics models using TDA [8]. We provide a brief overview of TDA and persistent homology in Section 2, for additional reviews see [2] [7].
In this paper we explore topics relating to horizontal gene transfer in bacteria and the emergence of antibiotic resistance in pathogenic strains. We show that TDA can not only quantify gene transfer events, but also characterize the scale of gene transfer. The scale of recombination can be measured from the distribution of birth times of the H 1 invariants in the persistent barcode diagram. It has been shown that recombination rates decrease with increasing sequence divergence [9]. We characterize the rate and scale of intraspecies recombination in several pathogenic bacteria of public health concern. We select a set of pathogenic bacteria that are of significant public health interest based on a recently released World Health Organization (WHO) report on antimicrobial resistance [26]. Using persistent homology, we characterize the rate and scale of recombination in the core genome using multilocus sequence data. To extend our characterization to the whole genome, we use protein family annotations as a proxy for sequence composition. This allows us to compute a similarity matrix between strains. Comparing persistence diagrams gives us information about the relative scales of gene transfer at arbitrary loci. The species selected for study and the sample sizes in each analysis are specified in Table 1. Next, we explore the spread of antibiotic resistance genes in S. aureus using Mapper, an algorithm for partial clustering and visualization of high dimensional data [19]. We identify two major populations of S. aureus, and observe one cluster with strong enrichment for the antibiotic resistance gene mecA. Importantly, resistance appears to be increasingly spreading in the second population. Finally, we consider the risk of horizontal transfer of resistance genes from the human microbiome into an antibiotic sensitive strain, using β-Lactam resistence as an example. In this environment, benign bacterial strains can harbor known resistance genes. We use a network analysis to visualize the spread of antibiotic resistance gene mecA into nonnative phyla. Each individual has a unique microbiome, and we speculate that microbiome typing of this sort may useful in developing personalized antibiotic therapies. These results demonstrate the important role the HCI-KDD approach can play in tackling the challenges of large scale -omics data applied to clinical settings and personalized medicine: Interactive visualization through graph and network construction, data mining global invariants with topological algorithms, and knowledge discovery through data integration and fusion [11][10].

Topological Data Analysis and Persistent Homology
Topological data analysis computes global invariants from point cloud data. These global invariants represent loops, holes, and higher dimensional voids in data. A topological representation of the data is constructed by building a set of triangulated objects representing the connectivity of the data at different scales, called a filtration. Various constructions exist for triangulating data. The most efficient approach for large scale data is the Vietoris-Rips complex, which associates a simplex to a set of points if they are pairwise connected. In this way, the complex is specified purely by its 1-skeleton, which can be efficiently computed.
Given a filtration, persistent homology is an algorithm to associate homology groups to each scale, which give information about the invariants in the data. H 0 gives information about the connectivity, H 1 about loops, etc. The output of the algorithm is a set of intervals corresponding to topological features present at different scales. The homology information can be compactly summarized as a barcode diagram, in which invariants are represented as horizontal line segments with the birth and death scales as the left and right edge of the line segment, respectively. Alternatively, the homology information can be represented as a persistence diagram, a 2-D plot in which intervals are represented as points on a (birth, death) plane.
For the purposes of studying biological sequence data, and horizontal evolutionary processes in particular, these approaches are widely useful, as it was shown in [3] that sequence datasets with treelike phylogeny will have vanishing higher homology. The observed homological features are therefore direct evidence for horizontal exchange amonst the sequences in a sample. These topological constructions become a natural way of reasoning about evolutionary relationships between organisms in cases where treelike phylogeny is not appopriate.

Evolutionary Scales of Recombination in the Core Genome
The genes comprising the bacterial genome can be largely paritioned into two groups: the core genome, consisting of those genes that are highly conserved and characteristic of a given species classification, and the accessory genome, consisting of those genes whose presence can be variable even within strains of the same species. We first sought to examine scales of recombination in the core bacterial genome using multilocus sequence typing (MLST) data. MLST is a method of rapidly assigning a sequence profile to a sample bacterial strain. For each species, a predetermined set of loci on a small number of housekeeping genes are selected as representative of the core genome of the species. As new strains are sequenced, they are annotated with a profile corresponding to the sequence type at each locus. If a sample has a previously unseen type at a given locus, it is appended to the list of types at that locus. Large online databases have curated MLST data from labs around the world; significant pathogens can have several thousand typed strains (over 10,000 in the case of Neisseria spp.).
Because different species will be typed at different loci, examining direct interspecies genetic exchange with this data is unfeasible, however MLST provides a large quantity of data with which to examine intraspecies exchange in the core genome. However, because the selected loci are generally all housekeeping genes, this type of recombination analysis will be only informative about genetic exchange in the core genome. Mobile genetic elements will have separate rates of exchange.
We investigate horizontal exchange in the core genome for twelve pathogens using MLST data from PubMLST [13]. For each strain, a pseudogenome can be constructed by concatenating the typed sequence at each locus. Using a Hamming metric, we construct a pairwise distance matrix between strains and compute persistent homology on the resulting metric space. Because of the large number of sample strains, we employ a Lazy Witness complex with 250 landmark points and ν = 0 [6]. The computation is performed using javaplex [21].
An example of our output is shown in Figure 1, where we plot the H 1 barcode diagrams for K. pneumoniae and S. enterica. The two species have distinct recombination profiles, characterized by the range of recombinations: K. pneumoniae recombines at only one short-lived scale, while S. enterica recombines both at the short-lived scale and a longer-lived scale. We repeat this analysis for each species, and plot the results as a persistence diagram in Figure 2  We define a relative rate of recombination by counting the number of H 1 loops across the filtration and dividing by the number of samples for that species. The results are shown in Figure 3, where we observe that different species can have vastly different recombination profiles. For example, S. enterica and E. coli have the highest recombination rates, while H. pylori is substantially lower than the others. Coupled with the smaller scale of recombinations suggests that the H. pylori core genome is relatively resistant to recombination except within closely related strains.

Protein Families as a Proxy for Genome Wide Reticulation
Protein family annotations cluster proteins into sets of isofunctional homologs, i.e., clusters of proteins with both similar sequence composition and similar function. A particular strain is represented as a binary vector indicating the presence or absence of a given protein family. Correlations between strains can reveal genome-wide patterns of genetic exchange, unlike the MLST data which can only provide evidence of exchange in the core genome. We use the FigFam protein annotations in the Pathosystems Resource Institute Center (PATRIC) database because of the breadth of pathogenic strain coverage and depth of genomic annotations [24]. The FigFam annotation scheme consists of over 100,000 protein families curated from over 950,000 unique proteins [14]. For each strain we compute a transformation into FigFam space. We transform into this space because the frequency of genome rearrangements and differences in mobile genetic elements makes whole genome alignments unreliable, even for strains within the same species. As justification for performing this step, it has been shown experimentally that recombination rates decrease with increasing genetic distance [9]. After transforming, we construct a strain-strain correlation matrix and compute the persistent homology in this space. In Figure 4 we show the persistence diagram relating the structure and scale between different species. We find that different species have a much more diverse topological structure in this space than in MLST space, and a wide variety of recombination scales. The large scales of exchange in H. influenzae suggest it can regularly acquire novel genetic material from distantly related strains.

Antibiotic Resistance in Staphylococcus aureus
S. aureus is a gram positive bacteria commonly found in the nostrils and upper respiratory tract. Certain strains can cause severe infection in high-risk populations, particulary in the hospital setting. The emergence of antibiotic resistant S. aureus is therefore of significant clinical concern. Methicillin resistant S. aureus (MRSA) strains are resistant to β-Lactam antibiotics including penicillin and cephalosporin. Resistance is conferred by the gene mecA, an element of the Staphyloccoccal cassette chromosome mec (SCCmec). mecA codes for a dysfunctional penicillin-binding protein 2a (PBP2a), which inhibits β-Lactam antibiotic binding, the primary mechanism of action [12]. Of substantial clinical importance are methods for characterizing the spread of MRSA within the S. aureus population.
To address this question, we use the FigFam annotations in PATRIC, as described in the previous section. PATRIC contains genomic annotations for 461 strains of S. aureus, collectively spanning 3,578 protein families. We perform a clustering analysis using the Mapper algorithm as implemented in Ayasdi Iris [1]. Principal and second metric singular value decomposition are used as filter functions, with a 4x gain and an equalized resolution of 30. This results in a graph structure with two large clusters, connected by a narrow bridge, as shown in Figure 5. The two clusters are consistent with previous phylogenetic studies using multilocus sequence data to identify two major population groups [4]. Of the 461 S. aureus strains in PATRIC, 142 carry the mecA gene. When we color nodes in the network based on an enrichment for the presence of mecA, we observe a much stronger enrichment in one of the two clusters. This suggests that β-Lactam resistance has already begun to dominate in that clade, likely due to selective pressures. More strikingly, we observe that while mecA enrichment is not as strong in the second cluster, there is a distinct path of enrichment emanating along the connecting bridge between the two clusters and into the less enriched cluster. This suggests the hypothesis that antibiotic resistance has spread from the first cluster into the second cluster via strains intermediate to the two, and will likely continue to be selected for in the second cluster.

Microbiome as a Reservoir of Antibiotic Resistance Genes
While antibiotic resistance can be acquired through gene exchange between strains of the same species, it is also possible for gene exchange to occur between distantly related species. It has been recognized that an individual's microbiome, the set of microorganisms that exist symbiotically within a human host, can act as a reservoir of antimicrobial resistance genes [20,17]. It is of substantial clini- To address this question, we use data from the Human Microbiome Project (HMP), a major research initiative performing metagenomic characterization of hundreds of healthy human microbiomes [22]. The HMP has defined a set of reference strains that have been observed in the human microbiome. We collect FigFam annotations from PATRIC for the reference strain list in the gastrointestinal tract. We focus on the gastrointestinal tract because it is an isolated environment and likely to undergo higher rates of exchange than other anatomic regions. Of the 717 gastrointestinal tract reference strains, 321 had FigFam annotations. We computed a similarity matrix as in previous sections, using correlation as distance. The resulting network is shown in Figure 6, where strains are colored by phyla-level classifications. While largely recapitulating phylogeny, the network depicts interesting correlations between phyla, such as the loop between Firmicutes, Bacteroides, and Proteobacteria. Next, we searched for genomic annotations relating to β-Lactam resistance. 10 strains in the reference set had matching annotations, and we highlight those strains in the network with green diamonds. We observe resistance mostly concentrated in the Firmicutes, of which S. aureus is a member, however there is a strain of Proteobacteria that has acquired the resistance gene. Transfer of β-Lactam resistance into the Protebacteria is clinically worrisome. Pathogenic Proteobacteria include S. enterica, V. cholerae, and H. pylori, and emergence of β-Lactam resistance will severely impact currently used antibiotic drug therapies.
The species composition of each individual's microbiome can differ substantially due to a wide variety of poorly understood factors [22]. In this case, an individual's personal microbiome network will differ from the network we show in Figure 6, which was constructed from the set of all strains that have been reported across studies of multiple individuals. The relative risk for acquiring self-induced resistance will therefore vary from person to person and by the infectious strain acquired. However, a network analysis of this type can be used to assess risk and give clues as to possible routes by which antibiotic resistance may be acquired. In the clinical setting, this could assist in developing personalized antibiotic treatment regimens. We propose a more thorough expansion of this work, examining the full range of antibiotic resistance genes in order to quantify microbiome risk factors for treatment failure. We foresee an era of genomically informed infectious disease management in the clinical setting, based on an understanding of a patient's personal microbiome network profile.
In this paper we have brought some ideas from topological data analysis to bear on problems in pathogenic microbial genetics. First, we used persistent homology to evaluate recombination rates in the bacterial core genomes using MLST profile data. We showed that different pathogens have different recombination rates. We expanded this to gene transfer across whole bacterial genomes by using protein family annotations in the PATRIC database. We found different scales of recombination in different pathogens. Next, we explored the spread of MRSA in S. aureus populations using topological methods. We identified two major population clusters of S. aureus, and noted increasing resistance in a previously isolated population. Finally, we studied the emergence of β-Lactam resistance in the microbiome, and proposed methods by which personal risk could be assessed by microbiome typing. Each stage of this analysis represents a successful application of the HCI-KDD approach to biomedical discovery. Our results point to the important role for graph mining and topological data mining in health and personalized medicine.