Aligning coding sequences with frameshift extension penalties

Background Frameshift translation is an important phenomenon that contributes to the appearance of novel coding DNA sequences (CDS) and functions in gene evolution, by allowing alternative amino acid translations of gene coding regions. Frameshift translations can be identified by aligning two CDS, from a same gene or from homologous genes, while accounting for their codon structure. Two main classes of algorithms have been proposed to solve the problem of aligning CDS, either by amino acid sequence alignment back-translation, or by simultaneously accounting for the nucleotide and amino acid levels. The former does not allow to account for frameshift translations and up to now, the latter exclusively accounts for frameshift translation initiation, not considering the length of the translation disruption caused by a frameshift. Results We introduce a new scoring scheme with an algorithm for the pairwise alignment of CDS accounting for frameshift translation initiation and length, while simultaneously considering nucleotide and amino acid sequences. The main specificity of the scoring scheme is the introduction of a penalty cost accounting for frameshift extension length to compute an adequate similarity score for a CDS alignment. The second specificity of the model is that the search space of the problem solved is the set of all feasible alignments between two CDS. Previous approaches have considered restricted search space or additional constraints on the decomposition of an alignment into length-3 sub-alignments. The algorithm described in this paper has the same asymptotic time complexity as the classical Needleman–Wunsch algorithm. Conclusions We compare the method to other CDS alignment methods based on an application to the comparison of pairs of CDS from homologous human, mouse and cow genes of ten mammalian gene families from the Ensembl-Compara database. The results show that our method is particularly robust to parameter changes as compared to existing methods. It also appears to be a good compromise, performing well both in the presence and absence of frameshift translations. An implementation of the method is available at https://github.com/UdeS-CoBIUS/FsePSA. Electronic supplementary material The online version of this article (doi:10.1186/s13015-017-0101-4) contains supplementary material, which is available to authorized users.


Background
Biological sequence alignment is a cornerstone of bioinformatics and is widely used in such fields as phylogenetic reconstruction, gene finding, genome assembly. The accuracy of the sequence alignments and similarity measures are directly related to the accuracy of subsequent analysis. CDS alignment methods have many important applications for gene tree and protein tree reconstruction. In fact, they are useful to cluster homologous CDS into groups of orthologous splicing isoforms [1,2] and combine partial trees on orthology groups into a complete protein tree for a gene family [3,4]. Aligning and measuring the similarity between homologous CDS requires to account for frameshift (FS) translations that cannot be detected at the amino acid (AA) level, but lead to a high similarity at the nucleotide level between functionnaly different sub-sequences.
FS translation consists in alternative AA translations of a coding region of DNA using different translation frames [5]. It is an important phenomenon resulting from different scenarios such as, insertion or deletion of a

Open Access
Algorithms for Molecular Biology *Correspondence: safa.jammali@USherbrooke.ca 1 Département d'informatique, Faculté des Sciences, Université de Sherbrooke, Sherbrooke, QC J1K2R1, Canada Full list of author information is available at the end of the article nucleotide sequence whose length is not a multiple of 3 in a CDS through alternative splicing [6,7] or evolutionary genomic indels [8,9], programmed ribosomal frameshifting [10], or sequencing errors [11]. Recent studies have reported the role of FS translations in the appearance of novel CDS and functions in gene evolution [6,12]. FS translation has also been found to be linked to several diseases such as the Crohn's disease [13]. The computational detection of FS translations requires the alignment of CDS while accounting for their codon structure. A classical approach for aligning two CDS used in most alignment tools [14,15] consists in a three-step method, where the CDS are first translated into AA sequences using their actual coding frame, then AA sequences are aligned, and finally the AA alignment is back-translated to a CDS alignment. This approach does not account for alternative AA translations between two CDS and it leads to incorrect alignment of the coding regions subject to FS translation. The opposite problem of aligning protein sequences while recovering their hypothetical nucleotide CDS sequences and accounting for FS translation was also studied in several papers [16,17].
Here, we consider the problem of aligning two CDS while accounting for FS translation, by simultaneously accounting for their nucleotide and AA sequences. The problem has recently regained attention due to the increasing evidence for alternative protein production through FS translation by eukaryotic gene families [18,19].
The problem was first addressed by Hein et al. [20,21] who proposed a DNA/protein model such that the score of an alignment between two CDS of length n and m is a combination of its score at the nucleotide level and its score at the AA level. They described a O(n 2 m 2 ) algorithm in [20], later improved to a O(nm) algorithm in [21] for computing an optimal score alignment, under the constraint that the search space of the problem is restricted. Arvestad [22] later proposed a CDS alignment scoring model with a O(nm) alignment algorithm accounting for codon structures and FS translations based on the concept of generalized substitutions introduced in [23]. In this model, the score of a CDS alignment depends on its decomposition into a concatenation of codon fragment alignments, such that a codon fragment of a CDS is defined as a substring of length w, 0 ≤ w ≤ 5. This decomposition into codon fragment alignments allows to define a score of the CDS alignment at the AA level. More recently, Ranwez et al. [18] proposed a simplification of the model of Arvestad limiting the maximum length of a codon fragment to 3. Under this model, a O(nm) CDS alignment algorithm was described and extended in the context of multiple sequence alignment [18]. In the models of Arvestad [22] and Ranwez et al. [18], several scores may be computed for the same alignment based on different decompositions into codon fragment alignments. The corresponding algorithms for aligning two CDS then consist in computing an optimal score decomposition of an alignment between the two CDS. This optimal score exclusively accounts for FS translation initiations, i.e a FS translation in an alignment is penalized by adding a constant FS cost, which only penalizes the initiation of the FS, not accounting for the length of this FS translation. However, taking account of FS translation lengths is important in order to increase the precision of CDS alignment scores, as these lengths induce more or less disruptions between the protein sequences.
In this paper, we propose the first alignment algorithm that accounts for both the initiation and the length of FS translations in order to compute the similarity scores of CDS alignments. The remaining of the paper is organized as follows. In the "Motivation", we illustrate the importance of accounting for FS translation length when aligning CDS. In the "Preliminaries", we give some preliminary definitions and we introduce a new CDS alignment scoring model with a self-contained definition of the score of an alignment penalizing both the initiation and the extension of FS translations. In the "Methods", a dynamic programming algorithm for computing an optimal score alignment between two CDS is described. Finally, in the "Results", we present and discuss the results of a comparison of our method with other CDS alignment methods for a pairwise comparison of CDS from homologous genes of ten mammalian gene families.

Motivation: importance of accounting for FS translation length
The two main goals of aligning biological sequences are to evaluate the similarity and to identify similar regions between the sequences, used thereafter to realize molecular analyses such as evolutionary, functional and structural predictions. In practice, CDS alignment can be used to exhaustively identify the conserved features of a set of proteins. Thus, the definition of CDS similarity must account for sequence conservation and disruptions at both the nucleotide and the protein levels. Figure 1 illustrates the importance of accounting for AA translations and FS translation length in order to compute an adequate similarity score for a CDS alignment. It describes an example of three CDS Seq1, Seq2 and Seq3. Seq1 has a length of 45. The CDS Seq2 has length 60 and is obtained from Seq1 by deleting the nucleotide 'G' at position 30 and adding 16 nucleotides at the end. The CDS Seq3 has length 60 and is obtained from Seq1 by deleting the nucleotide 'G' at position 15 and adding 16 nucleotides at the end.
When looking at the AA translations of Seq1, Seq2 and Seq3, we observe that the similarity between Seq2 and Seq1 is higher than the similarity between Seq3 and Seq1 at the protein level, because Seq1 and Seq2 share a longer AA prefix "M T E S K Q P W H" (amino acids in black characters in the alignments). However, the pairwise CDS alignment algorithms that do not account for the length of FS translations would return the same score for the two following optimal alignments of Seq1 with Seq2 and Seq1 with Seq3, penalizing only the initiation of one FS translation in both cases (positions marked with a "!" symbol in the alignments), and not penalizing the sequence disruptions at the protein level.
From an evolutionary point of view, a good scoring model for evaluating the similarity between two CDS in the presence of FS translations should then penalize not only the initiation of FS but also the length of FS translations extension (amino acids in gray characters in the alignments). The alignment of Seq1 with Seq2 would then have a higher similarity score than the alignment of Seq1 with Seq3.

Preliminaries: score of CDS alignment
In this section, we formally describe a new definition of the score of a CDS alignment that penalizes both the initiation and the extension of FS translations.

Definition 1 [Coding DNA sequence (CDS)]
A coding DNA sequence (CDS) is a DNA sequence on the alphabet of nucleotides N = {A, C, G, T } whose length n is a multiple of 3. A coding sequence is composed of a concatenation of n 3 codons that are the words of length 3 in the sequence ending at positions 3i, 1 ≤ i ≤ n 3 . The AA translation of a CDS is a protein sequence of length n 3 on the alphabet A of AA such that each codon of the CDS is translated into an AA symbol in the protein sequence.
Note that, in practice an entire CDS begins with a start codon "ATG" and ends with a stop codon "TAA", "TAG" or "TGA". The following notations and conventions are used in Definition 3 to denote the different sets of codons and nucleotides in A and B. The set of IM codons in A (resp. B) is denoted by IM A→B (resp. IM B→A ). The set of FSext codons in A (resp. B) is denoted by FSext A→B (resp. FSext B→A ). The set of InDel codons in A (resp. B) is denoted by InDel A→B (resp. InDel B→A ). The set of MFS nucleotides in A (resp. B) is denoted by MFS A→B (resp. MFS B→A ). In these sets, the codons of A and B are simply identified by the position (column) of their last nucleotide in the alignment. In this case, we always have IM A→B = IM B→A as in the example below. The MFS nucleotides are also identified by their positions in the alignment.
For example, for the alignment depicted in Fig. 2 In the alignment scoring model described in Definition 3, the substitutions of IM and FSext codons are scored using an AA scoring function s aa such that aligned codons with silent nucleotide mutations get the same score as identity. A fixed FS extension cost denoted by fs_extend_cost is added for each FSext codon. The insertions/deletions of InDel codons are scored by adding a fixed gap cost denoted by gap_cost for each InDel codon. The alignment of MFS nucleotides are scored independently from each other, using a nucleotide scoring function s an . The insertions or deletions of nucleotides in FSinit codons are responsible for the initiation of FS translations. They are then scored by adding a fixed FS opening cost denoted by fs_open_cost for each FSinit codon. Note that, by convention, the values of all penalty costs for gap and FS (gap_cost, fs_ open_cost, fs_extend_cost) are negative. Note also that the scoring scheme assumes that the AA and the nucleotide scoring functions, s aa and s an , are symmetric.
The table D F is filled as follows:

Methods
In this section, we describe a O(nm) time and space complexity algorithm that solves the problem of finding a maximum score alignment between two CDS A and B of lengths n and m. Similarly to other classical sequence alignment algorithms [24], we use dynamic programming tables that are indexed by the pairs of prefixes of the two CDS. The table D stores the maximum scores of the alignments between prefixes of A and B. The table D F is used to account for potential cases of FS extensions that are counted subsequently.
Proof of Lemma 2 The proof follows from Lemma 1.
The five cases follow from the application of Lemma 1, case 2 for computing D(i + 1, j + 1), and by keeping only the cases where A[i + 1] and B[j + 1] are aligned together (cases 1, 2, 3, 4, 5 among the 11 cases). However, in each of the cases, we must subtract half of the score of aligning , because this score will be added subsequently.
The proof of Lemma 1 is given in the Additional file 1. Figure 3 illustrates the configurations of alignment considered in Lemma 1 for computing D(i, j) for cases 1 and 2. (a) iii.
(c) (a) iv. The alignment algorithm using Lemmas 1 and 2 is described in the next theorem.   For nm 9 calls of the case 1 of Lemma 2 For 2 × nm 9 calls of the cases 2 and 3 of Lemma 2 For 2 × nm 9 calls of the cases 4 and 5 of Lemma 2 Total = 12.55 nm

Results and discussion
We implemented the present CDS alignment algorithm with an affine gap penalty scheme [25] such that the penalty for a concatenation of k inserted (resp. deleted) codons is gap_open_cost + k * gap_cost, such that gap_open_cost is a negative penalty cost for gap initiations. This was done by adding two dynamic programming tables G A and G B such that the cell G A (i, j) (resp. G B (i, j)) contains the maximum score of an alignment between the prefixes A

Data
We evaluated the algorithm through applications on a mammalian dataset containing CDS sequences from ten gene families obtained from the database Ensembl-Compara version 83 [26]. The first gene family named "FAM86" is such that three CDS from three of its paralogous human genes were shown in [6] to share a common FS region translated in three different frames in the three CDS (see Fig. 4 for an illustration of the multiple alignment of these three CDS). The nine other families are the nine smallest (in term of the overall length of CDS) of fifteen gene families listed in [12] where they were shown to display one FS translation region between some pairs of CDS. For each gene family, the CDS of all human, mouse and cow genes belonging to the family and satisfying Definition 1 were downloaded. The overall number of distinct pairs of CDS within the ten gene families is 4011. Table 1 gives the details about the content and size of the ten gene families (The CDS of the ten gene families are provided in the Additional file 2).

Evaluation strategies
We compared the accuracy of five pairwise global alignment methods, including the present method, for computing CDS alignments in the presence or absence of FS translation between the compared CDS. The five methods vary according to the alignment algorithm used, either the present CDS alignment algorithm called FsePSA allowing to penalize both FS translation initiation and extension, or the CDS alignment algorithm called MACSE [18] penalizing FS translation initiation, or the Needleman-Wunsch (NW) sequence alignment algorithm [24] penalizing neither. Table 2 summarizes the alignment algorithm and the values of parameters used for each of the five methods.

Table 1 Detailed description of the ten gene families of the mammalian dataset
For each gene family, the family identifier used in [6] or [12], the Ensembl identifier of a human gene member of the family, the number of human, mouse and cow genes in the family, the total number of CDS of these genes, the total sum of lengths of these CDS and the number of distinct pairs of CDS are given The present CDS alignment algorithm is used in two of the five methods, namely fse and fse0. These two methods differ according to the value given to the parameter fs_extend_cost, either fs_extend_cost < 0 (−1, −0.5 or −0.2) for the method fse penalizing FS translation extension, or fs_extend_cost = 0 for the method fse0 not penalizing FS translation extension. The pairwise version of MACSE [18] is used in the method called macse_p. The NW alignment algorithm is used in the last two methods, the method called needlenuc computing scores and alignments at the nucleotide level and the method called needleprot at the AA level. For all methods using both the amino acid and nucleotide scoring functions s aa and s an , s an was fixed to +1/-1 for match/mismatch, so that the overall score of 3 consecutive nucleotide identities in an alignment scores less than the smallest identity score in s aa . All other parameters shared by several methods were given the same value for all methods. In particular, for the three methods fse, fse0 and macse_p penalizing FS translation initiation, the parameter fs_open_cost was given the values −10, −20 or −30. All other parameters were fixed to the default values for the NW algorithm implementation of NCBI Blast at the nucleotide and AA levels [27].

Gene family Human gene # of genes # of CDS Length
We used the five methods to compute pairwise alignments between the pairs of CDS within each of the ten gene families of our dataset, yielding 4011 alignments in total for each of five methods. In the absence of available benchmarks for the direct evaluation of the accuracy of CDS alignments, we base our evaluation on four indirect strategies.
In the first strategy, we consider the CDS multiple alignment of each gene family obtained using MACSE [18] as a benchmark. This strategy exploits the fact that multiple alignments are usually more accurate than pairwise alignments. It then assumes that the MACSE multiple alignments are closer to the reality than the pairwise alignments obtained using the five methods. Note that all the pairwise alignment methods included in the comparison can be extended to multiple sequence alignment methods using classical strategies. Thus, the more accurate pairwise alignment methods should lead to more accurate multiple alignment methods. Here, we focus on the comparison of the pairwise versions of the methods. In the second strategy, we consider six composition criteria for a CDS pairwise alignment called Iden-tity_NT, Identity_AA, Gap_init, Gap_length, FS_init, FS_length. The definitions of these criteria are given below, and used to compare the five methods. In the third strategy, we manually build and use as a benchmark, the real multiple alignment of three CDS from three paralogous human genes of the gene family I (FAM86). In the fourth strategy, we generate and use a set of three CDS splicing orthology groups, each group containing seven existing or putative CDS from seven genes of gene family I (FAM86).
Based on the results of the large-scale experiments discussed in the following, the best compromise for default values of FsePSA parameters are −30 for fs_open_ cost and −1 for fs_extend_cost.

Discussion
First strategy: using MACSE multiple alignments as benchmark MACSE [18] was used with its default parameters (fs_open_cost = −30, stop_cost = −100, s aa = BLOSUM62 matrix, gap_open_cost = −7 , gap_cost = −1) to compute the CDS multiple alignment of each of the ten gene families. For each MACSE multiple alignment of N CDS, we consider the N (N −1) 2 induced pairwise alignments as a benchmark. In total, we then obtained a benchmark composed of 4011 pairwise alignments. In order to compare an alignment (A ′ , B ′ ) obtained with one of the five methods to the corresponding alignment (A ′′ , B ′′ ) in the benchmark, we computed the number of nucleotides aligned in (A ′ , B ′ ) with the same partner as in the benchmark alignment (A ′′ , B ′′ ). Table 3 shows the overall percentage of nucleotides aligned with the same partners as in the benchmark for  2). It shows that the different versions of the fse method and the fse0 method have the best scores greater than 79.4%, followed by the needleprot method with a score of 78.82%. On the opposite, the needlenuc and macse_p method with fs_open_ cost= −30 return the worst scores, respectively 50.95% and 47.35%. These results also show that the fse method is more robust to the fs_open_cost parameter changes as compared to the macse_p method, whose scores show a large variation between 47.35 and 78.29%. Note that the needlenuc and needleprot do not account for the fs_open_cost parameter.

Second strategy: using six composition criteria for CDS pairwise alignment
Six criteria were defined and used to compare the five pairwise alignment methods. Given a pairwise CDS alignment, the first criterion Identity_NT counts the number of gap-free columns in the alignment containing a nucleotide match. The second criterion Identity_AA counts the number of IM and FSext codons c in the alignment that are aligned with a triplet of nucleotides yielding the same amino acid as c. The third criterion Gap_init is the number of gap-containing columns in the alignment, either insertion or deletion columns that are preceded by a different type of column. The fourth criterion Gap_length is the overall number of gap-containing columns in the alignment. The fifth criterion FS_ init is the number of FS translation segments found in the alignment. The last criterion FS_length is the overall number of columns in the alignment intersecting a FSext codon.
Note that the definitions of the six criteria exploit the definitions of codon sets used in Definition 3 but they are independent of any alignment scoring scheme. For example, for the alignment depicted in Fig. 2 For each of the nine cases obtained by combining the values of the parameters fs_open_cost (−10, −20 or −30) and fs_extend_cost (−1, −0.5 or −0.2), we considered the 4011 pairs of CDS from the ten gene families dataset, and partitioned them into three sets. For each case, the first set called the noFS dataset is composed of the pairs of CDS for which the pairwise alignments obtained using the fse0, fse and macse_p methods all have the criteria FS_init = 0. The second set called the FS dataset is composed of the pairs of CDS for which the alignments obtained using the fse0, fse and macse_p methods all have the criteria FS_init > 0. The third set called the ambiguFS dataset is composed of the remaining pairs of CDS.
Note that, in all nine cases, the set of CDS pairs for which FS_init = 0 with the macse_p method was strictly included in the set of CDS pairs for which FS_init = 0 with the fse method. For each of the nine cases, we computed the overall value of the six criteria for each method (fse0, fse, macse_p, needlenuc and needleprot) and each dataset (noFS, FS and ambiguFS). Tables 4, 5 and 6 present the results.

Results for the noFS datasets
For the noFS datasets, we assume that the real alignments should not contain FS translations. So, the needleprot method most likely computes the more accurate alignments since it does not allow any FS translation in the alignments. Indeed, it computes a maximum score NW alignment at the AA level and back-translates this alignment at the nucleotide level. We then take the needleprot result as a reference for the noFS dataset, in all cases. By construction of the noFS dataset, for a fixed value of the parameter fs_open_cost, the fse0 and fse methods necessarily return two alignments with the same similarity

Table 3 Comparison with MACSE multiple alignments benchmark
Percentage of nucleotides aligned with the same partner as in the benchmark alignments induced by the MACSE multiple alignments, for each method for varying fs_open_cost (−10, −20 and −30) and fs_extend_cost (−1, −0.5 and −0.2). In each case, the number of CDS pairs with an alignment that presents the highest similarity with the corresponding benchmark alignment as compared to the other methods is given in parenthesis. score for each pair of CDS of the dataset. Indeed, we observed that, for each value of fs_open_cost (−10 , −20 or −30), the alignments obtained using the methods fse0 or fse with varying values of the parameter fs_extend_cost are unchanged. Table  4 summarizes the results for fs_open_cost = −10, −20 and −30, presenting the results of the varying versions of fse and fse0 in a single line in the three cases. It shows that the results of the fse and fse0 methods are the closest to the reference for all the six criteria in all cases. However, they slightly overestimate or underestimate the criteria. The tendency of overestimating the Identity_AA and all other criteria is particularly accentuated for the macse_p method as compared to the fse and fse0 methods, in all cases. On the opposite, the needlenuc method always largely underestimates the Identity_AA, while overestimating all other criterion.

Results for the FS datasets
For the FS datasets, we assume that the real alignments must contain FS translations. So, the needleprot method can no longer produce the most accurate results. On the contrary, it is most likely that it underestimates the Identity_AA criterion. Indeed, it correctly aligns AA in CDS regions that are free of FS translation, but in FS translation regions, it either leads to several AA mismatches in the case of high mismatches scores, or to an overestimation of the Gap_init criterion. As expected, we observed that the value of Identity_AA for the needleprot method was always the lowest (data shown in the Additional file 3). We focus on the four other methods. Table 5 summarizes the results for the nine cases considered. For the Identity_NT and Identity_AA criteria, the differences between the values for the four methods are negligible. The main differences between the results reside in the values of the Gap_init and FS_init criteria. In particular, the FS_init criterion is useful to compare the accuracy of the methods for correctly identifying real FS translation regions. In [6] (for family I) and [12] (for families II-X), at most one FS translation region was detected and manually validated for each pair of CDS of the ten gene families. So, the expected number of FS translation regions per alignment in the FS data is 1. In Table 5, we observe that, in all cases, the fse and fse0 methods are the only methods for which the average numbers of FS_init are close to 1 with +/− standard error values smaller than 1. The macse method and especially the needlenuc method overestimate the number of FS translation regions per alignment with large standard error values in all cases.

Results for the ambiguFS datasets
For the ambiguFS datasets, all methods do not agree for the presence or absence of FS translation regions between the pairs of CDS. Note that the needlenuc method reports FS translations for all pairs of CDS, with the highest average number of FS translation regions per alignment in all cases (data shown in the Additional file 3). As needlenuc is already shown to perform poorly in both the absence and the presence of FS translation regions, we focus on the four other methods. Table 6 summarizes the results. We observe that, for all criteria, macse_p has higher values than fse0, fse and needleprot that have similar values. The most significant difference between the results resides in the values for the FS_init and FS_length criteria. The fse method always reports a null or a very small number of FS regions with an average FS_init equals to 1 as expected. In all cases, the fse0 and macse methods overestimate the number of FS translation regions per alignment.

Third strategy: using a 3-CDS manually-built benchmark
We manually built the real pairwise alignments of three CDS from three paralogous human genes of gene family I, the CDS FAM86C1-002 coding for protein  Figure 4 also shows that each pair of CDS shares a single FS translation region. Table 7 shows the normalized pairwise similarity scores and the number of FS translation regions computed by the five alignment methods (the pairwise alignments computed by the five methods with varying fs_open_ cost and fs_extend_cost are given in the Additional file 4). It shows that needleprot and fse (in Table 6 Values of the six criteria for the ambiguFS dataset  Table 7 also illustrates the fact that needlenuc and macse_p strongly overestimate the number of FS translation regions per alignment in all cases. The fse method with the parameters fs_ open_cost= −10 and fs_extend_cost= −1 is the only method that allows to infer that FAM86C1-002 and FAM86B1-001 are the most similar and to detect a single FS translation region for each alignment.

Fourth strategy: inferring CDS splicing orthology groups and protein phylogenies
Based on the three CDS used in the previous strategy, CDS FAM86C1-002 from human gene ENSG00000158483, FAM86B1-001 from human gene ENSG00000186523 and FAM86B2-202 from human gene ENSG00000145002, we generated a dataset of three CDS splicing orthology groups composed of 21 homologous CDS. Each group contains one of the three initial CDS and its six splicing orthologs in the following set of seven genes from gene family I: human genes ENSG00000158483 denoted H1, ENSG00000186523 denoted H2 and ENSG00000145002 denoted H3, each containing one of the initial CDS, chimpanzee gene ENSPTRG00000007738 denoted Ch, mouse gene ENSMUSG00000022544 denoted M, rat gene ENSRNOG0-0000002876 denoted R and cow gene ENSBTAG00000008222 denoted Co. The CDS splicing orthologs were predicted based on the spliced alignment tool Splign [28] as follows: for each initial CDS A 1 of a gene A and each gene B different from A, A 1 was aligned to B and a putative or existing CDS of B ortholog to A 1 with the same splicing structure was inferred. The 21 resulting CDS are given in Additional file 5.
We computed the normalized pairwise similarity scores between the CDS, using the five alignment methods (the pairwise alignments computed by the five methods with varying fs_open_cost and fs_extend_cost are given in the Additional file 5). For each method, we constructed a phylogeny using an UPGMA and a Neighbor-Joining (NJ) algorithm, based on the computed CDS similarity matrix. The UPGMA algorithm was used to classify the CDS into three groups and infer the similarity relationships between the groups independently of any rate of evolution. The NJ algorithm was used to reconstruct the phylogeny inside each group. Table 8 summarizes the results. The three splicing orthology groups are denoted G1 (containing CDS C1-002), G2 (containing CDS B1-001) and G3 (containing CDS B2-202).
All methods allow to correctly classify the CDS into the three initial splicing orthology groups G1, G2, and G3. However, the needleprot and fse methods are the only methods that allow to infer the correct similarity relationships ((G1,G2),G3) between the groups, confirming the results of the third evaluation strategy. For all methods, the CDS phylogeny reconstructed inside the group G2 is (Co,((M,R),((H1,Ch),(H2,H3)))) inducing an evolution of the seven genes with a speciation event at the root of the gene tree. The phylogeny reconstructed for