FastTagger: an efficient algorithm for genome-wide tag SNP selection using multi-marker linkage disequilibrium

Multi-marker tagging rules

Most SNPs have only two alleles, so we consider only bi-allelic SNPs. Given a population, the allele with higher frequency in the population is called major allele, and the allele with lower frequency is called minor allele. We use uppercase letters to denote the major alleles of SNPs, and use lowercase letters to denote the minor alleles. SNPs that are far apart from each other usually are not linked. Here we require that the distance between every pair of SNPs in a rule must not exceed a predefined distance threshold max_dist.

Given k SNPs S = {SNP₁, SNP₂, ⋯, SNP _k }, there are 2 ^k possible haplotypes over the k loci. To calculate the r² statistic of rule S → SNP _x , we need to divide the 2 ^k haplotypes into two non-empty groups and map the two groups to the two alleles of SNP _x . MultiTag [19] and MMTagger [21] uses different methods to do the mapping.

The co-occurrence model

MMTagger does the mapping based on the co-occurrences of the alleles of the SNPs on the left hand side and the alleles of the SNP on the right hand side. Let H be a haplotype over the SNP set S on the left hand side, A and a be the two alleles of SNP _x on the right hand side, and f(H) be the frequency of H. We use f(HA) to denote the frequency of H and SNP _x = A occurring together, and f(Ha) to denote the frequency of H and SNP _x = a occurring together. If f (HA) > f (Ha), we map haplotype H to allele A of SNP _x , otherwise we map haplotype H to allele a of SNP _x . Let H _A be the set of haplotypes mapped to allele A, and H _a be the set of haplotypes mapped to allele a. We convert SNP set S to a bi-allelic marker with two "alleles" H _A and H _a . Then we can calculate the r² statistic between S and SNP _x as follows.

(1)

where P(H _A ), P (H _a ), P (A), P (a) and P (H _A A) are the relative frequencies of H _A , H _a , A, a and H _A A respectively.

We implemented both models in the FastTagger algorithm, and let users choose which model they want to use.

If the r² statistic between S and SNP _x is no less than a predefined threshold min_r 2, we say that SNP _x can be tagged by S, and R : S → SNP _x is a tagging rule. With the increase of the size of S, the haplotypes of S partition the whole dataset into finer and finer groups. In an extreme case, every haplotype of S occurs at most once. In this case, the association between haplotypes of S and alleles of SNP _x becomes unreliable. To prevent over-fitting, we put a constraint on the size of S. The size of S should not exceeds a predefined threshold max_size.

The r² statistics can be calculated from phased haplotype data directly. If the SNP data are in the form of unphased genotype data, we can use existing haplotype inference algorithms such as PHASE [22] to convert genotype data into phased haplotype data. We can also estimate k-marker haplotype frequencies directly from genotype data without phasing using the algorithms described in [23, 24]. The second approach is used in algorithm LD-select [9].

Generating tagging rules

To generate all the tagging rules, we need to enumerate all the SNP sets that satisfy the maximum distance constraint and maximum size constraint, and then calculate the r² statistics between these SNP sets and their nearby SNPs. The search space can be enormously large when the number of SNPs is large. We use several techniques to reduce the number of rules to be tested.

Selecting tag SNPs using a greedy approach

Finding the smallest set of tag SNPs is computationally expensive. FastTagger uses a greedy approach similar to the one proposed in [9, 19] to find a near optimal set of tag SNPs.

Let C be the set of candidate tag SNPs, T be the set of tag SNPs selected, and V be the set of SNPs not being covered. A SNP is covered if either it is a tag SNP or it can be tagged by some SNP set S such that S ⊆ T. Initially, C and V contain all the SNPs, and T is empty.

FastTagger first identifies those SNPs that do not appear at the right hand side of any tagging rules, and these SNPs must be selected as tag SNPs. FastTagger puts them into T and remove them from C. These SNPs are also removed from V. For the remaining SNPs in V, if they can be tagged by some SNP set S such that S ⊆ T, then they are removed from V too.

Next, for each SNP SNP _i ∈ C, FastTagger finds the set of SNPs in V that are covered by SNP _i . A SNP SNP _j in V is covered by SNP _i if SNP _j is not tagged by any subsets of T and there exists a subset S of T such that SNP _j is tagged by S ∪ {SNP _i }.

FastTagger then picks a SNP from C that covers the largest number of SNPs in V as a tag SNP. This newly picked tag SNP is put into T and removed from C. All the SNPs that are covered by it including itself are removed from V. This process is repeated until V is empty, that is, all the SNPs have been covered. In each iteration, in order to find the set of SNPs covered by every candidate tag SNP in C, FastTagger needs to keep the tagging rules in memory. However, the number of rules generated can be very large. It is possible that the total size of tagging rules is too large to fit into the main memory. To solve this problem, we can break the whole chromosome into several chunks such that the rules over every chunk can fit into the main memory. We then select tag SNPs within each chunk.

When selecting tag SNPs within each chunk, only those tagging rules whose SNPs all fall into this chunk are used. To also utilize the rules across chunks, we allow two adjacent chunks to have certain overlap. The length of the overlap is determined by the max_dist threshold. The SNPs in one chunk that are within max_dist bases away from the first SNP of the next chunk are included in the next chunk since they can tag or be tagged by SNPs in the next chunk. FastTagger finds tag SNPs from each chunk from left to right. The tag SNPs selected in the current chunk that also belong to the next chunk will be passed on to the next chunk as tag SNPs. Note that if the distance between two adjacent SNPs is larger than max_dist, then these two SNPs are used as a breakpoint even if there is enough memory. The reason being that if the distance between two adjacent SNPs is larger than max_dist, then the two SNPs cannot tag each other or each other's neighbors.

Using the above method, FastTagger can work on chromosomes containing more than 100 k SNPs with as less as 50 MB memory, while existing algorithm consumes more than 1 GB memory even on chromosomes containing around 30 k SNPs.

Comparison with other algorithms

The first experiment is to compare FastTagger with LRTag [13], MMTagger [21] and MultiTag [19]. LRTag uses only pair-wise LD to find tag SNPs, and it has been shown to outperform LD-select and FESTA. Hence we choose LRTag as a representative of the pairwise algorithms. MMTagger and MultiTag both use multi-marker LD to find tag SNPs. We obtained the programs from their respective authors. FastTagger used all the techniques described previously except the skipping rules technique. LRTag takes pre-computed pairwise r² statistics as input, so the running time of LRTag includes only tag SNP selection time. We report the results at min_r 2 = 0.95 here, results at min_r 2 = 0.9 and min_r 2 = 0.8 can be found in supplementary materials. For all the four algorithms, the selected tag SNPs can cover the whole region of interest.

Table 3 shows that the multi-marker model can reduce the number of tag SNPs significantly under the same min_r 2 threshold compared with the pairwise model (Table 2). The number of tag SNPs is reduced by more than 30% when max_size = 2. When max_size = 3, the number of tag SNPs is reduced by more than 40%. However, calculating multi-marker r² statistics is much more expensive than computing pairwise r². FastTagger is more than 10 times slower when max_size = 2, and hundreds of times slower when max_size = 3.

On ENCODE regions, FastTagger and MMTagger take similar time to finish when max_size = 2; when max_size = 3, FastTagger is 2-3 times faster than MMTagger. On the 6 chromosomes, FastTagger is 2-6 times faster than MMTagger. Both algorithms are orders of magnitude faster than MultiTag. The number of tag SNPs selected by FastTagger under the co-occurrence model is smaller than that selected by MMTagger and MultiTag.

The effectiveness of the techniques used in FastTagger

Portability and prediction accuracy

Multi-marker models group combinations of the alleles on the left hand side into two groups, and then map these two groups to the two alleles on the right hand side. Compared with pairwise model, multi-marker models are more prone to over-fitting. Here we use three populations in HapMap--the Han Chinese population (HCB), the Japanese population (JPT) and the Caucasian population(CEU)--to study the portability and prediction accuracy of tagging rules of different lengths. We use chr19 in this experiment. We first generate tagging rules from one population, and then calculate the r 2 statistics and prediction accuracy of these rules in the other populations. The prediction accuracy of a rule is defined as the proportion of alleles of the SNP on the right hand side that are correctly predicted by the alleles of the SNPs on the left hand side. The results reported below are results when rules are generated from individuals in the Han Chinese population and are evaluated using individuals in the other two populations. In all three populations, we consider only those SNPs with MAF ≥ 5%.

Figures 1, 2 and 3 show the distribution of the r² values of the rules generated from the Han Chinese population using the two multi-marker models in the three populations. Table 11 shows average r² of the rules in the three populations. For all the three lengths, the average r² of the rules in the Japanese population and the Caucasian population is lower than that in the Chinese population. The decrease of length-2 and length-3 rules is more significant than that of length-1 rules, which indicates that longer rules are more prone to over-fitting than shorter rules for both models. The r² values of the rules become much lower in the Caucasian population than that in the Japanese population, which is consistent with the genetic differences between the three populations.

		#rules			average r2			average accuracy
len	model	HCB	JPT	CEU	HCB	JPT	CEU	HCB	JPT	CEU
1	pairwise	85961	84123	69083	0.978	0.942	0.865	0.995	0.989	0.966
2	co-occurrence	1563176	1472654	1014934	0.965	0.878	0.745	0.993	0.977	0.938
2	one-vs-the-rest	1560181	1469765	1012699	0.965	0.881	0.753	0.993	0.977	0.940
3	co-occurrence	26182522	24495802	16064120	0.952	0.790	0.665	0.990	0.960	0.913
3	one-vs-the-rest	7074493	6269985	3955224	0.970	0.791	0.659	0.994	0.970	0.919

The rules are generated from Han Chinese population with min_r 2 = 0.9. Some rules may become invalid in the other two populations because the MAF of some SNPs in the other two populations may be smaller than 5%. When only pairwise LD is used, all algorithms generate the same set of rules. When multi-markers are considered, FastTagger-COOC and MMTagger generate the same set of rules using the co-occurrence model; FastTagger-avsR and MultiTag generate the same set of rules using the one-vs-the-rest model.

The same trend is observed on prediction accuracy (Figure 4, 5 and 6). Even though the rules are generated from the Chinese population, their accuracy in the Japanese population is always above 80%. Even for length-3 rules, 94% rules generated using the co-occurrence model have an accuracy no less than 90%, and 97.4% rules generated using the one-vs-the-rest model have an accuracy no less than 90% in the Japanese population. The average accuracy of length-3 rules is above 96% for both models in the Japanese population(Table 11). The average accuracy of the rules in the Caucasian population is lower than that in the Japanese population, but it is still above 91% even for length-3 rules. We believe that if we use individuals from the same population to do the testing, the average r² and accuracy should be even higher. As for the two models, the number of length-2 rules generated by the two models is similar, while the co-occurrence model generates about 3.5 times more length-3 rules than the one-vs-the-rest model. The average r² and accuracy of the length-3 rules generated using the one-vs-the-rest model is higher than that generated using the co-occurrence model on both populations. However, since much less rules are generated under the one-vs-the-rest model, the one-vs-the-rest model needs more tag SNPs to cover all the other SNPs than the co-occurrence model as shown in Table 3.