Syntenator: Multiple gene order alignments with a gene-specific scoring function
© Rödelsperger and Dieterich; licensee BioMed Central Ltd. 2008
Received: 20 June 2008
Accepted: 06 November 2008
Published: 06 November 2008
Identification of homologous regions or conserved syntenies across genomes is one crucial step in comparative genomics. This task is usually performed by genome alignment softwares like WABA or blastz. In case of conserved syntenies, such regions are defined as conserved gene orders. On the gene order level, homologous regions can even be found between distantly related genomes, which do not align on the nucleotide sequence level.
We present a novel approach to identify regions of conserved synteny across multiple genomes. Syntenator represents genomes and alignments thereof as partial order graphs (POGs). These POGs are aligned by a dynamic programming approach employing a gene-specific scoring function. The scoring function reflects the level of protein sequence similarity for each possible gene pair. Our method consistently defines larger homologous regions in pairwise gene order alignments than nucleotide-level comparisons. Our method is superior to methods that work on predefined homology gene sets (as implemented in Blockfinder). Syntenator successfully reproduces 80% of the EnsEMBL man-mouse conserved syntenic blocks. The full potential of our method becomes visible by comparing remotely related genomes and multiple genomes. Gene order alignments potentially resolve up to 75% of the EnsEMBL 1:many orthology relations and 27% of the many:many orthology relations.
We propose Syntenator as a software solution to reliably infer conserved syntenies among distantly related genomes. The software is available from http://www2.tuebingen.mpg.de/abt4/plone.
Whole genome sequencing has boosted our knowledge database on genome architectures. Identification of conserved genomic regions across species borders has drawn much attention to the field of comparative genomics [1, 2]. The identification of homologous regions between genomes supports genome annotation, function prediction and the study of evolutionary relationships between species. Depending on the level of divergence, homologous regions are usually defined by conserved orders of local genomic alignments , orthologous exons  or genes .
Conservation of gene order across multiple species is usually referred to as 'conserved synteny' or 'collinearity'. In the context of genome evolution, collinear blocks could be used to measure evolutionary distances between genomes in terms of genome rearrangement distances (GRD). The order of all collinear blocks in a genome can be represented as a sequence of signed integers, the GRD denotes the number of rearrangements to transform one such sequence into another .
The standard approach to reconstructing blocks of 'conserved synteny' is to first define a homolog assignment of gene copies. Subsequently, maximal blocks of collinearity are determined on the given homolog assignment and genomic gene orders in the compared genomes.
Traditionally, orthologs were defined by best-reciprocal BLASTP hits (BRH). For example, COGs (Cluster of Orthologous Groups, ) are built from cliques of size 3 in the graph of mutual best cross-species BLAST hits. These seed clusters are subsequently merged into bigger clusters provided that one side is shared between them. Other approaches (e.g.  or ) improve on this approach as they also take gene duplication and gene loss events into account.
The existence of gene families complicates homolog assignment based on protein sequence similarity. The genomic context of a gene copy might provide additional information as to the gene's evolutionary history. Gene copies that are surrounded by the same genes in different genomes are more likely to be true ancestral copies. Consequently, homolog assignment and conservation of gene orders are interlinked and should be jointly studied.
Boyer et al.  present a generic approach to merge information from two or more primary graphs. They explicitly discuss the problem of finding contiguous genes with conserved order across multiple genomes. Gene tuples (one gene per genome) are initially built from a set of orthology relations (protein sequence similarity and alignment coverage cutoff) and enter a multigraph as vertices. These vertices are connected by edge sets, which are defined by the gene order in the respective genomes. Subsequently, common connected components are searched that constitute blocks of conserved gene orders. The worst-case time complexity of the proposed algorithm for finding common connected components is O(n(e·n + m)) where n is the total number of nodes in the multigraph, e is the number of primary graphs and m is the total number of edges in the multigraph. Boyer et al.  noticed that this procedure could be too stringent and allow the insertion of additional edges in the primary graphs. We have re-implemented this method in our Blockfinder algorithm (Additional file 3).
Conceptually more advanced approaches consider all genes of the compared genomes simultaneously. In a partial order alignment approach, a score function is used to integrate protein sequence similarity over the genomic context. Previous work on pairwise gene order alignment has been presented by Haas et al.  and Wang et al. . Both methods resort to dynamic programming approaches that are closely related to the Smith-Waterman algorithm  and operate on directed acyclic graphs. Along these lines, we propose the Syntenator algorithm that facilitates multiple gene order alignments with a novel scoring function. In short, Syntenator is a hybrid approach that combines protein sequence similarity and genomic context dynamically.
Partial gene order alignment
These directed acyclic graphs are subsequently called partial order graphs (POGs). The simple instance of a POG represents a chromosome or genomic contig where each node (gene) has an in-degree = out-degree = 1 (except the start and end nodes). Two simple POGs (and DAGs in general, ) can be aligned by using an extension of the Smith-Waterman algorithm. Figure 1A shows the trace back matrix of a pairwise alignment of two simple POGs from species A and B. Several local alignments (above a user-defined threshold) are extracted from the trace back matrix. These alignments are ranked based on their score (see Figure 1B). This step necessitates the definition of a scoring function for genes and we will present one in the implementation section. The set of pairwise local alignments defines a gene-gene mapping or mapping of vertices of the two input graphs (Figure 1C). Three possible scenarios may occur in this simple example: two genes match, two genes do not match or one gene has no counterpart in the other graph (gap case). In the last step (Figure 1D), two POGs are merged to form a new and possibly more complex POG. All vertices of matching genes are fused into single vertices. The remaining vertices are retained as individual nodes. This merging step must yield a directed acyclic graph for the next (multiple) alignment step. This is obviously the case in our simple example but far from trivial for an alignment of two complex POGs. We will now turn to the actual implementation of Syntenator where we will discuss all relevant aspects of POG alignment.
Syntenator combines conservation of gene order and local sequence similarity to deduce gene orthology. Partial order alignments are represented by partial order graphs (POG). We present an implementation that operates on one POG and one simple chain graph, which is a representation of a linearly ordered gene set (e.g. a genome). An extension of the concept to the alignment of two arbitrary POGs will be discussed in detail.
Modifications to the recurrence relation
Each cell S(n, m) of the dynamic programming matrix is maximized over the four possibilities: match, insertion, deletion and starting a new alignment. The main difference to traditional pairwise local alignment are P and Q, the sets of predecessor nodes of n and m in the corresponding POGs. For complex POGs, we have to consider |P| × |Q| alternative candidates in case of a match. The most simple case is |P| = |Q| = 1 if we were to align two genomes. Our implementation operates on one POG and one simple chain. Consequently, we have either |P| = 1 or |Q| = 1. The expressions s(n, m) and Δ denote the match score for two nodes and the gap penalty, respectively.
Gene order alignment
Initially, all pairwise alignments between two POGs (e.g. G1 and G2) are computed in forward and reverse direction. An alignment in the reverse direction requires the reversal of all edges in one of the two POGs. For each comparison (in both directions), we consider all local (sub)optimal alignments above a certain threshold Θ. All alignments are ranked by their scores in descending order. Based on these alignments, we decide which vertices match and should be fused into a common vertex. We greedily assign vertex matches by traversing the ordered list top-down.
Algorithm 1 (see Appendix) shows the adaptations of the algorithm of Lee et al.  to produce a set of all suboptimal alignment paths P. Such a path consists of a tuple (s, L, r) where s denotes the score, L is a list of aligned node pairs and r indicates wether a gene order was aligned in its original or reversed orientation. The score is adjusted by subtracting the initial score s init which is defined as the last minimal score encountered during traceback before the score exceeds the final alignment score or 0 if no such minimum exists. This adjustment is necessary to prevent that alignments inherit scores from previous higher scoring alignments.
Merging genome graphs
In POA, two graphs, G1 and G2, are merged after each round of pairwise alignments. We have already discussed how to identify pairs of vertices (e.g. (v, w) with v ∈ G1 and w ∈ G2) that should be merged between both graphs. We denote this as 1:1 mapping M.
In the merging step, we iterate over all vertices w ∈ G2 and add a copy of w to G1 if w ∉ M. If (v, w) ∈ M we fuse v and w by copying the genes stored at w to v. If a G1-equivalent of the predecessor node of w exists, we connect this G1-equivalent predecessor node of w to v. All connections between nodes that were not fused, but simply added to the graph, are retained in the merged graph.
These particular problems did not arise in the original implementation for protein or EST sequence alignment (e.g. ) where DAGs are aligned in one defined orientation (e.g. N to C terminus for proteins, 5' to 3' end for ESTs) and just one optimal alignment is reported.
To resolve newly introduced cycles in scenario 1 (Figure 2A), we use a topological ordering of G1 and check at all branching points, whether a loop path consisting of new nodes from G2 induces a cycle in the merged graph G3. We have to test if the loop path returns to a node in G1 at an index which is less or greater in terms of the topological order than the index of the branching point from which we started off. If the path is a forward path and the index of the returning point is smaller than the index of the branching point, all edges within the path have to be reversed to keep the graph acyclic. This procedure leads to G4 in Figure 2A. The case for the backward path works analogously. If the newly added loop is part of a greater loop in G1, we have to search in both directions for the endpoints of the old loop to define an order relation on the newly added loop.
The second case (Figure 2B) emerges if local alignments of opposite orientations exist. In the given example, a cycle would be formed between nodes C and D as they are aligned in opposite orientation to A and B. This is circumvented by keeping the edge orientation of one graph (G1) for the reverse alignment. The "dashed" edges are added to preserve the original order relations of G2.
Repetitive regions that may result from duplication events do not introduce cycles into the merged POG since we greedily enforce a 1:1 mapping of gene nodes. Only the best matching repeat copies would be merged.
Our algorithm relies on BLASTP hits as general similarity measure. From the set of all-against-all BLASTP hits, we save a bitscore for each gene pair in a lookup table. In case of alternative transcripts the highest score between any two protein products is saved.
We chose a scoring function that allows us to order alignments according to the number of aligned pairs or to the sum of pairwise similarities in case of equal numbers of pairs.
For each pair of genes (A, B) a symmetric score function is given by Eqn. 2. The individual contributions are shown in Eqn. 3.
We require sbitscore to be ≥ 50. The match score is always < 2: .
This can be interpreted as summing up over the entries of a non-symmetric weighted adjacency matrix of all pairwise homology relationships. A mismatch score is assigned if the two genes under comparison either have no BLAST hit or if they are located on different strands.
nv, wdenotes the number of genes of nodes v and w, denotes the number of species in the graphs of v and w. The term Cv, win the denominator of Eqn. 4 is a scaling factor whose definition depends on the current alignment score. Cv, wis equal to the number of comparisons between either all species in nodes v and w or the number of all species in the graphs of v and w (Eqn. 5). This correction scheme was implemented because weak BLAST hits tend to appear in the set of genes of both vertices more often if the number of compared genes increases. As a consequence pairwise scores tend to be higher than the averaged scores of multiple comparisons. In order to equalize this effect, we replace nv, wby as soon as the alignment score σ exceeds the threshold Θ. This triggers a switch towards a more specific search for alignments containing genes from multiple species.
We applied both approaches to detect conserved syntenies in four mammalian species, namely human (NCBI 36), mouse (NCBI m36), rat (RGSC 3.4) and dog (CanFam 1.0). We computed all pairwise all-against-all BLASTP searches in advance. The BLASTP hit ranks and bitscores are subsequently used by Blockfinder and Syntenator. The number of genes with putative homologs at an E-value cutoff < 0.1 is shown in Additional File 1. Only these genes are considered in whole genome alignments. We contrasted our findings to the EnsEMBL compara database, which reports pairwise conserved synteny relations based on nucleotide alignments.
Application of Syntenator
We used the aforementioned data to construct POGs for all genomes. Classical methods like best reciprocal hits and COGs  select best BLAST hits to assign orthologs. We suggest that in order to maximize conserved synteny, non-best hits should be taken into account. Nevertheless highly abundant protein domains drastically increase the number of BLAST homologies for certain genes [9, 15] but these homologs are unlikely to be true 1:1 orthologs. In order to reduce the amount of data being passed on to Syntenator we apply certain filters: The BLAST similarity relations were filtered to contain only the 5 best hits per query. Hits were further removed if their bitscore dropped below 95% of the best score.
If we chose to include more BLAST similarity relations per gene, more alignments would pass the minimal threshold Θ. That is why, the actual choice of the BLAST similarity relations is a tradeoff between speed and sensitivity. Our filtering step cuts down on the number of candidate alignments that would have to be evaluated.
Syntenator was run on this data set using a linear gap score of - 2.0, a mismatch score of - 3.0 and a threshold of 2.0. These values were motivated by assuming that a complete loss of two genes is less likely than a mismatch between two diverging genes. A threshold of 2.0 requires that a pairwise ungapped alignment consists of at least two gene pairs.
Pairwise genome comparison
We compared the performance of gene order alignment approaches (Blockfinder and Syntenator) to the EnsEMBL compara synteny data set. Herein, Blockfinder utilized three homology data sets, which are all based on EnsEMBL release 46: 1) Ensembl orthologs (1:1, 1:many and many:many), 2) Best reciprocal BLASTP hits (BRH) and 3) 3-best-reciprocal BLASTP hits (BRH3). In the last case, a gene may have up to 3 hits. Generally, a gene node may have up to (g - 1) * n homology relations to other genes, where g is the number of species and n is the number of considered BLASTP hits.
How much conserved synteny information we lose as compared to the "gold" standard as given by the EnsEMBL compara data
How much we improve over simpler methods that define homology relations in advance (e.g. Blockfinder).
However, the full potential of our method unfolds when two remotely related species are aligned. We compared whole genome alignment coverage of Syntenator and UCSC blastz runs on the nucleotide level. Blastz  is a pairwise whole genome alignment method, which produces local nucleotide sequence alignments.
Multiple genome comparison
Syntenator multiple gene order alignments
Alignment order: MRHD
Number of blocks
Genes in alignments
Block size in Mb
Alignment order: HDMR
Number of blocks
Genes in alignments
Block size in Mb
In the last round of the multiple alignment, either the POG of human, dog and mouse is aligned to the rat genome or the POG of mouse, rat and human is aligned to the dog genome. It is apparent that the final multiple gene order alignment is sensitive to the order of alignment steps. After this last alignment round, Syntenator reports only genes that have been aligned with the genome that was added last (either dog or rat). This is also reflected in Table 1 where the genome that was added last shows the highest percentage of aligned genes. Please note that multiple genome alignments do not necessarily contain genes from all species. In total, there are 11,164 and 11,309 alignment nodes in the MRHD and HDMR POGs that consist of 4-tuples (nodes with one gene from each species). This is close to the lowest number of genes from a single species in the two multiple gene order alignments (see Table 1). Future work will address alternative scoring schemes as well as a more rigorous assessment of the impact of alignment orders.
Comparison of orthology prediction
Another potential application of Syntenator is its use to assign gene homology relations. To this end, we compared orthology predictions of Syntenator and the EnsEMBL system. Whole genome alignments with Syntenator were performed with an alignment score threshold of 1.0 so that a prediction is made for each gene with a homolog. The other parameters were set to -3 for a mismatch and -2 for a gap. Firstly, we checked how many EnsEMBL 1:1 orthologs are contained in the Syntenator gene pairs. Except for the man/mouse comparison (94%), ~97% of all EnsEMBL 1:1 orthologs are also predicted by Syntenator (see Additional File 2). Many instances of 1:many or many:many EnsEMBL orthologs could be formally resolved to 1:1 orthologs (see Additional File 1). Resolving 1:many and many:many relations bears the question of what a true ortholog is. We argue for defining orthology on protein sequence similarity and gene order. This argument has been previously made in the context of gene function prediction . Our Syntenator framework accomplishes this task. However, a good test set is not available to our knowledge and simulating whole genome evolution is beyond the scope of this manuscript.
We have established Syntenator as a new method to identify regions of gene order conservation over multiple genomes. Furthermore, we propose that our method could be used to resolve gene homologies. Instead of defining an orthologous group from sequence similarity alone, our method chooses the ortholog from a set of candidate genes according to available synteny information. This observation is necessary as relying on best reciprocal hits exclusively does not guarantee to find the 'true' ortholog. This circumstance might be explained by the weakened selective pressure on duplicated genes .
Blockfinder chooses orthologs from a set of candidate orthologous genes by maximizing collinearity across all species. The initial clique graph does not capture all existing BLAST homologies. Genes outside of cliques are excluded from the subsequent analysis. In general, this is a disadvantage but turns into an advantage when genomes with poor gene annotations are used.
Syntenator integrates all gene positions and complete BLAST data into the computation of collinear blocks. Herein, synteny information is used as the first criterion to define orthology, although substantial local sequence similarity as expressed by BLAST scores is still required.
In summary, our work extends existing methods for orthology prediction and provides new tools to compare local and global genome architectures of multiple species, especially for genomes that do not align on the nucleotide level.
Availability and requirements
Project name: Syntenator
Project home page: http://www2.tuebingen.mpg.de/abt4/plone/projects/syntenator
Operating system(s): Platform independent
Programming language: Java
Other requirements: Java 1.4.2 or higher
License: freely available to academia
Any restrictions to use by non-academics: license is needed
Algorithm 1: Computing a set of suboptimal gene order alignments
N, M are the number of nodes in both graphs. A is the dynamic programming matrix and T is the traceback matrix. Cells of T contain the index tuple of the predecessor cell pointing to any cell in the computed area. The indices i and j iterate over the topological orders of both graphs (line 5,6).
1: for i ← 0 to N do // initialize matrices
2: A(0,i) ← 0, T(0,i) ← (0, 0)
3: for j ← 0 to M do
4: A(j, 0) ← 0, T(j, 0) ← (0, 0)
5: for i ← 1 to N do // dynamic programming
6: for j ← 1 to M do
7: (A(j, i),T(j, i)) ← Score(j, i, A, T)
8: (p j ,p i ) ← T (j, i)
9: if A(j, i) > A(p j ,p i ) then
10: L ← L ∪ (A(j, i),j, i) // store each cell with increasing score
Score(j, i, A, T) fills cells A(j, i) and T(j, i) according to Section "Score function" and Eqn. 1 (line 7).
Subsequently the scores A(j, i) and A(p j ,p i ) are compared and cells with increasing score are stored as candidates in L (line 9,10). The candidate alignments in L are processed by decreasing score. An alignment path p is stored, if the difference s - s init exceeds the threshold Θ (line 14–16).
11: for k ← to |L| do
12: (s, j, i) ← L(k)
13: s init ← InitialScore(j, i, A, T)
14: if s - s init >Θ then
15: p ← Traceback(j, i, s init ,A, T)
16: P ← P ∪ p // update the sorted set of paths
Both authors acknowledge funding of the Max Planck Society.
- Miller W, Makova KD, Nekrutenko A, Hardison RC: Comparative genomics. Annu Rev Genomics Hum Genet. 2004, 5: 15-56. 10.1146/annurev.genom.5.061903.180057PubMedView ArticleGoogle Scholar
- Murphy WJ, Pevzner PA, O'Brien SJ: Mammalian phylogenomics comes of age. Trends Genet. 2004, 20 (12): 631-639. 10.1016/j.tig.2004.09.005PubMedView ArticleGoogle Scholar
- Stein LD, Bao Z, Blasiar D, Blumenthal T: MRB: The genome sequence of Caenorhabditis briggsae: a platform for comparative genomics. PLoS Biol. 2003, 1 (2): E45- 10.1371/journal.pbio.0000045PubMedPubMed CentralView ArticleGoogle Scholar
- Brudno M, Poliakov A, Salamov A, Cooper GM, Sidow A, Rubin EM, Solovyev V, Batzoglou S, Dubchak I: Automated whole-genome multiple alignment of rat, mouse, and human. Genome Res. 2004, 14 (4): 685-692. 10.1101/gr.2067704PubMedPubMed CentralView ArticleGoogle Scholar
- Haas BJ, Delcher AL, Wortman JR, Salzberg SL: DAGchainer: a tool for mining segmental genome duplications and synteny. Bioinformatics. 2004, 20 (18): 3643-3646. 10.1093/bioinformatics/bth397PubMedView ArticleGoogle Scholar
- Pevzner P, Tesler G: Genome rearrangements in mammalian evolution: lessons from human and mouse genomes. Genome Res. 2003, 13: 37-45. 10.1101/gr.757503PubMedPubMed CentralView ArticleGoogle Scholar
- Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR: BK: The COG database: an updated version includes eukaryotes. BMC Bioinformatics. 2003, 4: 41- 10.1186/1471-2105-4-41PubMedPubMed CentralView ArticleGoogle Scholar
- Alexeyenko A, Tamas I, Liu G, Sonnhammer ELL: Automatic clustering of orthologs and inparalogs shared by multiple proteomes. Bioinformatics. 2006, 22 (14): e9-15. 10.1093/bioinformatics/btl213PubMedView ArticleGoogle Scholar
- Goodstadt L, Ponting CP: Phylogenetic reconstruction of orthology, paralogy, and conserved synteny for dog and human. PLoS Comput Biol. 2006, 2 (9): e133- 10.1371/journal.pcbi.0020133PubMedPubMed CentralView ArticleGoogle Scholar
- Boyer F, Morgat A, Labarre L, Pothier J, Viari A: Syntons, metabolons and interactons: an exact graph-theoretical approach for exploring neighbourhood between genomic and functional data. Bioinformatics. 2005, 21 (23): 4209-4215. 10.1093/bioinformatics/bti711PubMedView ArticleGoogle Scholar
- Wang X, Shi X, Li Z, Zhu Q, Kong L, Tang W, Ge S, Luo J: Statistical inference of chromosomal homology based on gene colinearity and applications to Arabidopsis and rice. BMC Bioinformatics. 2006, 7: 447-10.1186/1471-2105-7-447PubMedPubMed CentralView ArticleGoogle Scholar
- Smith TF, Waterman MS: Identification of common molecular subsequences. J Mol Biol. 1981, 147: 195-197.10.1016/0022-2836(81)90087-5PubMedView ArticleGoogle Scholar
- Lee C, Grasso C, Sharlow MF: Multiple sequence alignment using partial order graphs. Bioinformatics. 2002, 18 (3): 452-464. 10.1093/bioinformatics/18.3.452PubMedView ArticleGoogle Scholar
- Grasso C, Lee C: Combining partial order alignment and progressive multiple sequence alignment increases alignment speed and scalability to very large alignment problems. Bioinformatics. 2004, 20 (10): 1546-1556. 10.1093/bioinformatics/bth126PubMedView ArticleGoogle Scholar
- Tatusov RL, Koonin EV, Lipman DJ: A genomic perspective on protein families. Science. 1997, 278 (5338): 631-637.10.1126/science.278.5338.631PubMedView ArticleGoogle Scholar
- Schwartz S, Kent WJ, Smit A, Zhang Z, Baertsch R, Hardison RC, Haussler D, Miller W: Human-mouse alignments with BLASTZ. Genome Res. 2003, 13: 103-107. 10.1101/gr.809403PubMedPubMed CentralView ArticleGoogle Scholar
- Notebaart RA, Huynen MA, Teusink B, Siezen RJ, Snel B: Correlation between sequence conservation and the genomic context after gene duplication. Nucleic Acids Res. 2005, 33 (19): 6164-6171. 10.1093/nar/gki913PubMedPubMed CentralView ArticleGoogle Scholar
- Kondrashov FA, Rogozin IB, Wolf YI, Koonin EV: Selection in the evolution of gene duplications. Genome Biol. 2002, 3 (2): RESEARCH0008- 10.1186/gb-2002-3-2-research0008PubMedPubMed CentralView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.