LINflow: a computational pipeline that combines an alignment-free with an alignment-based method to accelerate generation of similarity matrices for prokaryotic genomes

Overview

In short, to minimize the number of computationally expensive ANI computations when constructing a genomic similarity matrix, LINflow sequentially adds genomes to a dataset, at each step efficiently identifies the genome already in the dataset that is most similar to the newly added genome, precisely calculates the ANI value only between this pair of genomes, and assigns a LIN to the new genome based on this ANI value and the LIN of the most similar genome. The LINs are then used to further accelerate the identification of the most similar genome and, most importantly, to efficiently infer all remaining pairwise ANI values to construct the complete genome similarity matrix. In other words, the main purpose of LINs in LINflow is to reduce the number of necessary ANI computations to one per genome when computing a genome similarity matrix.

The LINflow flowchart is shown in Fig. 1. When a new genome is added to the dataset, LINflow identifies the genome that is most similar to this new genome among the genomes that were previously added using the computationally efficient alignment-free tool sourmash (Brown & Irber, 2016; Pierce et al., 2019). This is accomplished via a two-step procedure, which consists of first identifying the LINgroup to which the new genome belongs and then identifying the most similar genome in this LINgroup. By default, LINflow uses the 95%-level LINgroup in this step, but this can be modified by the user based on the expected genomic similarities in a specific dataset. The precise ANI value between the two genomes is then computed using the more computationally expensive, but precise, pyani tool (Pritchard, 2014). LINflow uses this ANI value to assign a LIN to the new genome based on the LIN of its most similar genome (which LIN was previously assigned based on the ANI value to its most similar genome when that genome itself was previously added to the dataset). The assigned LINs can then be used to infer all-against-all ANI values even though only a single ANI computation is performed for each genome.

To make the results reusable and easily accessible in terms of reading and writing, a relational database managed by MySQL (https://dev.mysql.com/) is used to store data, with the schema shown in Fig. 2. This relational database connects tables with primary and foreign keys, and the connections between tables are represented by arrows. The genome table stores the locations of the genome sequences. The taxonomy table stores the taxonomic information corresponding to each genome in the database. LIN schemes (i.e., the number of LIN positions and the corresponding ANI thresholds), based on which LINs are assigned, are kept in the Scheme table. Besides three default LIN schemes, new schemes can be added by users so that LINs can be assigned according to the users’ needs in resolution. LINs are assigned with the three default schemes and with one user-defined scheme if there is any.

The individual steps of the pipeline are described in detail below.

Generation and storing of signatures

LINflow uses sourmash version 2.0 (Pierce et al., 2019) to generate signatures for both k = 21 and k = 51 with n = 2,000 (i.e., each signature consists of 2,000 hashes) for all genomes and stores them as individual files. The first member of each 95%-level LINgroup is chosen as the representative genome and a copy of its signature file is saved in a second directory together with signature files of all other representative genomes.

Choice of LIN scheme

LINflow allows the user to choose from four default LIN schemes: (1) a 20-position LIN scheme that ranges from 70% ANI to 99.999% ANI to cope with genus to strain level differentiation and is currently used in LINbase (Tian et al., 2020) (Table 1), (2) a 300-position scheme with positions starting at 70% and increasing by 0.1, and reaching 99.9% (used in this manuscript with the datasets listed below), (3) a 3,000-position scheme starting at 70% and reaching 99.99% at 0.01% intervals (recommended for constructing highly resolved similarity matrices for strains belonging to the same species), (4) a 300,000-position scheme starting at 70% and reaching 99.99999% at 0.00001% intervals. The user can also define any custom LIN scheme with any number of positions using any ANI value of choice (however, it is not advised to use ANI values below 70% since ANI does not reflect evolutionary relationships when ANI falls below 70%).

Initiation of LINflow

LINflow, by default, includes the 20-position LIN scheme as one of the schemes during each run. An arbitrary genome will be selected to assign the first LIN with 0 in each position. Its signature will be generated and saved as the representative of the LINgroup 0_A0_B0_C0_D0_E0_F using the 20-position default scheme both in the directories for representative genomes and this LINgroup.

LIN assignment

The new genome’s (G_Query) signature S_Query will be first queried against the representative genomes of all existing 95%-level LINgroups (or F LINgroups) with k = 21, using sourmash. Based on the analysis of more than 6,000 bacterial genomes from different families with k = 21, the Jaccard similarity of 0.2475 was found to correspond to 95% ANI (data not shown). If the highest Jaccard similarity J to one of the representative genomes S_Query is above 0.2475, then the corresponding genome represents the 95%-level LINgroup (L_95%), belongs to. If 0.0025 < J <= 0.2475, then the corresponding genome is at least 70% similar to G_Query, which means they share at least the A position in the default LIN scheme. If J <= 0.0025, the corresponding genome is less than 70% similar to G_Query.

If $J>0.2475$, S_Query will be queried against all members of L_95% by sourmash with k = 51. The most similar genome G_Subject according to Jaccard similarity in L_95% is identified as the most similar genome to G_Query in the entire database.

If $0.0025<J\le 0.2475$, S_Query will be queried against all the members of L_95% by sourmash with k = 21. The most similar genome G_Subject according to Jaccard similarity in L_95% is identified as the most similar genome to G_Query in the entire database.

For the above cases, ANI between G_Query, and G_Subject, ANI_Query, will be calculated with pyani. To assign the LIN to the query genome, G_Subject’s LIN, LIN_Subject, will be used as the reference from A to the last position that the ANI threshold is lower than or equal to ANI_Query, the first position that ANI threshold is larger than ANI_Query will be assigned a number that has not been used with the prefix in the database, and the rest of the positions will be filled with zeros. For example, ANI_Query = 95.4575%, it is over 95% at the F position but lower than 96% at the G position, so it will use LIN_Subject from A to F as LIN_Query’s A to F. At LIN_Query’s G position, a number that has never been used together with LIN_Subject’s prefix from A to F will be assigned. Each of LIN_Query’s H to T positions will be assigned 0.

If $J\le 0.0025$, no genome in the current database has over 70% ANI to G_Query, so that a new number that has never been used in the A position before will be assigned to LIN_Query’s A position, and the rest of LIN_Query will be filled with zeros.

Update of database and signature file system

G_Query, G_Subject, ANI_Query and LIN_Query will all be written to the database. If LIN_Query creates a new 95%-level LINgroup, a new directory for this LINgroup will be created and S_Query will be saved in this directory as a member and as a representative genome with other representative genomes, otherwise, S_Query will be only saved in the existing 95%-level LINgroup it belongs to.

Datasets

We compared the performance of LINflow with sourmash, pyani, and FastANI for four datasets of 484 genomes in total. Dataset A includes 248 genome sequences belonging to the genus Pseudomonas. Among the 248 genomes, 222 are from 46 named species. The remaining 26 genomes are from yet unnamed and undescribed Pseudomonas species. Dataset B consists of 43 Xanthomonas perforans genomes. Dataset C includes 140 genomes of the genus Lysinibacillus, whereby 96 of them belong to 27 named species. Dataset D includes 53 genomes from the species Xylella fastidiosa and two genomes from the species Xylella taiwanensis. A separate dataset E with the genome sequence of Pseudomonas caeni DSM 24390 was used to assess the computational speed of the above tools when adding a new genome to the already-analyzed dataset A. Table S1 lists the genomes included in each dataset and the actual genome sequences can be accessed directly in this repository: https://code.vt.edu/linbaseproject/linflow_datasets.

Comparison of tools

LINflow was compared with pyani (blast option), sourmash and FastANI in regard to speed, memory usage and accuracy. Parameters used when running these programs are listed in Table 2. Note that for each pair of genomes A and B, pyani computes ANI by calculating both the ANI of A to B and the ANI of B to A. We used the average of the two pairwise ANI values. The calculations were executed on a 2.4 GHz CPU on Cascades, an Advanced Research Computing system at Virginia Tech and the execution time and memory usage were monitored by the job scheduler built in the system.

Hierarchical clustering using the complete linkage method was applied to the similarity matrices of each dataset calculated by the four software suites using custom R scripts. Heatmaps were generated based on the hierarchical clustering results to investigate whether FastANI, LINflow and sourmash are able to classify bacteria into species as accurate as pyani. Rows and columns of the heatmaps were reordered to be in the same order as the heatmap generated by pyani.

Finally, pairwise Mantel tests (Mantel, 1967) in combination with Pearson correlation coefficients were computed using custom R scripts to determine how well the similarity matrices (ANI for LINflow, FastANI and pyani and Jaccard similarity for sourmash) obtained with the different tools correlated with each other.

Computational speed and memory usage

The CPU time needed to analyze datasets A, B, C and D by each software is shown in Table 3. For sourmash, two separate times are indicated since sourmash runs two separate commands: “sourmash compute”, which computes signatures and “sourmash compare”, which computes the pairwise Jaccard similarity between signatures. Although the FastANI workflow is also split into two phases, the indexing phase and the compute phase, the two phases cannot be executed separately. Therefore, only the total execution time is shown.

Sourmash was the fastest out of the four tools, being 1,000–4,000 times faster than pyani, depending on the dataset. FastANI and LINflow took a similar time to execute compared to pyani, being 84–253 times and 20–150 times faster, respectively.

The memory usage of LINflow was lower than pyani and lower than FastANI when analyzing datasets A and C, which have relatively larger numbers of genomes compared to datasets B and D (Table 4). Furthermore, the time cost for adding a new genome to an existing dataset by LINflow did not increase as significantly as for pyani and FastANI. This can be seen by comparing the time cost for adding a new genome to dataset A (Table 5) with the average processing time for each genome in dataset A (Table 3).

Accuracy

Similarity matrices obtained with sourmash, FastANI and LINflow were compared to those obtained with pyani, which we considered the gold standard, since it is exclusively based on BLAST reflecting the original description of ANI (Konstantinidis & Tiedje, 2005a, 2005b). For sourmash, we determined its performance separately for k = 21 and k = 51. For LINflow, we used two of the four default schemes: the 20-position scheme used in LINbase (Tian et al., 2020) (see Table 1 for the LIN scheme and Table S1 for the result) and the 300-position scheme (with ANI values increasing from 70% at the left-most LIN position to 99.9% at the right-most LIN position with 0.1% intervals between neighboring positions). The LINbase scheme was used to assign LINs to the genomes and classify them as LINgroups. The 300-position scheme was used to determine the ANI similarity matrix. After similarity matrices were computed with all tools, heatmaps were generated to visualize the genomic relatedness among the analyzed genomes.

Heatmaps derived from the ANI matrices obtained with pyani (Fig. 3A), FastANI (Fig. 3B) and LINflow (Fig. 3C) show the same species level (ANI ≥ 95%) clustering of the 247 Pseudomonas genomes of dataset A visible as red blocks along the diagonal. Five major clusters are easily visible. Cluster 1 consists of genomes belonging to P. aeruginosa, cluster 2 represents the species P. chlororaphis, clusters 3, 4 and 5 constitute the P. syringae species complex and related genomes. Note that LINflow not only classified the genomes as species but also distinguished intraspecific groups as LINgroups. The LIN prefixes denote LINgroups and show both intergroup and intragroup relationships.

Sourmash was able to perform species-level clustering with k values of both 21 (Fig. 3D) and 51 (Fig. 3E). The obtained results suggest that Jaccard similarity calculated with k = 51 only weakly correlates with ANI for low ANI values compared to k = 21, for example, clusters 2, 4 and 5. Instead of genomes that are highly similar to each other, for example, genomes in cluster 5, k = 51 provides higher resolution than k = 21.

To further evaluate the accuracy of LINflow compared to the other tools, the complete similarity matrices obtained for all four datasets with all tools were compared with each other using the Mantel test (Mantel, 1967) using Pearson’s correlation coefficients. Results are reported in Fig. 4. One can easily see how the results obtained with LINflow are highly correlated with those obtained with pyani for datasets A, C and D (Pearson correlation coefficient of 0.99 or 1) but not for dataset B, which has a Pearson correlation coefficient of only 0.78. Since many pairs of genomes in set B have ANI values above 99.8% and differ from each other by less than 0.1%, we computed ANI values for set B also using the 3,000-position and a 300,000-position LIN scheme hypothesizing that the lower correlation was due to rounding of ANI values to the first decimal place when using the 300-position LIN scheme. However, switching to the 3,000-position and 300,000-position LIN scheme only increased the correlation between LINflow and pyani slightly to a Pearson correlation coefficient of 0.82.

When comparing Pearson Correlation coefficients for the FastANI vs. pyani comparison and the sourmash vs. pyani comparison with the LINflow vs. pyani comparison, LINflow shows the same or higher correlation with pyani for all datasets with the exception of dataset B, for which LINflow shows the lowest correlation with pyani.

Finally, to determine how LINflow and FastANI correlate with pyani from low to high pairwise ANI values, we plotted all pairwise ANI results obtained with LINflow and FastANI against pyani results similar to what was done by Jain et al. (2018c). Figure 5 shows how FastANI and LINflow both correlate very well with pyani for ANI values above 85% but deviate from pyani at ANI values from 85% to 70%. While ANI values computed with FastANI are higher than the corresponding pyani ANI values in this range, ANI values inferred by LINflow merge into a relatively small number of ANI values for many different pyani ANI values (see “Discussion” for an explanation of this phenomenon).