Finding all maximal perfect haplotype blocks in linear time.

Recent large-scale community sequencing efforts allow at an unprecedented level of detail the identification of genomic regions that show signatures of natural selection. Traditional methods for identifying such regions from individuals' haplotype data, however, require excessive computing times and therefore are not applicable to current datasets. In 2019, Cunha et al. (Advances in bioinformatics and computational biology: 11th Brazilian symposium on bioinformatics, BSB 2018, Niterói, Brazil, October 30 - November 1, 2018, Proceedings, 2018. 10.1007/978-3-030-01722-4_3) suggested the maximal perfect haplotype block as a very simple combinatorial pattern, forming the basis of a new method to perform rapid genome-wide selection scans. The algorithm they presented for identifying these blocks, however, had a worst-case running time quadratic in the genome length. It was posed as an open problem whether an optimal, linear-time algorithm exists. In this paper we give two algorithms that achieve this time bound, one conceptually very simple one using suffix trees and a second one using the positional Burrows-Wheeler Transform, that is very efficient also in practice.

Introduction and background

As a result of the technological advances that went hand in hand with the genomics efforts of the last decades, today it is possible to experimentally obtain and study the genomes of large numbers of individuals, or even multiple samples from an individual. For instance, the National Human Genome Research Institute and the European Bioinformatics Institute have collected more than 3500 genome-wide association study publications in their GWAS Catalog [1].

Probably the most prominent example of large-scale sequencing projects is the 1000 Genomes Project (now International Genome Sample Resource, IGSR), initiated with the goal of sequencing the genomes of more than one thousand human individuals to identify 95% of all genomic variants in the population with allele frequency of at least 1% (down toward 0.1% in coding regions). The final publications from phase 3 of the project report about genetic variations from more than 2500 genomes [2, 3].

Recently, several countries announced large-scale national research programs to capture the diversity of their populations, while some of these efforts started already more than 20 years ago. Since 1996 Iceland’s deCODE company is mining Icelanders’ genetic and medical data for disease genes. In 2015, deCODE published insights gained from sequencing the whole genomes of 2636 Icelanders [4]. Genome of the Netherlands (GoNL) is a whole genome sequencing project aiming to characterize DNA sequence variation in the Dutch population using a representative sample consisting of 250 trio families from all provinces in the Netherlands. In 2016, GoNL analysed whole genome sequencing data of 769 individuals and published a haplotype-resolved map of 1.9 million genome variants [5]. Similar projects have been established in larger scale in the UK: Following the UK10K project for identifying rare genetic variants in health and disease (2010–2013), Genomics England was set up in late 2012 to deliver the 100,000 Genomes Project [6]. This flagship project has by now sequenced 100,000 whole genomes from patients and their families, focusing on rare diseases, some common types of cancer, and infectious diseases. The scale of these projects is culminating in the US federal Precision Medicine Initiative, where the NIH is funding the All of Us research program1 to analyze genetic information from more than 1 million American volunteers. Even more extreme suggestions go as far as to propose “to sequence the DNA of all life on Earth”2.

The main motivation for the collection of these large and comprehensive data sets is the hope for a better understanding of genomic variation and how variants relate to health and disease, but basic research in evolution, population genetics, functional genomics and studies on demographic history can also profit enormously.

One important approach connecting evolution and functional genomics is the search for genomic regions under natural selection based on population data. The selection coefficient [7] is an established parameter quantifying the relative fitness of two genetic variants. Unfortunately, haplotype-based methods for estimating selection coefficients have not been designed with the massive genome data sets available today in mind, and may therefore take prohibitively long when applied to large-scale population data. In view of the large population sequencing efforts described above, methods are needed that—at similar sensitivity—scale to much higher dimensions.

Only recently a method for the fast computation of a genome-wide selection scan has been proposed that can be computed quickly even for large datasets [8]. The method is based on a very simple combinatorial string pattern, maximal perfect haplotype blocks. Although considerably faster than previous methods, the running time of the algorithm presented in that paper is not optimal, as it takes time in order to find all maximal perfect haplotype blocks in k genomes of length n each. This is sufficient to analyse individual human chromosomes on a laptop computer, for datasets of the size of the 1000 Genomes Project (thousands of genomes and millions of variations). However, with the larger datasets currently underway and with higher resolution it will not scale favourably. More efficient methods are therefore necessary and it was phrased as an open question whether there exists a linear-time algorithm to find all maximal perfect haplotype blocks.

In this paper we settle this open problem affirmatively. More specifically, after some basic definitions in “Basic definitions” section we present in “Linear-time method I: based on suffix trees” and “Linear-time method II: based on the positional BWT” sections two new algorithms for finding all maximal perfect haplotype blocks in optimal time. The latter of these two algorithms is then experimentally compared to the one from [8] in “Empirical evaluation” section, proving its superiority in running time by a factor of about 5 and memory usage by up to two orders of magnitude for larger data sets. “Conclusion” section concludes the paper.

This paper is an extended version of the preliminary work presented in [9]. Source code and test data are available from https://gitlab.com/bacazaux/haploblocks.

Basic definitions

The typical input to genome-wide selection studies is a set of haplotype-resolved genomes, or haplotypes for short. Clearly, for a given set of haplotypes only those sites are of interest where there is variation in the genomes. Therefore, formally, we consider as input to our methods a k × nhaplotype matrix where each of the k rows corresponds to one haplotype and each of the n columns corresponds to one variable genetic site.

Most methods distinguish only between ancestral and derived allele, reflecting the fact that most sites are biallelic. Therefore the entries in a haplotype matrix are often considered binary where the ancestral allele is encoded by 0 and the derived allele is encoded by 1. However, the computational problem and its solutions considered in this paper do not depend on this restriction and instead are applicable to any type of sequence over a constant-size alphabet .

The concept of a maximal perfect haplotype block as defined in [8] is the following, where s[i, j] denotes the substring of a string s from position i to position j and denotes the elements of an ordered set S restricted to index set K:

Definition 1

Given k sequences of the same length n (representing the rows of a haplotype matrix), a maximal perfect haplotype block is a triple (K, i, j) with , and such that

for all (equality),

or for some (left-maximality),

or for some (right-maximality), and

with such that for all (row-maximality).

Definition 1 is illustrated in Fig. 1.

Fig. 1

Illustration of Definition 1. A binary haplotype matrix with three maximal perfect haplotype blocks , and highlighted. (The example contains additional maximal perfect haplotype blocks that are not shown.)

In Cunha et al. [8] it was shown that the number of maximal perfect haplotype blocks is O(kn), while the algorithm presented there takes time to find all blocks. It is based on the observation that branching vertices in the trie of the suffixes of the input sequences starting at position p correspond to right-maximal and row-maximal blocks, while left-maximality can be tested by comparing and . In the next two sections we show how this running time can be improved.

Linear-time method I: based on suffix trees

In this section, we present our first algorithm to find all maximal perfect haplotype blocks in linear time. This solution is purely theoretical, it would likely require large amounts of memory while being slow in practice. However, it demonstrates the connection to the concept of maximal repeats in strings. We recall from [10, Section 7.12] that a maximal repeat is a substring occurring at least twice in a string or a set of strings and such that it cannot be extended to the left or to the right without losing occurrences.

Let , with the being k different characters absent from the original alphabet . The key point is that any maximal perfect haplotype block in S is a maximal repeat in . The opposite is not true: In a maximal perfect haplotype block, all occurrences of the repeat are located at the same position of each sequence of S (equality condition in Definition 1), while this constraint does not exist for maximal repeats in .

Nevertheless, finding all maximal perfect haplotype blocks in S can be performed by computing all maximal repeats in , while keeping only those whose occurrences are located at the same positions over all in which they occur. This can be done by performing the following procedure3: The key idea here is that the index characters impose that each maximal repeat occurrence starts at the same position in all sequences and, as a consequence, ensure that all occurrences occur in distinct sequences from S.

“Decorate” each sequence to create , where the index characters are symbols from an alphabet , disjoint from the original alphabet .

Find in all maximal repeats.

Any maximal repeat in with corresponds to a maximal perfect haplotype block of length , starting at position in the input sequences from S.

Hence any maximal repeat defines a unique maximal perfect haplotype block . The value |K| is the number of occurrences of r. Also the set K can be derived from occurrence positions of r in , as any position in r corresponds to a unique position in . We prefer to omit useless technical details here.

The maximal repeat occurrences in may be found using a suffix tree, constructed in time linear with respect to the size of the input data O(kn), even for large integer alphabets [12], as we have here. The maximal repeat detection is also linear with the size of the input data [10, Section 7.12.1]. Therefore the overall time complexity is O(kn).

Linear-time method II: based on the positional BWT

Here we present our second algorithm to find all maximal perfect haplotype blocks in linear time. It works by scanning the haplotype matrix column by column while maintaining the positional Burrows–Wheeler Transform (pBWT) [13] of the current column. For simplicity of presentation we assume that all rows of the haplotype matrix S are distinct. Recall that the pBWT of S consists of a pair of arrays for each column of S: For each l, , we have arrays and of length k such that the array is a permutation of the elements in the set with colexicographically (i.e. right-to-left lexicographically) sorted, and the array indicates the index from which the current and the previous rows coincide. Formally, and for all r, , we have Further let us denote by the inverse permutation of . For readers familiar with string processing terminology, the arrays and are analogous to the suffix array and the inverse suffix array, respectively, while the arrays are analogous to the LCP array.

Conditions 1, 2 and 4 (equality, left-maximality and row-maximality) of Definition 1 can be stated in terms of the arrays and as follows.

Definition 2

A quadruple (i, j; x, y) with and is called an available block if the following holds:The interval [x, y] of an available block (i, j; x, y) is called the colexicographic range of the block.

for all (equality),

there exists at least one such that (left-maximality), and

( or ) and ( or ) (row-maximality).

Lemma 1

Suppose we have a maximal perfect haplotype block (K, i, j), then the setmust be a contiguous range [x, y] of indices such that (i, j; x, y) is an available block.

Proof

This necessary condition follows immediately from Definitions 1 and 2 and the definition of the pBWT (arrays and ).

Let us consider the set of available blocks ending at column l. We have that because each available block corresponds to a distinct branching node in the trie of the reverses of , and the number of branching nodes in the trie is bounded from above by the number of leaves k. The branching nodes of the trie can be enumerated in O(k) time by using a standard algorithm [14] for enumerating LCP intervals of the LCP array of the trie, . This gives us the colexicographic ranges [x, y] of all available blocks in . An example is shown in Fig. 2.

Fig. 2

Available blocks. Left: an example of a haplotype matrix up to column 6 with the two arrays and on the right. Center: the colexicographically sorted rows and the array listed on the right. Right: the trie of the reverses of the rows of the matrix. For example, the block is available because , , , is the consecutive range , we have for all with , and we have and . The repeat in the block is 00, and we see it is a branching node in the trie on the right

The only thing left is to show how to check the right-maximality property of an available block. The following lemma gives the sufficient condition for this:

Lemma 2

An available block (i, j; x, y) corresponds to a maximal haplotype block (K, i, j) if and only ifor.

If , right-maximality according to Definition 1 holds trivially. If , right-maximality requires that there are two rows for which . Since all rows s, t qualifying for this condition are within the colexicographic range [x, y] of our available block, the statement follows immediately.

To check the condition of Lemma 2 in constant time for , we build a bit vector such that and if and only if . Now the block is right-maximal if and only if contains at least one 1-bit. We can build a vector of prefix sums of to answer this question in constant time.

Time and space complexity

We assume the column stream model, where we can stream the haplotype matrix column by column. We can thus build the arrays , and on the fly column by column [13], and also easily build the required prefix sums of arrays from these. The time is O(nk), since each of the n columns takes O(k) time to process. The algorithm needs to keep in memory only the data for two adjacent columns at a time, so in space O(k) we can report the colexicographic ranges of all maximal blocks ending in each column . If the colexicographic range of a block at column l is [x, y], then the rows in the original haplotype matrix are . There are O(nk) blocks and O(k) rows per block, so the time to report all rows explicitly is . In fact, a sharper bound that can also easily be achieved is where is the size of the output. Alternatively, we can store a complete representation of the answer taking O(nk) space by storing all the arrays and the colexicographic ranges of the maximal perfect blocks for each column, from which we can readily report all rows in any maximal perfect block in constant time per row.

Empirical evaluation

Since the algorithm of “Linear-time method I: based on suffix trees” section is mostly of theoretical interest, we evaluate only the pBWT-based algorithm presented in “Linear-time method II: based on the positional BWT” section. The source code is available from https://gitlab.com/bacazaux/haploblocks. As a baseline for comparison we use the implementation of the trie-based algorithm by Cunha et al. [8], available from the same gitlab site. The experiments were run on a machine with an Intel Xeon E5-2680 v4 2.4 GHz CPU, which has a 35 MB Intel SmartCache. The machine has 256 gigabytes of memory at a speed of 2400MT/s. The code was compiled with g++ using the -Ofast optimization flag.

Our test data consists of chromosomes 2, 6 and 22 from phase three of the 1000 Genomes Project [2], which provides whole-genome sequences of 2504 individuals from multiple populations worldwide. We preprocessed the data by extracting all biallelic SNPs from the provided VCF files4 and converting them to a binary haplotype matrix using our own program vcf2bm, also available from https://gitlab.com/bacazaux/haploblocks.

Our implementation has a user-defined parameter allowing to adjust the minimum size of a reported maximal perfect haplotype block (K, i, j), where size is defined as the width () times the number of rows (|K|) in the block. Table 1 shows the running times and memory usage of our implementation on the different chromosomes and for different settings of the minimum block size parameter. The larger the minimum block size, the faster the algorithm is, because there are less blocks to report. In general, it takes only a few minutes to process a complete human chromosome. Locating all 323,163,970 blocks of minimum size in all 22 human autosomes (non-sex chromosomes) took in total 4 h and 26 min with a memory peak of 12.8 MB (data not shown).

Table 1

Running times and memory usage of our pBWT-based implementation

Data set	#lines	#columns	Min block size	Time	Memory (MB)	#blocks
chr. 22	5008	1,055,454		4 min 54 s	12.8	148,613,645
chr. 22	5008	1,055,454	500,000	3 min 50 s	12.8	16,076,453
chr. 22	5008	1,055,454	1,000,000	3 min 40 s	12.8	2,228,762
chr. 22	5008	1,055,454	2,000,000	3 min 43 s	12.8	4779
chr. 6	5008	4,800,101		19 min 42 s	12.8	624,689,548
chr. 6	5008	4,800,101	500,000	17 min 20 s	12.8	89,840,467
chr. 6	5008	4,800,101	1,000,000	16 min 30 s	12.8	11,388,982
chr. 6	5008	4,800,101	2,000,000	16 min 36 s	12.8	5585
chr. 2	5008	6,786,300		31 min 57 s	12.8	946,717,897
chr. 2	5008	6,786,300	500,000	25 min 06 s	12.8	160,094,115
chr. 2	5008	6,786,300	1,000,000	23 min 24 s	12.8	25,533,314
chr. 2	5008	6,786,300	2,000,000	23 min 18 s	12.8	120,243

Note that in our streaming implementation the memory usage is dominated by the number of haplotypes times the buffer size, and therefore is essentially constant in this study

Table 2 shows a comparison of our implementation to the trie-based implementation from [8]. Our implementation is about 5 times faster on all datasets, and the memory consumption is up to 93 times smaller.

Table 2

Comparison of the trie-based implementation from [8] and our pBWT-based implementation with minimum block size

Data set	trie		pBWT
	Time	Memory	Time	Memory (MB)
chr. 22	17 min 08 s	927.8 MB	3 min 40 s	12.8
chr. 6	1 h 34 min 34 s	3.23 GB	16 min 30 s	12.8
chr. 2	2 h 07 min 21 s	4.46 GB	23 min 24 s	12.8

It is now easy to apply the method for estimating a local selection coefficient from the size of maximal perfect haplotype blocks covering a certain genomic region presented in [8]. This method estimates the likelihood of observing a haplotype block for a given selection coefficient s and the time t since the onset of selection following an approach presented by Chen et al. [15]. Therefore, chromosome-wide selection scans indicating the loci of maximum selection, as shown in Fig. 3 for the complete human chromosome 2 (size parameter ), can now be generated in less than half an hour.

Fig. 3

Selection scan for human chromosome 2. Shown is for each position of the chromosome the largest maximum likelihood estimate derived from any maximal perfect haplotype block overlapping that locus. It is easy to spot potential regions of high selection. The centromere, located around 93 Mbp, shows no signal as sequencing coverage is low here and no SNPs could be called

Conclusion

In this paper we presented two algorithms that are able to find all maximal perfect haplotype blocks in a haplotype matrix of size in linear time O(kn). In particular the second method, based on the positional Burrows–Wheeler Transform, performs also extremely well in practice, as it allows for a streaming implementation with extremely low memory footprint.

While an initial implementation of the method is available from https://gitlab.com/bacazaux/haploblocks, a user-friendly software combining the algorithm presented here with the computation of the selection coefficient suggested in [8] remains to be developed.