Statistical inference of chromosomal homology based on gene colinearity and applications to Arabidopsis and rice

Background The identification of chromosomal homology will shed light on such mysteries of genome evolution as DNA duplication, rearrangement and loss. Several approaches have been developed to detect chromosomal homology based on gene synteny or colinearity. However, the previously reported implementations lack statistical inferences which are essential to reveal actual homologies. Results In this study, we present a statistical approach to detect homologous chromosomal segments based on gene colinearity. We implement this approach in a software package ColinearScan to detect putative colinear regions using a dynamic programming algorithm. Statistical models are proposed to estimate proper parameter values and evaluate the significance of putative homologous regions. Statistical inference, high computational efficiency and flexibility of input data type are three key features of our approach. Conclusion We apply ColinearScan to the Arabidopsis and rice genomes to detect duplicated regions within each species and homologous fragments between these two species. We find many more homologous chromosomal segments in the rice genome than previously reported. We also find many small colinear segments between rice and Arabidopsis genomes.


Background
Exploration of homology between chromosomes facilitates the understanding of the structure, function and evolution of genomes. Extensive synteny and colinearity have been detected between chromosomes in different species of cereals [1], mammals [2] and yeasts [3] providing a deep insight into the evolution of ancient chromosomes.
Special care should be taken to reveal chromosomal homology due to numerous genomic changes such as chromosomal rearrangements, gene inversions and gene loss [13][14][15]. Many approaches have been developed to identify chromosomal homologues [16] based on genetic maps [17], sequence alignment [18,19], gene synteny [10] and gene colinearity [20][21][22][23]. By detecting the density and order of homologous gene pairs between chromosomes, colinearity approach can reveal reliable homologous regions and requires less computational resources. This approach also enables us to develop reasonable statistical tests to evaluate the significance of predicted homologous regions.
The typical implementations of the colinearity strategy are ADHoRe [20], FISH [24] and DiagHunter [25]. The implementations of these approaches have limitations in some aspects, though they have been widely adopted. Firstly, the gap size between neighboring homologous genes which is essential to define and detect true colinearity needs further evaluation [12,[20][21][22][23]26]. Secondly, statistical tests to assess predicted homologous regions are mainly based on a prerequisite of balanced gene loss rates between homologous regions. Finally, compositional and structural differences, especially gene density and repetition in genome-wide and local chromosomal regions, have not been fully addressed.
Here we describe a new colinearity approach characterized by improved statistical inference, flexibility and computational efficiency. Firstly, the selection of parameter values is theoretically explored, especially that of the gap length between neighboring genes. Secondly, the statistical test has been substantially strengthened with a mathematical deduction to evaluate the significance of the predicted homologous regions. Finally, the compositional and structural heterogeneity of chromosomes has been considered.
Using a dynamic programming algorithm, we developed ColinearScan and scanned the Arabidopsis and rice genomes to detect duplicated regions in each species and homologous chromosomal regions between these two species. We found 75.0% of Arabidopsis genes and 76.2% of rice genes were in duplicated regions. Moreover, we identified homologous fragments between these two species, in 32.9% of Arabidopsis and 16.8% of rice. Nearly all homologous segments were shorter than 0.6 Mb, indicating massive chromosomal rearrangements after the monocot-dicot divergence [27].

Algorithm
Gene homology matrix The first step in our colinearity approach is the construction of the gene homology matrix. To find homologous gene pairs between two chromosomes denoted as A and B, protein sequences encoded by genetically or physically positioned genes are used to perform an all-against-all BLAST search [28]. A gene homology matrix (GHM, denoted as H) is then constructed using the homology information from BLAST results. Chromosome A and B, represented by the positioned genes are arranged along H horizontally and vertically (Fig. 1A). A cell of H is filled with "1" if the corresponding genes on chromosome A and chromosome B are homologous, otherwise with "0". Tandem and other repetitive genes are widely distributed in chromosomes showing many "1"s in horizontal or vertical straight lines in the dot matrix map (Fig. 2) and causing problems in revealing true homology. Therefore, we used a general approach, masking the genes appearing more than 10 times in both chromosomes.

Dynamic programming algorithm
To reveal the homologous genes in colinearity between two chromosomes, we implemented a dynamic programming approach based on the well-known Smith-Waterman algorithm [29]. Using this approach, we can discover the longest putative sister regions represented by several proximal points of colinear homologous gene pairs in nearly diagonal orientations. The points may not be in close proximity due to large-scale gene loss, insertion and translocation. The extent of the proximity of the points is essential to reveal and evaluate the colinear sister segments. Lines forming by the points corresponding to the true colinear segments are either nearly parallel to the main diagonal line or the anti-diagonal line due to DNA segmental inversion. Homologous genes in colinear segments should all have the same or inverse transcriptional directions if no single gene inversions occur. We scan H in two directions, starting from the upper-left and upperright. Here, we describe the procedure starting from the upper-left, which also applies in the other direction. Transcriptional orientations of the genes are recorded but not used when performing the colinearity search.
To reveal the colinearity represented by the proximal points in H, we introduce a parameter mg (the maximum gap length) between two neighboring points. Then we define another matrix S (the scoring matrix) with the same size as H (Fig. 1B). A cell in matrix S represents the extension of a colinearity path, i.e., the value of each cell is the number of collinear gene pairs in the path accumulated from its starting point. The path extends and the value of the cell increases by 1 if there is a "1" in lower-A modified Smith-Waterman algorithm to locate colinearity Figure 1 A modified Smith-Waterman algorithm to locate colinearity. (A) A simplified gene homology matrix (GHM, denoted as H). Genes A1, A2, ..., A18 on chromosome A are arranged horizontally, and genes B1, B2, ..., B14 on chromosome B are arranged vertically. Each cell of the matrix is filled with "1" or "0" based on the homology information from BLASTP search, e.g., gene A1 and gene B1 are homologous, and gene A2 and B2 are non-homologous. (B) A modified dynamic programming procedure. A scoring matrix S is constructed recursively based on H, with mg set to 2 genes apart. The distance criterion demands that neighboring genes in colinearity are no more than 2 genes apart. Pointers are shown by dark or grey arrow lines. Two collinear paths containing 9 and 5 genes are shown by dark arrow lines reflecting the same colinear relationship between the corresponding chromosomal regions.
right neighborhood, and both vertical and horizontal distances are less than mg.
Initially, S is identical to H. We rebuild the matrix S recursively using a dynamic programming procedure: where S(i, j) is the score computed, H(i, j) is the homology information, pre(i, j) is the cell leading to the maximum score at the cell p(i, j), and a pointer (denoted by dark or gray arrow lines in Fig. 1B) is created from the cell p(i, j) to pre(i, j), dist(p(i, j), p(a, b)) is the distance between the cells p(i, j) and p(a, b). Eventually, the maximum score in S corresponds to the longest putative collinear segments. The longest colinearity path formed by dots of the homologous gene pairs is revealed by a trace-back procedure according to the pointers created. After the homologous genes in putative colinearity are recorded, we mask these putative colinear segments by setting H(i, j) to 0, rebuild the matrix S, and scan for other putative colinear segments till no sister regions containing more colinear genes than a threshold r could be found.

Maximum gap length
In colinearity methods, mg is the most important parameter which determines the length, quality and extensiveness of the predicated colinearity. The frequency of gene deletion in duplicated chromosomal segments is high and only a small fraction of homologous genes remains in colinearity. A small value of mg will result in finding many small colinear segments, and increase the difficulty of interpreting possible evolutionary events. On the other hand, a large mg value will surely result in high false positives. In fact, mg is dependent on the density of homologous gene pairs between the chromosomes. When the two chromosomes A and B are from the same species, homologous genes between them are mainly distributed at a similar density. On the other hand, when we compare two chromosomes from different species, the density of homologous genes may more divergent. Therefore we adopt two parameters to define the maximum distance between the neighboring dots, the maximum gap between genes in chromosome A (mgA) and the maximum gap between genes in chromosome B (mgB). When the chromosomes are from the same species, we set mgA = mgB.
Under the assumption that homologous genes are uniformly distributed in chromosomes, we explore the possibility of finding sister segments containing equal to or more than r genes by chance. However, this uniform distribution assumption is not very strict since we only need reasonable rather than optimal values of mg. Suppose the length of chromosome A and B are lenA and lenB, and the number of homologous gene pairs between two chromo-somes is pnum. The location of a gene is a random variable with a probability density , the joint probability density of the locations of r genes on chromosome A is , and the joint probability density of the locations on chromosome A and B of r homologous pairs is . Therefore, the probability p that r homologous gene pairs are in colinearity by chance can be evaluated by where the multiple integral field D is and x 1 ,x 2 ,...,x r ;y 1 ,y 2 ,...,y r are the positions of the genes on the chromosomes, is the number of permutations of homologous genes: . When mgA <<lenA and mgB <<lenB, the integral in the above formula can be approximated by Then we can estimate p by If we set the significance level of colinearity to α, then we have When the two chromosomes are from the same species, we assume mgA = mgB and denote them as mg, and the value of mg can be evaluated by When the two chromosomes are from different species, the gene densities are often different and thus we adopt , pre(i-1,j-1)) mg and S(i-1,j-1) S(pre(i-1,

Colinearity shadow
Many genes have homologues in their neighborhood as indicated by the distance distribution of homologous genes (Fig. 3). Neighboring homologues may result in many colinear segments parallel to each other in sister regions within the same chromosome, or between the different chromosomes ( Fig. 2A, 2C). If for convenience we call the longer ones the 'true' colinearity, then the shorter ones are the 'shadows'. The colinearity shadows may also reflect colinearity, but they do not provide further information and greatly increase the computational time. Thus after a putative colinearity in certain sister regions is found, we mask the neighboring homologous pairs within the corresponding the entire rectangular regions to avoid possible colinearity shadows.
When candidate colinearity is found in a specific region of GHM from upper-left to lower-right, colinear shadows perpendicular to this path might also be found in the same region scanning from upper-right to lower-left. These shadows may also reflect actual homology between two sister chromosomal segments, and they may occur when gene rearrangements are frequent in specific regions. A large amount of perpendicular shadows may affect efficiency of the algorithm. However, they do not occur very often so we do not mask these shadows.

Statistical test
Many genes are in multi-gene families and repetitive genes are extensively found in plant genomes. Putative colinear regions might be a reflection of extensive occurrences of single gene duplication or translocation. Therefore, it is critical to correctly detect the colinearity pattern by evaluating the significance of the putative colinear regions gen- erated by the above approach, and to develop an effective statistical method to assess whether the putative colinearity represents true homology or is produced by chance. Moreover, the distribution of homologous gene pairs is far from uniform. We use a statistical test considering divergent distribution of homologous gene pairs in different regions, rather than assuming a uniform distribution throughout the genome.
Given that a DNA segment resides in chromosome B, a corresponding colinear segment in chromosome A is generated due to many independent single-gene duplication events. Under the assumption of a uniform distribution of collinearly paired genes in the local chromosomal regions, we obtain the probability epvA to display the significance of the homology, to measure the possibility of the colinearity appearing by chance in the sister regions. If it is below a threshold, we take the putative colinearity as significant.
However, if we apply the above test directly to the putative different homologous regions, we cannot distinguish between patterns with similar numbers of homologous gene pairs in segments of different size. For example, the putative colinear regions found in small sister segments (Fig. 4B) more likely indicate true segmental homology than that in large sister segments (Fig. 4A). Dense coline- arity possibly generated by recent duplications should be evaluated as significant.

Different distribution pattern of homologous genes in sister regions
Rather than using pssA and pssB, we try to estimate epv on expected scales of the sister regions (essA and essB) determined by the numbers of colinear genes and all homologous genes. Assuming that paired genes are uniformly distributed throughout the chromosomes, we estimate the expected sizes for each pair of sister segments as follows where cnum is the number of colinearly paired genes in the putative sister segments, pnum is the number of all homologous gene pairs between chromosome A and chromosome B, λ is a coefficient to normalize the size of different putative colinear segments. For each pair of putative sister regions we calculate a preliminary coefficient number λ i by .
where pssA and pssB are the original sizes of the putative ith sister segments. These preliminary coefficients are averaged throughout the putative homologous sister pairs to obtain the estimate of λ. Finally, we redefine and by substituting the original size pssA and pssB with the expected size essA and essB. Thus, we assign different significance to the two patterns in Fig. 4, the sister regions with denser collinear genes have a smaller epv value.

Assessment of mg estimation
To test the applicability of the criteria defining mg, we performed a computational simulation test on rice chromosome 1 and Arabidopsis chromosome 1. Second, we shuffled the positions of the homologous gene pairs on both chromosomes and scanned the longest homologous segments occurring by chance, then checked whether it had 4 or more collinear genes. We repeated this process 1000 times and found collinear segment with 4 or more collinear genes 9 times.
We calculated mgA and mgB between every pair of chromosomes of rice and Arabidopsis, and applied the largest values (mgA = 160 Kb and mgB = 154 Kb) to all chromosome pairs. The candidate homologous segments can be verified by a statistical test. To check how it affects the scanning process between rice chromosome 1 and Arabidopsis chromosome 1 when using larger mgA and mgB, we found colinear segments with 4 or more collinear gene pairs in 33, out of 1000, simulations.

Application to rice and Arabidopsis
We explored the colinear segments within and between the chromosomes in Arabidopsis and rice. The genomic sequences of Arabidopsis thaliana were from GenBank (Accession NC_003070, NC_003071, NC_003074, NC_003075, NC_003076) [30]. Oryza sativa L. ssp. indica genomic sequences were downloaded from RiceGD [31] and the rice genes were predicted using the software BGF [32] from the Beijing Genome Institute. By performing all-against-all BLASTP, we revealed homologous gene pairs within and between Arabidopsis and rice (BLASTP score > 100).
For Arabidopsis, we searched for the duplication regions with the parameter value r = 4, mg = 116 Kb (~25 intervening genes) and λ = 2.6. We used the maximal mg values estimated between each pairs of Arabidopsis chromosomes. At the significant level of epv <= 0.01, we discovered 203 duplicated sister segments out of 350 candidates, among them 3 were possible perpendicular shadows (Supplementary table 1 [see Additional file 1]). About 75.0% of the genome is in duplicated segments and 22.4%, 1.8% of the genes are in segments with a multiplication level > 2 and > 4, respectively ( Table 1). The detected coverage of duplicated regions is a little more than that (71%) found by Blanc et al. [18], less than that (89%) reported by Bowers et al. [33]. The longest sister segments contain 106 colinear genes and extend more than 1.8 Mb in chromosome 2 and 1.55 Mb in chromosome 3 (  6 Mb in length and contain ~14 colinear genes, but most segments are much shorter, indicating extensive independent chromosomal rearrangements and gene loss or gain in each genome. The sister copy in rice is always 1-4 times longer than that in Arabidopsis, implying a possible chromosomal expansion in the rice genome (Table 4).  Table 4 and 5) reflect the chromosome and the gene order, e.g. '2_3518' represents the 3518th gene on chromosome 2.

Discussion
Identification of the duplicated segments, especially their distribution pattern in a genome, is essential for further inference on when and how the duplication or species divergence occurred, and whether or not recurrent duplication events happened. The selection of parameter values, in particular the maximum gap length between the neighboring genes, is critical to detect chromosomal homology. However, the selection of maximum gap length in previous reported studies was mainly empirical, which might fail to detect authentic duplicated segments [20,22,23]. Many fewer and shorter duplicated segments are discovered when a smaller gap length is adopted, such as in the case of rice [20], whereas more and longer duplicated regions can be found if a larger gap length is adopted. Moreover, different gap length should be used in different genomes such as Arabidopsis and rice, since the density of colinear genes varies due to DNA loss and insertion. By considering the difference in gene density, especially the density of homologous genes in different genomes, we determined the maximum gap length based on statistical analysis. For example, when the duplicated regions in Arabidopsis and rice are detected, the maximum gap lengths were estimated to be 116 Kb and 334 Kb, respectively.
The input data of our approach can be any type of genetic markers such as sequences, genetic markers. Various measurements can be used to represent the distance between markers, such as physical or genetic distances, or gene numbers. In most previous studies, the significance of the predicted colinear regions was evaluated by a permutation test, which is rather time-consuming [20,23]. We estimate the significance of the predicted colinear segments through statistical inference. The statistical inference has the advantage over the permutation test in terms of computational efficiency. It takes only 2 minutes to calculate the epv to evaluate their significance on a personal computer (AMD AthlonXP 2000+, 512 MB RAM) while running a permutation test takes several hours on the same machine. The massive gene duplications and translocations in its proximal regions will lead to many colinearity shadows, decreasing the computational efficiency. We include a neighborhood masking procedure in Colin-earScan to remove colinearity shadows in our algorithm, which dramatically improves the efficiency of detecting duplicated segments in the rice genome.
ADHoRe adopts linear regression analysis to infer duplicated chromosomal segments [20,21]. The underlying assumption is that gene loss rates have been balanced between sister segments, resulting in a straight line in the dot map. The colinear homologues in a chromosomal segment might be interspersed by individual genes that have no homologues at the corresponding position in its sister segment. At the very beginning of divergence of the sister segments, there should be one-to-one gene homology. Thereafter, massive gene deletions, translocations and chromosomal rearrangements occur and the initial pattern eventually becomes obscured [25]. The homologues with the conservative orders would appear in a straight line in the dot map if gene deletion or insertion had been balanced in different regions of the sister segments, otherwise in a curvy line. Wang et al. [15] explore the gene loss rates in the sister segments in rice and find that nearly straight lines are obtained for some sister segments, e.g., in chromosomes 11 and 12, and in chromosomes 2 and 4. However, curvy lines are also found for some sister segments, e.g., in chromosomes 1 and 5, and in chromosomes 8 and 9. A linearity assumption might fail to detect true duplicated segments. In FISH, Calabrese et al. [24] also adopt a colinearity strategy and develop a different statistical approach to evaluate the extension of collinear points, referred as clump in GHM. However, the value of key parameter p, reflecting the probability that a point occurs in the neighborhood of the former point, is artificially defined, and the maximal gap is deduced from p in their approach. DiagHunter [25] adopts a colinearity method similar to our approach, and the maximal length of the path is predefined. The program stops extending the current path until it reaches the maximal length threshold, or other neighboring points cannot be found.
Polyploidy has been supposed to be prevalent in plants.
Recently, genome-wide studies further suggest the ubiquity of polyploidy, even in genomes which have not been considered to undergo genomic duplication [35]. The small genome of Arabidopsis has been reported to have undergone at least one round of duplication by different groups [12,18]. Here using a different method, we discover that 75.0% of the Arabidopsis genome sequences are in duplicated regions and a significant portion of sequences have multiple copies. The previous studies in the rice genome have been focused on the large obvious duplicated segments, produced by the relatively recent duplication events [15,36]. Here, we detect 76.2% of rice sequences in duplicated regions, and 42.9% have multiple copies.
The possibility of constructing the monocot-dicot comparative genetic map has been discussed [37] based on the comparison of Arabidopsis and rice sequences [38,39]. However, a comprehensive detection of homologous regions between these two genomes has not been available. Based on gene colinearity, we detected homologous regions between Arabidopsis and rice, accounting for 32.9% and 16.9% of each genome. All homologous segments were shorter than 0.6 Mb in length, indicating the massive genome rearrangements in both genomes after the monocot-dicot divergence. Though the short homologous segments make it difficult to construct the comparative genetic map between monocot and dicot, the homologues in colinearity found in this study may provide clues for further work in comparative genomics.

Conclusion
We develop an algorithm to detect homologous chromosomal segments with conserved gene order, and we propose a statistical approach to estimate parameters and evaluate the significance of potential homology. We apply this approach to rice and Arabidopsis with high efficiency to detect potential colinear regions and evaluate their significance. We find many more homologous chromosomal segments in rice genomes than previously reported, which consolidated the inference that a polyploidy had occurred in the common ancestor of grasses. We also find many small colinear segments between rice and Arabidopsis genomes, providing clues to the evolutionary history of monocots and dicots.