Fast estimation of genetic relatedness between members of heterogeneous populations of closely related genomic variants

Background Many biological analysis tasks require extraction of families of genetically similar sequences from large datasets produced by Next-generation Sequencing (NGS). Such tasks include detection of viral transmissions by analysis of all genetically close pairs of sequences from viral datasets sampled from infected individuals or studying of evolution of viruses or immune repertoires by analysis of network of intra-host viral variants or antibody clonotypes formed by genetically close sequences. The most obvious naïeve algorithms to extract such sequence families are impractical in light of the massive size of modern NGS datasets. Results In this paper, we present fast and scalable k-mer-based framework to perform such sequence similarity queries efficiently, which specifically targets data produced by deep sequencing of heterogeneous populations such as viruses. It shows better filtering quality and time performance when comparing to other tools. The tool is freely available for download at https://github.com/vyacheslav-tsivina/signature-sj Conclusion The proposed tool allows for efficient detection of genetic relatedness between genomic samples produced by deep sequencing of heterogeneous populations. It should be especially useful for analysis of relatedness of genomes of viruses with unevenly distributed variable genomic regions, such as HIV and HCV. For the future we envision, that besides applications in molecular epidemiology the tool can also be adapted to immunosequencing and metagenomics data.

sampled from infected individuals [1,2]. Viral populations, for which the minimal inter-sample distance does not exceed the threshold, are considered to be potentially linked by transmission [1], while the number of pairs in P may suggest the time since a transmission event [3]. The related genetic network construction problem aims to build a graph with vertices corresponding to sequences from a given dataset T and edges corresponding to all pairs of sequences with an edit or Hamming distance less than the threshold t. This problem arises in studying and analysis of viral populations [4] or antibody repertoires [5]. Similar problems also emerged under different names in various areas of computer science [6][7][8][9][10].
The edit distance between a pair of sequences can be calculated in time O(L 2 ) using dynamic programming [11]. If only distances below a desired threshold t which is small relative to L are desired. The distance calculation can be carried out with a small subset of diagonals neighboring the main diagonal of the dynamic programming matrix, leading to O(tL) time algorithm [12]. In this case a naïve algorithm for the similarity join problem requiring pairwise comparison of all sequences has an asymptotic running time O(tLN 2 ), which is still impractical for more than several thousand sequences.
Several filtering-based approaches have been put forward to improve the efficiency of the similarity join-type problems by reducing the number of pairs to be compared. Note that while fast heuristic and approximate methods exist such as Shingling [13], LSH [7], or BLAST [14], this paper focuses on the problem of exact distance calculation.
The common filtering approach is based on on the fundamental idea that related sequences should share long k-mers (substrings of length k) [15]. Several existing methods rely on signature schemes to quickly locate feasibly linked pairs [6] by assuming that pairs with an edit or Hamming distance which does not exceed a threshold t will share at least a certain number of k-mer-based signature keys. However, straightforward application of this technique to viral sequencing data is not sufficiently efficient, since mutations are not distributed uniformly along viral genomes, but tend to concentrate in short hypervariable regions [16]. As a result, many viral sequences share k-mers, thus significantly reducing the efficiency of filtering. The same effect has been observed for immunosequencing data [5], where all antibodies originating from the same V gene often share a k-mer from that gene.
In this paper, we describe a tool which uses k-mer-based signature filtering scheme optimized for viral data to solve the following problems: • Sample pair filtering: given two NGS sequence samples T 1 and T 2 , quickly determine whether the distances between all inter-sample pairs of sequences are greater than the threshold t. • Inter-sample sequence retrieval (similarity join): given two NGS sequence samples T 1 and T 2 , find all inter-sample pairs of sequences at edit distance or hamming distance below the threshold t. • Intra-sample sequence retrieval (or genetic network construction): given an NGS sequence sample T 1 , find all pairs of sequences at edit distance or hamming distance below the threshold t.
The tool was validated using Hepatitis C Virus (HCV) data in the settings used for detection of viral transmissions and outbreaks [1,2].

Notation
In the methods description, we assume that input sequence samples T 1 and T 2 both contain N sequences of length L, which cover the same genomic region. From here onwards we will use k as a fixed predefined parameter.
Further we will use the following notation: . . s j -subsequence of S starting at position i and ending at position j. • k-mer -any subsequence of length k • k-segmentk -mer that starts at a position 1 + ik, i = 0, 1, 2, . . .. • K(S) -the set of all k -mers of the sequence S.
• R(S) -the family of all k -segments of the sequence S (possibly with repetitions). • h(S, Q) -Hamming distance between two sequences S and Q • l(S, Q) -edit distance (Levenshtein distance) between two sequences S and Q distance, as mentioned above, could be calculated using dynamic programming [12] Main Data Structure Our signature-based filtering scheme is based on the following simple observation: Proof If S and Q differ by an edit distance of t, then by the pigeon hole principal at most t k-segments differ between the sequences S and Q. So at least m − t k-segments must be the same.
Thus we need a fast way to calculate the number of common k-segments and k-mers for a given pair of sequences. To do it we introduce a hash function: where f : {A, C, G, T} → {0, 1, 2, 3} is an arbitrary bijection. For k-mers with k < 32, this hash function allows us to store them as 64-bit integers and can be quickly recursively calculated as follows: In addition, the hash can be inverted and so only the hash values of k − mers need to be stored. In the proposed framework, each sample T is stored using a data structure further referred to as a T-dictionary and denoted by dict(T), which consists of the following fields: • dict(T).HM -an inverted index of T [17], i.e. a hash table, where each key is a k -mer hash and its value is a set of all sequences from T that contain this k -mer.

Algorithm Description
We will first describe the approach for the sample pair filtering problem. Building a simple and fast filter for unrelated samples T 1 and T 2 is easy by applying Proposition 1.1 to whole samples as follows. Recall that T 1 and T 2 are considered to be genetically related, if the minimal edit distance between their sequences does not exceed the threshold t. Given two dictionaries dict(T 1 ) and dict(T 2 ), the necessary condition for their genetic relatedness is an existence of at least The sample pair filter pseudocode is presented at Algorithm 1.

Algorithm 1 Simple filter for unrelated samples
for lSegmentHashes ∈ dict1.SC do 4: for hash ∈ lSegmentHashes do 5: if hash ∈ dict2.KM then 6: coincidences ← coincidences + return coincidences ≥ m − t 12: end function 13: Assuming that membership verification for a hash set dict(T).KM can be performed in time O(1), the worstcase running time of the filter is O(NL). In real settings, samples with genetically related sequences produce significantly smaller maps dict(T).SC, thus leading to a lower average running time than in the worst case.
The algorithms for inter-sample sequence retrieval and intra-sample sequence retrieval problems are very similar, so we will describe the approach for the former problem. As before, let T 1 and T 2 be two samples. The algorithm first constructs the set of candidate neighbors CN S ⊆ T 2 for every sequence S ∈ T 1 . This procedure (the filtering), is followed by the verification procedure, which calculates actual neighbors of all sequences S ∈ T 1 by calculating distances between S and all sequences S ∈ CN S . The pseudocode for inter-sample sequence retrieval algorithm is presented in Algorithm 2.
The basic filtering strategy utilizes Proposition 1.1, with the following features aiming at improvement of the running time. For each S ∈ T 1 , the set CN S can be implemented as a hash table, with keys being sequences S ∈ T 2 and values CN S (S ) being numbers of matches between k-segments in S and k-mers of S . Let L S be the number of k-segments of S that occur as k-mers in T 2 , and I = (i 1 , i 2 , . . . .i L S ) be the list of starting positions of these k-segments. To calculate the number of matches between k-segments in S and k-mers of S we may iterate over the list I and increment the current value of CN S (S ), when necessary. If after j iterations the inequality does not hold, then S cannot accumulate the required number of matches with the remaining iterations, and therefore the sequence S can be filtered out right away. These considerations imply that the order in which starting positions of k-segments are examined is important in determining the algorithm's running time.
The order of k-segment starting positions is determined heuristically as follows. For each position i let k S (i) = |dict(T 2 ).HM(S[ i : i+k−1] )| be the number of sequences from T 2 that contain the i-th k-segment from S. If we sort positions by ascending order of the numbers k S (i) it usually leads to faster pruning of sequence pairs as this order minimizes the size of the candidate set that must be examined at each iteration.
Another simple adjustment could be implemented using the fact that the hamming distance is an upper bound for an edit distance, while the calculation of the former is significantly faster. Therefore if h(S, Q) ≤ t, then Q can be added to the list of neighbors of S without the edit distance calculation.

Hamming distance adjustment
The filtering strategy described above can be further improved, if the input sequences are aligned to a reference. In this case the samples can be compared using Hamming distance instead of an edit distance. For Hamming Algorithm 2 Signature-based filter to find sequence pairs closer than threshold 1 order ← indexes of segments sorted in ascending order by their frequencies in T 2

5:
L s ← number of k-segments from s that appear in T2 6: for s' ∈ keys of CN s do 15: if h(s, s') ≤ t or led(s, s') = −1 then 16: s and s' are related 17: end if 18: end for 19: end for 20: end function 21: 22: function FILL adds all sequences from T2 that share the same k-mer with s 23: function FILL (CN s , s, dictT1, dictT2,  Furthermore, genomic heterogeneity is distributed highly irregularly along the genomes of species of interest. For example, Fig. 1 illustrates the distribution of nucleotide entropy for a particular intra-host population along the 264bp-long genomic HCV region at the junction of envelope glycoproteins E1 and E2, which is often used in epidemiological and immunological studies [1,18,19]. It should be noticed that k-segments from conserved regions are significantly less useful for Fig. 1 Distribution of nucleotide entropy along the E1/E2 region of HCV for a population of 469 unrelated genotype 1a sequences obtained from NCBI the filtering as we want to maximize detectable differences between tested sequences. The non-uniformity in genomic heterogeneity can be taken into account by switching to the framework with k-segments of unequal size. By selecting k-segment boundaries that contain roughly equal amounts of average information entropy over the dataset, the filtering speed and quality could be significantly improved. Figure 2 provides an example, when entropy-based k-segments length allows more accurate filtering than uniform k-segments length. Formally, let H i = − 4 j=1 P(x i j )log 2 (P(x i j )) be the sample nucleotide entropy at position i, where P(x i j ) is a frequency of nucleotide x i j on i-th position of the alignment. The segments are selected in such a way that for every

Validation Data
The developed tool was validated using NGS datasets of intra-host HCV populations sampled from infected individuals. Each dataset contains the E1/E2 junction of the HCV genome of length 264nt, which contains the Hyper Variable Region 1 (HVR1) region. Each sample was processed by error correction and haplotyping tools, and as a result we receive as an input datasets consisting of unique HCV haplotypes.
We used a set of 413 samples from [2] with 501.5 haplotypes per sample in average produced by NGS; 8 datatsets d 1 , . . . , d 8 with 1000, 2000, . . . , 128 000 sequences produced by random sampling from NGS dataset with sequences sampled from chronically infected individuals and one additional NGS dataset m 1 consisting of 10 467 sequences. The data are available in tool's repository.
In all tests, the threshold t = 3.77% ≡ 10nt was used. This value is derived in [1] as empirically validated recommended threshold for separation between epidemiologically related and unrelated intra-host HCV HVR1 populations.
All tests were run on server with 128 Intel Xenon E7-4850 2.1GHz cores and 1.4Tb RAM. For Inter-sample sequence retrieval desktop PC was used with 4 Intel(R) Core(TM) i7-5500 2.4GHz cores and 8Gb RAM. All code is written on Java to provide a threaded, platform independent solution.

Sample pair filtering and Inter-Sample Sequence Retrieval validation
For Sample pair filtering and Inter-Sample Sequence Retrieval problems, we validated the tool using HCV datasets from [2]. The proposed approach has been compared with the Filter Composition pipeline proposed in [2]. Both methods were run on a desktop computer, as in the original paper [2]. The results are reported in Table 1.
Here we show the result of comparison of all pairs of samples and all inter-sample pairs of sequences. The proposed sample pair filtration algorithm removed 92% of all possible samples pair comparisons, and sequence pair filtering algorithm managed to filter out 99.996% of all possible sequence pairs. The latter means that only 888,914 out of 22,037,502,011 sequence pairs passed from filtering to verification stage of the algorithm. As a result, the proposed approach significantly outperforms the Filter Composition Pipeline in filtering quality and in running time.
We studied how the filtering quality is affected by different optimization subroutines (Table 2). Disabling sample pair filtering increases the running time for comparison of all samples by 42%, while the impact of sorting of ksegment starting positions is even higher, with disabling of this step slowing down the comparison by 254%.
Preprocessing and dictionary building can take up a significant portion of the total running time of a signaturebased filtering algorithm, when samples are distant and few distance calculations are required. For the given collection of 413 samples, preprocessing of all samples takes ∼ 4840ms, which constitutes approximately 1/3 of the total running time of the algorithm. Note that in the case when significant number of closely related sequence pairs is present, the situation is different (see the next section).  The algorithm performance depends on the size of the k-mers and k-segments. Small k leads to larger number of matches between k-segments and k-mers of distant sequences, which can cause extra sequences to be added to the candidate lists thus leading to decrease in filtering quality. Larger k leads to fewer false matches but unfortunately also a larger k-mer dictionaries. We examined different k-mer sizes to determine the optimal size for our datasets and found that k = 11 gives the best performance.

Intra-sample Sequence Retrieval Validation
For Intra-sample Sequence Retrieval Problem, we validated the proposed approach using datasets d1,...,d8,m1. First, for it was compared with a single-thread, brute force method with the worst-case complexity O(N 2 Lt), which performs pairwise comparison of all sequences and calculates limited edit distance using dynamic programming as described in [12]. The results are presented in Table 3.
The running time of the proposed tool was also compared with the running time of a recently published method from [5], which was originally designed for the analogous problem for immunosequencing data. Figure 3 illustrates that signature-based filtering approach demonstrates the significant advantage. Figure 4 demonstrates that for aligned sequences in most cases the adjustment utilizing entropy-based variable-size k-segments allows to achieve a significant speedup with respect to a constant-size k-segment.
The speedup described above is achieved by the combination of the several features. The first feature is the quality of filtering, which is analyzed in Tables 4 and 5. On average, only ∼ 10% of sequence pairs that pass filtering step ("false positives") are not genetically related. As expected, most of the false positive pairs were very close to the threshold (Fig. 5). With the threshold set at t = 10, pairs with an edit distance of l(S, Q) = 11, 12, 13 represent up to 75% of all false positives. Pairs that are so close to the threshold are difficult to filter out.   Another important feature is the fact that as the input increases in size the runtime of the algorithm is dominated by the edit distance calculations (Fig. 6). However, the filtering and the Hamming distance shortcut reduces the number of edit distance calculations that must be performed. As a result, the actual edit distance is only calculated on small portion of the total pairs from the dataset (Table 4).
We attempted to improve the filtering performance using other methods such as k-mer similarity [20], true matches [6], Hamming radius filter [2]. However, the overhead of these methods was greater than any runtime savings.

Discussion
In this paper we presented an efficient signature-based tool to solve problems of edit or Hamming distance sequence retrieval for NGS data obtained from heterogeneous viral populations. It outperforms other approaches to this problem by including several data-specific steps and filters. The proposed approach was designed having problems of computational molecular epidemiology in mind. Until recent years, genomic analyses of viral transmissions and epidemic spread used a single viral sequence per infected individual. The advent of sequencing technologies now allows to analyze thousands of viral haplotypes per patient. Furthermore, just in the United States, from 2.7 million to 3.9 million people are infected with HCV [21], while ∼ 1.1 million people are infected with HIV [22]. These numbers put an immense computational burden on real-time advanced molecular surveillance systems, such as Global Hepatitis Outbreak Surveillance Technology (GHOST) [23], which is currently being deployed by Centers for Disease Control and Prevention. When deployed, such system should have computational capacity to identify, whether a query set of viral samples is genetically related with any sample from a database consisting of hundreds of thousands of samples each consisting of thousands of sequences. The proposed approach aim to allow to process such queries efficiently. It builds on the general idea proposed in [6], which is heavily optimized by utilization of efficient data structures, such as inverted indexes and hash maps, and introduction of running time-improving procedures, such as efficient hash values calculation and determination of optimal order of k-mers processing. The proposed optimization steps allows for more than 2.5-fold running time decrease in comparison with the non-optimized filtering (Table 2). Furthermore, the proposed method takes into account uneven distribution of heterogeneous position along viral genomes by using variable entropybased k-mers. It allows to improve both filtering quality (Fig. 2) and speed (Fig. 4). In general, for viral samples comparison the proposed filtering approach allows to eliminate the overwhelming majority of sequence comparisons and achieve a substantial running time decrease (Tables 1, 2, 3, 4 and 5).

Conclusion
The proposed tool allows for efficient detection of genetic relatedness between genomic samples produced by deep sequencing of heterogeneous populations. The tool is freely available for download at https://github.com/ vyacheslav-tsivina/signature-sj. It should be especially useful for analysis of relatedness of genomes of viruses with unevenly distributed variable genomic regions, such as HIV and HCV. For the future we envision, that besides applications in molecular epidemiology the tool can also be adapted to immunosequencing and metagenomics data.

Funding
Publication of this manuscript was sponsored by IH grant 1R01EB025022-01

Availability of data and materials
The datasets generated and analyzed during the current study are available in the github repository https://github.com/vyacheslav-tsivina/signature-sj/tree/ master/data

About this supplement
This article has been published as part of BMC Bioinformatics Volume 19 Supplement 11, 2018: Proceedings from the 6th Workshop on Computational Advances in Molecular Epidemiology (CAME 2017). The full contents of the supplement are available online at https://bmcbioinformatics.biomedcentral. com/articles/supplements/volume-19-supplement-11.