Cgaln: fast and space-efficient whole-genome alignment

Overview

Figure 1 shows the flow of Cgaln. First, Cgaln divides the input sequences into "blocks" with a fixed length (Figure 1A). These blocks are taken as "elements" to be aligned at the block level. Each cell of the meshlike structure is associated with a block-to-block similarity score. The similarity between two blocks, each from the two input sequences, is evaluated by the frequency of seeds (k-mers) commonly found in the blocks. Similar methods based on seed counts have been used to quickly estimate the degree of similarity between two protein sequences [22, 23]. Based on these similarity scores, the block-level alignments are obtained by a local DP algorithm (Figure 1B). For the local DP, we apply the Smith-Waterman algorithm [24] modified so that sub-optimal similarities are also reported [25]. The size of a block is arbitrary within the limits that both the number of blocks and the size of each block should be < 2¹⁶ in order to be represented by unsigned short integers.

The nucleotide-level alignment is conducted on the restricted regions included in the block-level alignment found in the first stage (Figure 1C). We adopt a seed-extension strategy widely used in homology search programs such as Blast [26]. The obtained seed matches are integrated into gapless HSPs. The HSPs are then filtered and chained to conform to coherent alignment. If necessary, the chained HSPs may serve as anchor points, the subsequences between which are aligned by a standard DP algorithm or recursive search for shorter seeds.

The above-mentioned algorithm implicitly assumes that genomic rearrangements smaller than the block size of Cgaln (about 10 Kb) are rare, and they are not searched for with the default settings. To discover rearrangements smaller than the block size, we prepared the "-sr" option that considers inversions as small as a few hundred bases at the recursive phase applied to inter-HSP regions (see below).

Cgaln accepts two single- or multi-fasta files. When either or both files are in multi-fasta form, Cgaln aligns every pair of single-fasta entry sequences and outputs the united results. In this section, we assume for simplicity that both input sequences are single-fasta files. By setting the "-r" option, Cgaln examines both orientations of a sequence in a single run.

Seed designing

Cgaln uses the "spaced seed" proposed in PatternHunter [27] for fast and sensitive seed matching. A k/w spaced seed is a discrete series of nucleotides of length w in which k <w positions are examined for nucleotide matching. Cgaln uses an 11/18 spaced seed by default, with the same pattern as that used in PatternHunter (expressed as 111*1**1*1**11*111, where "1" and "*" respectively indicate the positions to be examined or ignored). This pattern was designed to be most sensitive for a pair of randomly generated sequences with 70% nucleotide identities. Although some other seed patterns are more sensitive than this pattern under some conditions (e.g., for coding regions), we chose this seed design so that Cgaln could be generally applicable to whole genome sequences of various species. Optionally, a 12/19 or 13/20 spaced seed can be used, as suggested by Mak and Benson [28] (expressed as 1111*1*1**11**1*111 and 1111*1**11**11*1*111, respectively). We use the term k-mer hereafter to refer to these k/w spaced seeds, as the value for k is most important. A larger k is appropriate for longer sequences but requires more memory, and hence there is a compromise in the choice of the best value for k. The "code" of each seed is its quaternary expression calculated from the k examined positions by converting nucleotides A, C, G, and T into numerals 0, 1, 2, and 3, respectively. Thus, the code of a k/w seed is between 0 and 4 ^k - 1.

Block-level alignment (BA)

In this subsection, we describe block-level alignment (abbreviated as BA hereafter).

Measuring the score between two blocks

Let us denote the given input genomic or chromosomal (single-fasta) sequences G _a and G _b . Let L _a and L _b be the lengths of G _a and G _b , respectively. First, Cgaln divides G _a and G _b into blocks with the fixed length of J, except that the last block may be shorter than J. The numbers of blocks in G _a and G _b are denoted as m _a and m _b , and thus m _a = ⌈L _a /J⌉ and m _b = ⌈L _a /J⌉. Let be the x-th block of G _a and be the y-th block of G _b (1 ≤ x ≤ m _a , 1 ≤ y ≤ m _b ).

The score M _x,y , the measure of similarity between and , is defined as the summation of the scores of seeds commonly found in both and . The score of seed k _i is evaluated by the probability of the chance of finding one or more matches of k _i in a block pair. Let be the total number of k _i in G _a , then the expected mean number of occurrences of k _i in a block is

(1)

We assume the Poisson distribution for the probability of k _i occurring f _i times in a block , i.e.,

(2)

p _b (f _i ) is obtained analogously. When k _i occurs f _i times in and h _i times in , k _i matches f _i h _i times between and . It is clearly an overestimation to use this number as the measure of similarity between and , because true homologous matches should line up around a single diagonal. In such a case, we regard min(f _i , h _i ) as the number of matches as suggested by [23]. Thus, the score of seed k _i is

(3)

where the probability of seed k _i occurring more than or equal to x times in a block is

(4)

Summing up s(k _i ) for all k-mers, we obtain the similarity score as:

(5)

Local DP alignment at block level

After obtaining M _x,y for (1 ≤ x ≤ m _a and 1 ≤ y ≤ m _b ), the local alignment of BA is conducted by using local DP as follows:

(6)

where d is the gap penalty. In practice, M_{x, y}is calculated in the process of this DP procedure to save the required memory from O(m _a m _b ) to O(min(m _a , m _b )). To obtain the optimal and suboptimal locally best-matched alignments, we use the algorithm proposed by Gotoh [25]. This algorithm assigns each block cell to a "colony", which is a candidate for local alignment. To define the colony borders, we apply the X-drop-off approach [29]. At a certain block cell, if the score F_{x, y}increases from 0 to a positive value, a new colony starts. A colony is "significant" if F_{x, y}> T _col somewhere in the colony, where T _col denotes a predefined threshold. A colony ends when F_{x, y}becomes ≤ 0 or falls by more than T _col from the maximum score of the colony, then F_{x, y}is reset to 0. Only significant colonies are retained for further analysis. This method can greatly reduce the storage requirement, while the computational time remains O(m _a m _b ).

It should be noted that M_{x, y}is always ≥ 0 even between nonhomologous blocks, while it is a prerequisite that, on average, the similarity score must be less than zero for the local DP algorithm to work properly [30]. Consequently, we subtract a constant bias B from each M_{x, y}, so that , i.e., the X-drop-off should finish after passing through, on average, a unrelated cells (a = 5 by default). If there are repetitive sequences that escaped masking, BA may output them many times. Cgaln has an option for suppressing such repetitious outputs. If there are inconsistent colonies that overlap each other, the "-fc" option filters all of them out except for the one with the highest scored.

Nucleotide-level alignment (NA)

Cgaln applies the nucleotide-level alignment (abbreviated as NA hereafter) within the restricted areas that are composed of block cells included in the local alignments of BA. However, the area covered by the set of block cells thus obtained is insufficient for NA, because the expected results of NA may shift slightly away from the aligned cells. Hence, we extend the area to be searched by NA to the "envelope" of the block cells found in the first stage (see Figure 2A).

Evaluation method

It is difficult to evaluate the accuracy of genomic alignment because of the lack of "true" alignment data [34]. In this paper, we focused our attention on the orthologous protein-coding genes and corresponding coding exons. This approach has obvious drawbacks, as most genomic alignment programs, including Cgaln, are designed to find not only orthologs but also other homologs. Moreover, alignment of non-coding genes and intergenic regions can be misinterpreted in our procedure. However, we expect that this kind of imprecision would not much affect the evaluation of the relative performance of various methods.

Figure 3 schematically illustrates a case of alignment results. For a bacterial genome pair, we considered how many homologous base pairs in the reference alignment of orthologous gene pairs are correctly aligned by a given method for calculating sensitivity and specificity. A true positive (TP) indicates the number of genomically aligned bases that coincide with those in the reference alignment, a false negative (FN) indicates the number of bases in the reference alignment not observed in the genomic alignment, and a false positive (FP) indicates the number of genomically aligned bases outside the reference alignment. Thus, the sensitivity (Sn) and specificity (Sp) are defined as:

(7)

For the example shown in Figure 3,

(8)

and

(9)

For a mammalian chromosomal pair, sensitivity is evaluated as before, in which homologous exons are considered homologous regions. However, it is dificult to evaluate specificity because the biologically significant regions that hold large amounts of the entire genome are not clear. In this examination, we counted the total length of generated HSPs ((b) + (c) + (e) + (f) in Figure 3) as an indicator of specificity, for if the total HSP length from one alignment program is considerably greater than that from other programs with the same sensitivity, the program may be judged to have low specificity. We also evaluated specificity with a human chromosome 20 - mouse chromosome 2 pair. Human chromosome 20 is considered completely orthologous to mouse chromosome 2 [35], and the alignment result of human chromosome 20 and mouse chromosome 2 should coincide with that of human chromosome 20 and the whole genomic sequence of mouse. Therefore, specificity can be evaluated as

(10)

where len_20-2 is the total length of HSPs for human chromosome 20 vs. mouse chromosome 2, while len_20-allis that for human chromosome 20 vs. whole mouse genome.

Preparation of data

We obtained the whole-genome sequences of Escherichia coli CFT073 (Accession No. NC_004431.1, 5,231,428 bp), E. coli K12 (NC_000913.2, 4,639,675 bp), and Salmonella typhimurium (NC_003197.1, 4,857,432 bp) from NCBI http://www.ncbi.nlm.nih.gov/. Before applying these genomic sequences to alignment programs, we masked repetitive sequences by MaskerAid [36] with default parameters. We also prepared the whole-genome sequences of human (build hg19) and mouse (build mm9) from the UCSC Genome Browser http://genome.ucsc.edu/. These human and mouse genome sequences were already soft-masked, and we did not use any masking tool. We also downloaded the latest version (KOREF_20090224) of Korean individual (SJK) genomic sequence [37] from the FTP site ftp://ftp.kobic.kr/pub/KOBIC-KoreanGenome/.

Our examination requires reference data on the locations of homologous regions on the input sequences. To obtain such a reference dataset, we first collected sets of orthologous gene pairs. For bacterial genomes, we used MBGD (Microbial Genome Database for Comparative Analysis) [38] for this purpose. MBGD is a database for comparative analysis of microbial genomes, and possesses data on orthologous gene clusters of bacteria. While COG [39] is well known as an orthologous gene database, we preferred MBGD because MBGD is better than COG at clustering orthologous genes in more detail. Data on orthologous gene pairs between human and mouse were obtained from RefSeq http://www.ncbi.nlm.nih.gov/RefSeq/ and Ensembl http://www.ensembl.org/. Next, we aligned corresponding cDNA sequences of all gene pairs by the standard local DP algorithm [24]. Finally, we located each gene on the reference genome. For bacterial pairs, the information on location was obtained from GenBank. For mammalian pairs, we mapped the cDNA sequences on the respective genomic regions by Spaln [40, 41]. From 10042 original gene pairs, 9373 pairs (= 18746 genes) could be mapped on the reference genomes for the homologous genomic regions. As 97% (18252) cDNAs were mapped on the reference genomes with 100% identity and the other genes were mapped with more than 95% identity, we regarded these data as valid. By combining the cDNA alignment and the mapping coordinates, we obtained the homologous exonic regions on the genomes.

Programs used for tests

To compare the performance of our algorithm with other leading programs presently available, we developed the Cgaln program that implements the above-mentioned algorithm in C on a Linux platform. We report on comparisons of the accuracy and computational speed of Cgaln with those of Blastz, AVID, and NUCmer. Blastz is one of the principal pairwise alignment programs for long sequences, and is used as an internal engine of several multiple genomic sequence alignment programs, such as MultiPipMaker [42], TBA, and MultiZ [43]. We examined Blastz with two sets of parameter values; with the default parameter set and with the tuned parameter set (T = 2, C = 2). The option "T = 2" disregards transitions as matches; this speeds up computation but slightly reduces sensitivity. The option "C = 2" directs "chain and extend", which helps to reduce false positives. AVID is also a fast and accurate genomic alignment program, but it is a global aligner and not suitable for the alignment of genomes with large rearrangements. We applied AVID only to bacterial genome pairs. NUCmer is a variant of MUMmer 3.0. It clusters the matches of MUMmer 3.0 and tries to align the non-exact regions between the matches by DP.

Although there are other genomic alignment tools, such as CHAOS, LAGAN, SLAGAN, DIALIGN, GLASS, and WABA, they are either too slow to execute or their source codes are not available. All experiments were performed on a 2.0 GHz Core2Quad (64-bit CPU) with 8 GB memory.

Parameter tuning

Preliminary examinations indicated that the performance of Cgaln depends strongly on the outcome of BA, and hence a proper choice of parameter values at this level is essential, especially for distantly related genome pairs. We tested Cgaln with various block sizes and k-mer sizes, and found that, if appropriate threshold values for "significant" colonies T _col are used, there is no remarkable difference in accuracy or computational time in a wide range of block sizes. However, if a block is too large, computation will be slow because a lot of spurious HSPs tend to be generated within block cells. Thus, block sizes between 5,000 bp and 30,000 bp are proven to be almost equally appropriate for both bacterial and mammalian genomes. It is desirable to automatically estimate the optimal T _col value for a given set of block sizes, k-mer size, and overall divergence between the sequences under comparison. In this study, we first tuned T _col for default block size and k-mer size (expressed by T_col-default), and then derived an empirical rule applicable to other block sizes:

(11)

where r _b is the ratio of the given block size to the default block size.

For a given k-mer size, we use the relation:

(12)

where k _{dif f} is the difference between the given k-mer size and the default k-mer size. Although these empirical rules work satisfactorily, a theoretically supported tuning algorithm remains to be established. We set the default value as follows: block size = 10000, T_col-default= 3000, gap penalty at BA d = T_col-default/15. These parameters can be changed by selecting various options.

Alignment of bacterial genomes

We first examined the performance of Cgaln in comparison with Blastz, AVID, and NUCmer for two sets of pairwise alignments of bacterial genomes: (E-E) E. coli CFT073 vs. E. coli K12 (2531 gene pairs), and (E-S) E. coli CFT073 vs. S. typhimurium (3385 gene pairs).

Blastz with the default parameters produced a lot of spurious alignments, which resulted in high sensitivity but low specificity. With tuned parameters, Blastz improved the specificity and computational time but decreased sensitivity. AVID, a global alignment tool, shows high sensitivity and specificity in the intra-species comparison (E-E). However, it cannot consider inversion and is proven to be insensitive and not specific in interspecies comparisons (E-S). Moreover, AVID consumes a lot of memory, about 2.2 GB. NUCmer is fast and more memory-efficient than Cgaln but is much less sensitive. In (E-E), Nucmer could identify about 82% of all genes, but the coverage of almost all genes was low (~50%), possibly due to the use of long seeds (20-mer), and so it is applicable only to closely related sequence comparisons. Table 1 and Table 2 show that Cgaln is the fastest and second most memory-efficient among these programs, with high accuracy. Compared with Blastz with tuned parameters, Cgaln is twice as fast and more space-efficient, comparably sensitive, but less specific. This inferiority in specificity is due to overalignment at BA. Cgaln with the "-fc" option (see Methods) could improve specificity to a level comparable to that of tuned Blastz. However, it is not always better to use this option, as there is a trade-off between sensitivity and specificity.

Alignment of mammalian chromosomes

We also compared the performance of Cgaln with that of Blastz and that of NUCmer on two kinds of mammalian homologous chromosome pairs: (H20-M2) human chromosome 20 (63,025,520 bp) vs. mouse chromosome 2 (181,748,087 bp), and (X-X) human chromosome X (155,270,560 bp) vs. mouse chromosome X (166,650,296 bp). The numbers of orthologous gene pairs on (H20-M2) and (X-X) are 278 and 341, respectively.

The results are summarized in Table 3 and Table 4. We applied Blastz to two types of input sequences: (i) intact chromosomal sequences, and (ii) chromosomal sequences split into chunks of 10 Mb with 10 Kb overlaps. The third rows of Table 3 and Table 4 labeled with "*" show the results obtained with the split sequences. The time for splitting and uniting is not included. Blastz with the default parameters resulted in nearly the same accuracy in the cases of (i) and (ii) above. This is not surprising, because Blastz with the default parameters does not chain HSPs. As in the case of bacterial genomes, the results of Blastz with the default parameters were sensitive but not specific, and the computation was slowest among all of the programs and settings examined. The low specificity (56.2%) and the large HSP length (especially in (X-X)) indicate that Blastz generates a lot of repetitious HSPs, although we used masked sequences. The "C = 2" option improved the specificity of Blastz drastically, and the "T = 2" option reduced the computational time with a slight decrease in sensitivity. However, when applied to the intact sequence pair of chromosome X, the tuned Blastz missed some long homologous regions (shown in Figure 4A), which resulted in poor exon coverage (36.5%). With the split sequences, this deficit was not observed (Figure 4B), while specificity declined because filtering by the chaining process was not sufficient. In the case of the chromosomal pair (H20-M2), such a big deficit was not observed, possibly because of the high similarity and the small number of rearrangements between the two chromosomes. The consumed memory of Blastz was large, especially for (X-X) (2.8 GB) when intact chromosomes were examined, regardless of the choice of options. NUCmer was as fast and memory-efficient as Cgaln in examination (H20-M2), but consumed twice as much memory as Cgaln in examination (X-X). Moreover, NUCmer was shown not to be sensitive in the pairwise alignment of mammalian chromosomes.

	length (bp)	Sn (%)	Sp mam (%)	time	memory (GB)
Blastz (default)	18,598,895	85.9	57.5	66 m 40 s	1.6
Blastz (T = 2 C = 2)	16,353,601	85.1	95.5	9 m 45 s	1.6
Blastz* (T = 2 C = 2)	16,665,937	85.3	80.8	12 m 20 s	0.4
NUCmer	1,118,494	5.8	75.2	3 m 55 s	1.0
Cgaln (-X2500 k = 11)	13,964,626	79.2	92.3	2 m 14 s	0.8
Cgaln (-X2500 k = 12)	15,154,530	81.4	90.6	2 m 16 s	1.1

Performance of alignment programs examined on human chromosome 20 (chr20) vs. mouse chromosome 2 (chr2). The row of Blastz with the symbol "*" refers to the examination on split chromosomal sequences. The threshold T _col was set to "-X2500". Sp _mam is evaluated as the ratio of the coverage with human chr20 - mouse chr2 to that with the human chr20 - whole-genome sequence of mouse.

	length (bp)	Sn (%)	time	memory (GB)
Blastz (default)	41,649,203	58.9	117 m 03 s	2.8
Blastz (T = 2 C = 2)	21,094,947	36.5	15 m 21 s	2.8
Blastz* (T = 2 C = 2)	32,684,714	58.0	23 m 22 s	0.4
NUCmer	1,299,148	5.4	4 m 01 s	2.3
Cgaln (-X2500 k = 11)	25,675,191	54.4	5 m 56 s	0.9
Cgaln (-X2500 k = 12)	27,776,861	56.0	4 m 21 s	1.2

Performance of alignment programs examined on human chromosome X (chrX) vs. mouse chromosome X. The threshold T _col was set to "-X2500".

We examined Cgaln with two kinds of k-mer spaced seed, 11-mer and 12-mer. Cgaln with 12-mer required more memory than that with 11-mer, but much less than Blastz. Cgaln performed better with the 12-mer spaced seed than with the 11-mer. In fact, several genes were missed with the 11-mer but identified with the 12-mer. With the 11-mer, the number of occurrences of each seed in a block often exceeds the proper range suitable for scoring by Poisson distribution, which we consider the main reason for the lower performance with the 11-mer. It should be noted that with either k-mer, Cgaln does not generate much noise like Blastz with default parameters, nor does it cause a big deficit like Blastz with tuned parameters for the intact chromosomal pair (X-X).

However, Cgaln was slightly less sensitive than Blastz because the nucleotide-level coverage of gene pairs in the former was slightly worse in the latter. Scrutinizing the computational processes, we found that the differences in gene coverage between Cgaln and Blastz originate mainly when HSPs are extended with gaps. In fact, we confirmed that the sensitivity of Cgaln was slightly improved when the X-drop-off threshold for gapped extension of HSP or the length-threshold for DP (T _dp ) was augmented. The slightly smaller total length of HSPs compared with those of Blastz also indicated that Cgaln "underaligns" the sequences. However, we question the significance of the difference in coverage, because Blastz might "overalign" sequences, as the X-drop-off parameter of Blastz is so large that it may improve nominal sensitivity but may align some non-homologous regions [44].

Application of Cgaln to human and mouse whole genomes

We also applied Cgaln to whole genomic sequences of human and mouse to investigate how many homologous genes Cgaln had caught. In this examination, we set a threshold parameter of T _col = 3000 and used the 12-mer seed. As the average nucleotide identity levels vary considerably among chromosomal pairs between human and mouse, it might be preferable to change T _col values depending on the chromosome pairs. However, for simplicity we used a fixed value in this examination. The computation took 750 m with 1.8 GB of memory on a desktop computer. Cgaln generated HSPs totaling 693,583,236 bp with 73.6% average nucleotide identities. A total of 8897 gene pairs (95%) were identified, including 7229 pairs (77%) with coverage of more than 80% of their entire exon lengths. The total coverage for all genes at the nucleotide level was 70.2%. Dot plots are shown in Figure 5, and indicate that Cgaln can be conveniently used to draw a general view of a homology map between mammalian whole genomes.

Identification of human individual differences