PanDelos: a dictionary-based method for pan-genome content discovery

A pan-genome can be abstractly considered as a structure defined on a set of genomes. The structure is built by identifying groups of homologous genes [1]. Two genes are homologous if they share a common ancestral gene. Homologous genes can be distinguished into paralogous, when homology occurs within the same genome, or orthologous, when homology occurs between different genomes. We call pan-genome content discovery the determination of homologous groups within a collection of genomes.

Different mechanisms are involved in gene transmission. Paralogy is linked to sequence duplication within the same genome. Orthology is associated to a “vertical” transmission. It happens among genomes in the same lineage and involves most of the genetic contents. On the contrary, “horizontal” transmission occurs between genomes of organisms of different lineages, involving one or few genes. Genes present in every genome are core genes of the pan-genome and they may be involved in essential living functionalities. Sequences shared by a subset of genomes are referred as dispensable and they represent variable features. Singleton genes are present only in one genome and represent some genome-specific functionality. The collective analyses of all the genes is developed for many specific interests, for example, for the study of a bacterial strain of a given species [2, 3]. Pan-genome analyses found many application in clinical studies [4, 5], for example they help in identifying drug-target genes in clinical studies [6, 7], or in exploring phylogenetic lineages of bacteria [8] that can be linked to strain-specific disease phenotypes [9].

Approaches to pan-genome content discovery need to take into account that gene duplication and transmission may introduce sequence alterations [10–13]. The variations make the task of recognizing homologous genes difficult, especially when ancestor genomes are no more available. Core genes are often under strong evolutionary selection, thus their sequences are transmitted almost without any alteration. The amount of variations affecting dispensable genes varies and the similarity between homologous sequences tends to decrease according to their phylogenetic distance. When closely related organisms are analyzed, reasonable thresholds on sequences similarity are applied to recognize gene families. However, when distant genomes are compared, global thresholds result less feasible. Suitable notions are needed to define adaptive thresholds especially when they present non-uniform phylogenetic distances.

The discovery of a pan-genome content is an NP-hard problem [14], and the complexity of the analysis is proportional to the number of input genomes. This is mainly due to the fact that all-against-all comparisons between gene sets are required to solve the task. State of art tools for pan-genome analysis are Roary [15] and EDGAR [16]. They use some heuristics to reduce the computational requirements in the definitions of thresholds for sequence alignments and the number of comparisons necessary in their procedures. Both approaches are based on a largely-used strategy that searches for reciprocally most similar genes between compared genomes [17, 18].

Roary combines an approach for clustering gene sequences (CD-HIT) with a procedure based on reciprocal BLAST alignments. CD-HIT [19] clustering counts the presence of k-mers, substrings of length k, among the analyzed sequences at different values of k. The results of CD-HIT are merged with normalized BLAST scores [20] and clustered via the MCL algorithm [21]. The Roary’s procedure requires intensively tuning of user-defined parameters to set the thresholds for discarding low homology values. Parameters are set globally, making Roary best performing on closely related genomes.

EDGAR uses adaptive thresholds depending on the distribution of BLAST gene scores. The retrieval of a distribution is made feasible by the normalization of alignment scores. The normalization is performed by fixing the self-alignment score of a sequence as the maximum one. The approach results suitable for comparing genomes with a considerable phylogenetic distance, but some disadvantages arise. It requires an expensive amount of sequence alignments, in fact for each 1-vs-1 genome comparisons, the complete gene sets of the two genomes must be cross-aligned. EDGAR chooses the threshold on the minimum feasible score by computing the distribution of the normalized gene blast of all scores. Scores are summed up and represented in a histogram, and a beta distribution is calculated from the mean and standard deviation of the observed values. A 97% quantile of the density function is used as a cutoff to asses orthology. The quantile has been identified by manual inspection of hundreds of histograms from real cases.

Roary and EDGAR are based on sequence alignment, however alternative strategies can be used for retrieving domain architecture between homologous genes [22] or for the detection of horizontal gene transfer [23], by exploiting alignment-free techniques.

We present PanDelos ¹, a methodology to discover pan-genome content in phylogenetically distant organisms based on information theory and network analysis. It is parameter-free, the thresholds are automatically deduced from the context. PanDelos avoids sequence alignment by introducing a similarity measure based on k-mers multiplicity, rather than simple presence/absence of mers. The strength of the strategy is supported by a non-empirical choice of the best appropriate k-mer length. Moreover, the selection of minimum similarity for which two sequences are eligible to be homologs is inspired by the knowledge coming from read mapping in next-generation sequencing and sequence reconstruction processes. Reciprocal best hits in 1-vs-1 genome comparison, aimed at discovering orthologous genes, are used as a basis to infer thresholds for paralogs discovery. Homology relations are integrated into a global network and groups of homologous genes are extracted from it by applying a community detection algorithm. PanDelos outperforms in terms of running times and quality discovery contents the existing approaches, Roary and EDGAR, in real applications and in synthetic benchmarks, that represent fully trusted golden truth.

The detection of gene homology performed by PanDelos is divided into 5 main steps that combine a candidate selection based on k-dictionaries, with a refinement procedure, developed by means of network analysis. Firstly, an optimal value of word length k is chosen according to properties of the input collection of genomes. Consequently, genes are compared and candidate homologous pairs are selected. The selection is firstly applied by setting a minimum amount of intersection between the k-dictionaries of two genes. Then, the generalized Jaccard similarity is used to measure the similarity between genes in order to extract bidirectional best hits. The extraction produces a homology network from which, at the end of a refinement procedure, gene families are retrieved. Figure 1 gives an overview of the overall schema.

In what follows, we first describe the details of PanDelos together with the engineering and extension of existing data structures that allows PanDelos to reach high performance and efficiency.

Data structures engineering for fast similarity computation

Limitations of enhanced suffix arrays for pan-genome computations

Given a string s, a suffix array (SA) [27] reports the lexicographically ordered suffixes of s equipped with their start position in s. Substring search by means of SA can be sped up by performing binary searches. An enhanced suffix array (ESA) [28] is a combination of SA with the LCP (Longest Common Prefix) array giving the length of the longest common prefix of a suffix with that one lexicographically preceding it. An ESA allows for efficient recovery of the k-mers multiplicities [29]. The values of the LCP array define contiguous regions of the ESA array, called LCP-intervals, which identify all the occurrences of k-mers. Additionally, an array of length N reports for each suffix the distance from its start to the first forward occurrence of a N symbol [30]. The N is used to represent positions in s that must be discarded in dictionary operation. The ESA structure performs k-mer enumeration in linear time by just doubling the memory requirement of simple SA. Since each k-mer must be checked for N inclusion, this verification increases the time complexity by a factor of k. However, with the additional N array, the complexity remains linear. Figure 2a gives an example of ESA+N structure that has been built for the string WLLPPP, and illustrates LCP intervals of 1-mers and 2-mers of the string.

Methodologies for efficient similarity calculation in sets of strings have been developed by means of suffix trees [31]. This type of data structure inspired the develop of suffix arrays as a memory efficient implementation that does not increase time requirements. However, enhanced suffix arrays are useful for single reference string analysis, but result less suitable to be applied to sets of strings.

Given two sequences, s and t, the comparison of their k-dictionaries can be performed by listing the k-mers of s and searching for them in t. Taking into account an ESA+N structure, the k-mers listing is performed in |s| time, and the search of each k-mer in t takes \(\mathcal {O}(k \cdot log(|t|))\) time. Since at most |s| distinct k-mers are in s, the overall time is \(\mathcal {O}(|s| \cdot k \cdot log(|t|))\).

The search described above takes into account only D_k(s), but t may contain k-mers not listed in the dictionary of s, thus the time requirement is doubled because D_k(t) must also be scanned. In 1-vs-1 genome comparison, the process must be repeated for every pair of genes belonging to the two genomes, resulting in a highly expensive procedure. In the next section, we address a way to efficiently improve the procedure.

PanDelos data structure engineering

We extend the ESA+N structure [30] in order to speed up the comparison of k-dictionaries when multiple sequences are taken into account. The goal is to compute the generalized Jaccard similarities between a set of sequences simultaneously.

The generalized Jaccard similarity between s and t can be expressed as:

$$J_{k}(s,t) = \frac{a}{b + c}, $$

where a is the sum of the minimum multiplicities of k-mers shared by the two sequences, b is the sum of the maximum multiplicities of the shared k-mers, and c is the sum of the multiplicities of k-mers appearing only in one of the two sequences.

Given a and b for every pair of sequences, then c is obtained as

$$c = (|s| + |t| - 2k + 2) - (a + b), $$

where (|s|+|t|−2k+2) is the sum of multiplicity of all k-mers in s and t. Therefore it can be rewritten as:

$$J_{k}(s,t) = \frac{a}{|s| + |t| - 2k + 2 - a}. $$

Given two genomes, \(\mathbb {G}^{1}=\{s_{1}, \dots, s_{n}\}\) and \(\mathbb {G}^{2}=\{t_{1}, \dots, t_{m}\}\), genes are concatenate in a single global sequence s₁·N·⋯·N·s_n·N·t₁·N·⋯·N·t_m. An ESA+N structure is built, and the concatenation by N symbols ensures that extracted k-mers do not cross between gene sequences. The data structure is extended with a further array, called SID. Given an LCP-interval, that represents a specific k-mer, the content of the corresponding interval in the SID array reports the identifiers of the sequences in which the k-mer is present. Moreover, the number of time a sequence identifier is repeated within the interval corresponds to the multiplicity of the k-mer within the specific sequence.

For each pair of sequences involved in the interval, we compute the sum of minima (the a term in the generalized Jaccard formula) by computing the partial of such sums and storing them in a matrix M

$$M[i,j] = \sum\limits_{w \in D_{k}\left(s_{i} \cup t_{j}\right)} min\left(c_{s_{i}}(w), c_{t_{j}}(w)\right), $$

for 1≤i≤n and 1≤j≤m.

Once the matrix is filled, the similarities are finally computed by the formula:

$$J\left(s_{i},t_{j}\right) = \frac{ M[i,j] }{ \left|s_{i}\right| + \left|t_{j}\right| -2k +2 - M[i,j]}. $$

In this way, we avoid comparisons between sequences that do no share any k-mers, and eliminate the logarithmic factor of searching k-mers into multiple suffix arrays, one for each gene sequence.

Similarly, two additional matrices are stored in order to calculate candidate homologous genes: P₁[i,j] reporting the percentage of multiplicities of k-mers in s_ishared with t_j, and P₂[i,j] storing the vice versa:

$$P_{1}[i,j] = \sum\limits_{w \in D_{k}(s_{i} \cap t_{j}) }\frac{c_{s_{i}}(w)}{|s|-k+1} $$

$$P_{2}[i,j] = \sum\limits_{w \in D_{k}(s_{i} \cap t_{j})} \frac{c_{t_{j}}(w)}{|t|-k+1}. $$

Figure 2 reports an example for four sequences that are concatenated into a single one. Then 2-mers, and their associated sequences identifiers, are retrieved from the data structure. In the example, two genomes are compared. The first genome contains the sequences s₁ and s₂, and the second genome contains the sequences t₁ and t₂. The word length 2 is chosen as best k value, thus sequences are compared by means of the multiplicities of the 2-mers they contain. Ideally, the matrix (c) has to be computed in order to calculate the Jaccard similarity between the sequences. For higher values of k, the storage and update of such matrix may require high computational efforts, thus its rows are computed on the fly by identifying k-mer intervals along the indexing structure and by iterating them. A linear iteration over the structure lists the 2-mers, together with the number of times they appear within each original sequence. During the iteration the three matrices M, P₁ and P₂, in Fig. 2d, e and f, are updated. After that every 2-mer have been iterated, Jaccard similarity are computed by means of the M matrix, while coverage percentage are computed by means of the P₁ and P₂ matrices.