Here we analyze four datasets, three of which we introduced earlier (Cong et al., 2016b): 20 Escherichia coli and seven Shigella genomes (ECS dataset), 110 enteric bacterial genomes (EB) and 143 genomes from BA. To these we now add a dataset of 144 bacterial genomes (BAC) purpose-built for this analysis. When this latter dataset was constructed, 24 bacterial orders in 12 classes were represented by at least one genus from which at least six genomes had been sequenced to high quality. Within each of these orders we selected one genus at random, and if that genus was represented by more than six genomes we chose six of them at random, thereby constituting BAC with 144 genomes in 24 genera. In this way we attempt to achieve as broad and balanced selection of genomes across Bacteria as possible, given the available data, a synthetic classification (NCBI) and the underlying biology.

As noted above, TF-IDF infers transfers from identified groups of donor genomes into a single recipient genome. It is therefore necessary to delineate groups prior to analysis. Here we recognize groups within the ECS dataset according to multi-locus sequence type (MLST; Gordon et al., 2008); within EB by genus, sometimes combining Escherichia and Shigella genomes into a single group; within BA by phylum, or alternatively by class; and within BAC by order. Other approaches to grouping are possible, some of which we explored earlier (Cong et al., 2016b).

The TF-IDF method per se proceeds in four steps, as follows: for each dataset we (A) extract unique k-mers and construct a k-mer dictionary; and (B) build a relationship matrix R in which rows represent individual genomes, columns represent the identified groups of genomes, and elements count the number of identical k-mers present in a genome and in each group other than its own. These counts are then normalized, and the mean element value computed over R. Unless indicated otherwise, this mean is used as the threshold for recognizing that a genome in the dataset may contain k-mers donated by a group in the dataset. (C) Within each genome, we then construct segments from neighboring k-mers that are present in the same donor group. We further merge these segments if they are separated by less than a gap threshold G, yielding potential lateral segments. (D) If the average frequency of k-mers in a potential lateral segment is less than that of all k-mers in the target genome’s own group, then we consider it an inferred lateral segment. Step B implements the IDF component, and step D the TF component (Cong et al., 2016a,b). Pseudocode is available in the Supplementary Material to Cong et al. (2016a), and the TF-IDF source code at https://github.com/congyingnan/TF-IDF.git.

Because the TF component requires a potential lateral segment to be infrequent in genomes of its own group, TF-IDF is expected to identify recent LGT events, i.e., those affecting one or a very few genomes in a target group. By contrast, k-mers descendant from transfers more ancient than the common ancestor of a target group would tend to occur widely within that group, and thus fail to be inferred as lateral.

We consider a gene to be lateral if it contains, or is overlapped by, at least one inferred lateral segment such that two distinct length thresholds are met. The inferred lateral segment must itself contain at least a specified minimum number of k-mers (including k-mers in any intervening gaps up to G = 2k); this minimum number is 10 for the BA dataset, 100 for EB, 500 for ECS (Cong et al., 2016b) and 10 for BAC. These values approximate the average length of all LGT detections in each dataset, thereby controlling in part for differences in sequence diversity among the datasets. In addition, the overlap must extend for at least a specified minimum number of k-mers (again including k-mers in gaps up to G = 2k); this minimum number is 10 for BA, 100 for ECS and EB (Cong et al., 2016b) and 10 for BAC.

As introduced above, GECs might variously be described as paths, transitively closed sets, paracliques or cliques (Skippington and Ragan, 2011). The first two structures fail to capture the density of connectivity, and many such structures of nearly equivalent size or value can often be found in relatively highly connected graphs such as the LGT networks we derive above. Paracliques differ from the corresponding cliques by relaxing the strict requirement that all edges be present, and in this way might better ameliorate the effects of incomplete or imperfect data. In the absence of theory or established practice, paraclique parameters would have to be optimized for each dataset, requiring intense computation. Constrained by these considerations, here we adopt the strictest yet clearest definition of GEC, as a set of vertices that share (donate or receive) genetic material from all other nodes within this set. That is, there must be at least one direct path between each node and every other node. Using this definition, GECs correspond to cliques in the LGT network.

The discovery and analysis of such cliques proceeds in four main steps:

From our previous work (Cong et al., 2016a,b) we know that k can strongly affect the detection of LGT, hence potentially the topologies of LGT networks. For that reason, we explore different values of k to test the stability of clique topology. In step (a) we explore values of k from 20 to 40. Depending on the data, false positives can predominate at k ≤ 20, while at k ≥ 40 shared k-mers become too rare, resulting in diminished performance. For consistency with earlier studies on the ECS, EB and BA datasets (Cong et al., 2016b), gap size G was fixed at 2k. The step size is 10 for the ECS and EB datasets, while for BA (where LGT signal is much weaker) and BAC (not studied heretofore) we set step size as 5 for improved resolution against k.

Our TF-IDF procedure infers LGT from a donor group to a recipient sequence, so at this point the vertices in our inferred networks are of two types: individual genomes when they are recipients of LGT, and groups of genomes when they are donors. Of course, members of a group may individually be (and often are) recipients of LGT from outside that group. Edges are directional, so we depict them using an arrow from donor to recipient (Figure 1).

We aim to delineate GECs, introduced above as sets of nodes that have both donated genetic material to, and received genetic material from, each other. However, it is unclear how to extract these relationships when nodes are of two types, individual sequences and groups. For this reason we take only groups as nodes in our network analysis. In step (b) we subsume each genome into its respective own-group identity, and merge all directed edges into those sequences into a single directed edge from the donor group to the recipient group (Figure 1). The integer weight on each edge gives the total number of genes inferred in this way to have been affected by LGT from donor to recipient groups.

Extracting cliques [step (c)] is known to be NP-hard (Karp, 1972), although a parameterized complexity approach (Abu-Khzam et al., 2006) is possible whenever the clique size can be bounded (Downey et al., 1999). For this we used the Graph Algorithms Pipeline for Pathway Analysis (GrAPPA) software suite (Langston Lab, the University of Tennessee¹). GrAPPA integrates multiple graph-theoretical tools for biological data analysis, including those designed to find cliques and paracliques. It implements tools to extract patterns efficiently from graphs, but deals only with undirected and unweighted networks. For this reason, we reformulated our directed networks as undirected graphs (i.e., we disregarded the arrowhead), deleted all weight annotations (number of inferred LGT genes) on each edge, and merged edges between pairs of nodes that are both donors and recipients. Such reformulation does not make full use of the LGT information (e.g., directionality) provided by TF-IDF, but nonetheless it preserves connectivity information sufficient for discovery of GECs as currently defined.

Here we use GrAPPA to find, for each dataset, one maximum clique (a clique with the greatest number of vertices) and all maximal cliques (cliques that are not included in a larger clique). We report only those maximal cliques with at least three vertices.

In step (d) we determine the biological processes that are enriched in the cliques previously extracted. As we show below (Results), the cliques inferred for the ECS and EB datasets encompass almost all the respective lateral genes, so the biological process enrichments are essentially the same as described previously (Cong et al., 2016b). Most LGT inferred for BAC involves the EB genera. Thus, here we report biological process enrichment only for cliques inferred, at different values of k, for the BA dataset. These genes are extracted from the GenBank genome record using GI numbers and coordinates, and collected as a test set. All genes in the dataset form the respective reference set. The enrichment statistic is a Fisher’s exact test, for which we set false discovery rate FDR = 0.05 as the threshold for selecting over- and under-represented Gene Ontology (Ashburner et al., 2000; Gene Ontology Consortium, 2004) terms.

We divide the ECS dataset (20 Escherichia coli and seven Shigella genomes) into six groups according to MLST (Gordon et al., 2008). In an earlier analysis of this dataset (Skippington and Ragan, 2012), lateral events identified by topological incongruence between trees inferred from individual putative orthogroups and an MRP (Ragan, 1992) reference supertree were shown to be biased more by phylogeny than by environment or lifestyle; concern was also expressed that defining GECs as cliques or paracliques might be too rigorous a standard.

Here, we use our TF-IDF method to infer LGT networks (Figure 2). For all k examined here, all six phyletic groups belong to a single clique, so the whole dataset forms one large GEC. Indeed, at k = 30 or 40, topologies of the two networks are identical (as before: Cong et al., 2016b). There is a clear trend overall of more detections on each edge as k increases, but with some exceptions: at k = 20 we find three edges not seen at k = 30, from group D to B2 (257 transfers), from B2 to S (443) and from B2 to E (3574), while transfers from D to B1 decrease from 4659 at k = 20 to 3200 at k = 30. For all other edges, more genes are affected by LGT at k = 30 than at k = 20. Likewise, when k is increased from 30 to 40, three edges show fewer detections (D to B1, 3200 to 1842; B1 to B2, 3658 to 3563; E to B2, 3804 to 3363) but all others have more.

Although clique topology is stable for 20 ≤k ≤ 40, the total number of lateral genes underlying each clique increases with k (Figure 3). This increase might appear to contradict our earlier finding that when k increases, the total number of detections and detection length should remain the same or decrease (at G = 2k). However, when k is small, more short segments tend to be detected as lateral (Supplementary Table 1). For example, at k = 20, 26% of lateral segments are ≥500 k-mers in length, our threshold for selecting the segments for mapping to genes. This proportion increases to 31% at k = 40. Thus, we infer more lateral segments of ≥500 k-mers at k = 40, which leads to more genes being inferred as affected by LGT.

The enteric bacteria dataset contains 110 genome sequences from five genera: Escherichia, Shigella, Salmonella, Klebsiella and Yersinia. Because the delineation of groups affects the detection results (Cong et al., 2016b), we infer LGT networks and extract cliques from five variants of this dataset: all genera present (referred to as EB-1); all genera except Shigella (EB-2) or alternatively, all except Escherichia (EB-3); with Escherichia and Shigella combined into a single group (EB-4); and with both Escherichia and Shigella removed (EB-5). These are five of the six variants we examined earlier (Cong et al., 2016b).

If we keep all 110 sequences and group them by genus (EB-1 dataset), the LGT network topologies change as k steps from 20 to 30 to 40. At k = 20, Escherichia, Shigella, Klebsiella constitute a single clique. At k = 30, we find two cliques, one consisting of Escherichia and Shigella, the other of Escherichia and Klebsiella. At k = 40 only one clique is found, consisting of Escherichia and Shigella. We infer many more LGT events between Escherichia and Shigella than between any other pair of genera. As Escherichia and Shigella are present in the clique across the examined range of k, we can say that they are the core nodes of this GEC.

Because genomes from Escherichia and Shigella share many more identical k-mers than do other groups, the lateral signal between these genera can drown out weaker lateral signal from or between other genera. This happens because the IDF values (elements of the R matrix) for these genomes are much higher than for the others (Cong et al., 2016a). This pushes up the IDF threshold, with the consequence that few lateral events are detected involving the other genera. To explore this effect, we also analyzed variant datasets which are modified so that Escherichia and Shigella do not both appear in the dataset as separate genera.

We first removed the Shigella genomes from the dataset while retaining those from Escherichia, thereby eliminating the effect of Shigella (EB-2 dataset). We now infer additional lateral events in both directions between all pairs of Escherichia, Salmonella, and Klebsiella. Thus we find a GEC composed of Escherichia, Klebsiella and Salmonella that remains stable with respect to k. We find similar results when we instead retain Shigella sequences while removing those of Escherichia (EB-3 dataset); the GEC here is Shigella, Klebsiella and Salmonella. We do infer lateral events between Klebsiella and Yersinia in EB-3, but these are not sufficient for Yersinia to join the GEC. In EB-4 we combine Escherichia and Shigella into a single group (ES); more lateral events were inferred from Salmonella to ES, but the GEC membership remains ES, Salmonella and Klebsiella. Lastly, to eliminate completely the effects of Escherichia and Shigella on LGT inference, we use only Klebsiella, Salmonella and Yersinia as input (EB-5). At k = 20 the sole clique contains all three genera, but at k = 30 or 40 the previous clique is split into two, one containing Klebsiella and Salmonella and the other Klebsiella and Yersinia. We thus conclude that Escherichia, Shigella, Klebsiella and Salmonella are all members of a larger GEC. Details are provided in Table 1.

The BA dataset, 143 genome sequences across BA, has been studied in our group using classical alignment-based and other computational methods for more than a decade (Beiko et al., 2005; Chan et al., 2009a,b). Like many empirical datasets it is unbalanced, with many more genomes representing some taxa (e.g., Proteobacteria, Firmicutes) than others. We group these genomes into fifteen phyla or, alternatively, into 31 classes. With more nodes than in the two previous datasets, there is potential for inferred LGT networks to be more complex. On the other hand, these genomes are more dissimilar to each other (Cong et al., 2016b), so fewer k-mers are shared and fewer instances of LGT are inferred.

When groups are delineated by phylum, the number of total LGT detections decreases significantly as k increases (Figure 4), and this causes edges in the LGT network to vanish and the cliques to shrink. At the smallest value of k = 20 six maximal cliques are found, each with five phyla. Five of these contain the High G+C Firmicutes, Proteobacteria and Low G+C Firmicutes, which together represent 14797 lateral genes, 95.5% of the total inferred over the entire network. Thus these phyla form the core of the inter-phylum GEC. We also observe a smaller GEC of Nanoarchaeota, Euryarchaeota and Crenarchaeota; although based on only 10 lateral genes, it is notable for showing potential GECs among Archaea. In addition, the Thermus/Deinococcus phylum contributes 244 lateral events, 1.5% of the total; as our dataset contains only one strain in this phylum, this particular genome appears to be more LGT-active than many other bacterial genomes.

The number of detections drops sharply at k > 20; recall that our earlier simulations (Cong et al., 2016a) indicate potential false positives at k ≤ 20, presumably due to identical k-mers shared between sequences and groups simply by coincidence. As k increases and LGT detections decrease in number, some edges in the LGT network vanish, but the core nodes – the High G+C Firmicutes, Low G+C Firmicutes and Proteobacteria – remain as members of the maximal cliques. Thermus/Deinococcus also remains active in sharing LGT with Proteobacteria for all investigated k.

When these genomes are alternatively grouped by class, the LGT networks are more complex. Again we see a sharp drop in detections for k ≥ 20. At k = 20, all but one of the 31 classes are involved in LGT (30696 genes), and we observe 23 maximal cliques (≥3 nodes) in the LGT network; however, five classes form core members of the GEC, with each being present in 17 maximal cliques (≥5 nodes) and in the maximum clique. These classes are the Actinomycetales (5377 genes with lateral origin), Bacillus/Clostridium (2277) and the α- (5944), β- (7322) and γ-Proteobacteria (8596). Together they contain 77.7% of all genes that contain regions of inferred lateral origin.

Since the sequences within BA are relatively dissimilar from each other, many fewer k-mers are shared between sequences than in the ECS and EB datasets. Thus the LGT detections are very sensitive to k (Figure 4). At k = 30 the γ- and β-Proteobacteria, Actinomycetales and Bacillus/Clostridium are hubs and play key roles in most cliques; at k = 40 fewer genes are inferred as lateral, and only the former two classes remain as the core.

Deinococcus is inferred to exchange genetic material with β- and γ-Proteobacteria at 20 ≤k ≤ 40. Lateral events are also inferred between Deinococcus and Actinomycetales, and between Deinococcus and Chlorococcales, at k < 40.

With the BAC dataset we again explore a broad phyletic range (24 orders representing 12 classes across Bacteria); but unlike the situation with BA (above), with BAC we maintain numerical balance (six genomes per order) and a comparable degree of local sequence diversity (each set of six genomes represents a single genus) to the extent possible, given the underlying biology and the availability of high-quality genome sequences. As simulations (Cong et al., 2016a) indicate a high likelihood of false positive detections at k = 20, here we vary k from 25 to 40 in steps of 5. As above, the total number of genes inferred to be affected by LGT events decreases with increased k (Figure 5).

At the smallest value of k = 25 we infer 81 edges in the LGT network, connecting 17 nodes corresponding to orders of Proteobacteria (12), Low G+C Firmicutes (four) and High G+C Firmicutes (one). The largest clique inferred contains seven orders (Supplementary Figure 16 and Supplementary Table 23). At k = 30 only four of these orders remain (Neisseriales, Enterobacteriales, Pasteurellales and Lactobacillales) in the maximum clique, and this clique persists through to k = 40 (Supplementary Figures 17–19 and Supplementary Tables 25, 27, and 29). Thus these four orders form the core nodes of the GEC for the BAC dataset. This is in complete agreement with our results from the BA dataset.

At k = 25, we also infer Enterobacteriales (represented here by six E. coli genomes) to have donated via LGT to 125 genes in other orders and to have accepted LGT from other orders into 111 genes, together 34% of all affected genes across this dataset. These results support the developing themes of LGT being more successful among more-closely related genomes, with enteric bacteria and Firmicutes particularly active.

Other maximal cliques containing more than three orders are also found in the BAC dataset (Supplementary Tables 24, 26, 28, and 30) at different k. Some are subsets of the maximum clique, and reflect the fractions of the maximum clique in specific parts of the LGT network, while others are independent of the maximum clique and indicate other regions of dense connection among bacterial orders.

In addition to clique membership and topology, we are also interested in the biological processes enriched among the genes affected by inferred lateral events, as these may point to physiological, ecological and other processes that help to construct and maintain bacterial communities in nature. Analyses of the LGT networks inferred for the ECS and EB datasets reveal that more than 90% of the genes affected by LGT are represented in the corresponding cliques. In the ECS dataset, all vertices are in the maximum clique. In such cases there is no need to carry out enrichment tests: biological processes contributing to clique formation will be indistinguishable from those of the whole LGT network to which the cliques belong, i.e., the total LGT edge sets (Cong et al., 2016b). For the BAC dataset (Supplementary Table 31), genes annotated for involvement in metabolic processes (e.g., small-molecule and amino acid biosynthesis) are about twice as numerous as those in the next most-numerous category, ribosomal proteins: see the Supplementary Material for Cong et al. (2016b), particularly Section 4.2.

For the BA dataset, however, clique topologies change significantly with k. Few LGT events are detected at k > 30, particularly when sequences are grouped by phylum (Figure 4). For optimal comparison, we carried out enrichment tests on lateral genes of maximum cliques in each network at k = 20, 25 or 30, with genomes grouped either by phylum or by class. These tests identify biological processes related to metabolism, transport and regulation as over-represented when sequences are grouped by phylum. The term translational elongation (GO:0006414) ranks in first position at k = 20, and seventh at k = 25, among over-represented terms. The most significantly under-represented biological processes relate to transposition and to RNA modification at k = 20, and to RNA processing and biosynthetic processes at k = 25. The only term under-represented at k = 30 describes the modification of macromolecules.

When the genomes are grouped instead by class, the main categories of GO terms significantly over-represented remain those describing metabolism, transport and regulation. Those most under-represent relate to transposition, RNA metabolisms and regulation at k = 20 and 25; at k = 30, processes of protein modification are under-represented.

In general, the patterns of over-representation are similar between analyses at phylum and class levels. Interestingly, translation elongation is significantly over-represented at phylum level, but much less so at class level. Transposition (GO:0032196) is significantly under-represented in most cases.