{'fraction_non_alphanumeric': 0.028811754103285358, 'fraction_numerical': 0.007455553431466259, 'mean_word_length': 4.444246208742195, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving spliced alignment for identification of ortholog groups and multiple CDS alignment

The Spliced Alignment Problem (SAP) that consists in finding an optimal semi-global alignment of a spliced RNA sequence on an unspliced genomic sequence has been largely considered for the prediction and the annotation of gene structures in genomes. Here, we re-visit it for the purpose of identifying CDS ortholog groups within a set of CDS from homologous genes and for computing multiple CDS alignments. We introduce a new constrained version of the spliced alignment problem together with an algorithm that exploits full information on the exon-intron structure of the input RNA and gene sequences in order to compute high-coverage accurate alignments. We show how pairwise spliced alignments between the CDS and the gene sequences of a gene family can be directly used in order to clusterize the set of CDS of the gene family into a set of ortholog groups. We also introduce an extension of the spliced alignment problem called Multiple Spliced Alignment Problem (MSAP) that consists in aligning simultaneously several RNA sequences on several genes from the same gene family. We develop a heuristic algorithmic solution for the problem. We show how to exploit multiple spliced alignments for the clustering of homologous CDS into ortholog and close paralog groups, and for the construction of multiple CDS alignments. An implementation of the method in Python is available on demande to SFA@USherbrooke.ca. Keywords: Spliced alignment, CDS ortholog groups, Multiple CDS alignment, Gene structure, Gene family.


Introduction
Spliced alignment consists in aligning spliced gene products against unspliced DNA sequences [8]. It is a key step for genome annotation, gene prediction, identification of gene structures and alternative splicing studies [7,17]. The quality of the transcriptome annotation is greatly improved with the use of spliced RNA data of related genes or genomes. In the past two decades, several spliced alignment tools have been developed for aligning RNA or protein data on genomic DNA ( [9,16] for example). Such methods look for an alignment that maximizes the sequence similarity and splice site consensus signals between a short RNA sequence and a large genome sequence. They are specifically developed for genome annotation purpose, without any assumption on the homology of the compared sequences. So, in order to achieve high sensitivity, they look for high-identity alignments in order to avoid as much false positive alignments as possible. As a consequence, these methods perform accurately for the comparison of sequences from closely related genomes but often misses the alignment of homologous sequences from distant genomes.
Although spliced alignment methods are classically used for the prediction and annotation of genes in genomes, they can also be specifically designed for the purpose of aligning sequences from a priori homologous genes. In this context, they can be used as a preliminary step for the identification of orthologous transcripts or proteins isoforms [17], the transfer of transcript annotations between homologous genes [1] or for studying the evolution of RNA isoforms in a gene family [4,10,12]. Spliced alignment can also improve the accuracy of multiple homologous CDS alignment by making use of one or several reference genes. Current spliced alignment methods make use of the gene structure by accounting for possible splice sites predicted in the gene sequence. However, accounting for the exact exon-intron structure of the gene, when it is known, can help improve the alignment further. Moreover, current methods do not make use of the exon structure of the query RNA sequence because they were mainly designed for aligning EST and RNA-seq data, for which the exon structure is unknown, on genomic sequences. However, when the RNA sequence is given with its corresponding gene sequence, it is possible to easily recover its exon structure and use this information for the spliced alignment against an other gene.
In this paper, we are interested in the Spliced Alignment Problem (SAP) that consists in finding an optimal spliced alignment between a CDS and a gene that captures all sequence similarities including those between phylogenetically distant homologous sequences. For aligning a CDS against a gene, under the assumption of homology between the sequences, it is possible to make use of the splicing structure of the sequences in order to accurately detect sequence homologies even in the case of phylogenetically distant sequences. First, in Section 3, we introduce a constrained version of the SAP problem and we present a CDS-gene spliced alignment method that exploits full information on the structure of the input CDS and the input gene in order to compute high-coverage accurate alignments. Second, in Section 4, we show how this method can be directly used for clustering the set of CDS of a gene family into ortholog and close paralog groups, allowing one-to-one as well as one-to-many orthology relations. Third, in Section 5, we introduce an extension of the Spliced Alignment Problem called Multiple Spliced Alignment Problem (MSAP) that consists in finding an optimal multiple spliced alignment between a set of CDS and a set of genes. The MSAP problem also extends a homonym problem introduced in [8] for the simultaneous spliced alignment of a set of CDS on a single gene. We describe an algorithm for the MSAP problem that consists in combining pairwise spliced alignments into a multiple spliced alignment, by merging progressively the pairwise alignments while filtering out false-positive sub-alignments. We show that the multiple spliced alignments allow to identify homologous segments (exons) across a set of homologous genes and to compute accurate CDS ortholog groups and multiple CDS alignments. Finally, Section 6 is devoted to the evaluation of our method by comparing it to other methods based on the application to the analysis of set of homologous CDS and genes from the Ensembl-Compara database.
2 Preliminaries: genes, CDS, splicing and orthology In this section, we give some formal definitions that will be useful for the remaining of the paper. Given a set S, |S| denotes the size of S, and given a sequence or an interval T , length(T ) denotes the length of T . Gene, Exon and CDS: A gene is DNA sequence on the alphabet of nucleotides Σ = {A, C, G, T }. Given a gene G of length n, an exon of G is a pair of integers (a, b) such that 1 ≤ a ≤ b ≤ n. The sequence of an exon (a, b) of G, denoted by G[a, b], is the segment of G identified by its start and end positions a and b in G.
Thus, the exons of C are non-overlapping and totally ordered by increasing location on the gene. We denote by C[i] the i th element of C. We denote the set of introns induced by C by The sequence of a CDS C of G, denoted by G C , is the concatenation of the sequences of gene exons composing C in the order in which they appear in C. So, the length of the CDS sequence G C equals the sum of the length of the gene exon sequences composing G C . By extension of the definition of a gene exon, an exon of a CDS sequence G C of length m is a pair of integers (k, l) such that 1 ≤ k ≤ l ≤ m and the segment G C [k, l] of G C , identified by its start and end positions k and l in G C , is exactly one of the gene exon sequences composing G C . We denote by E(G C ) the set of CDS exons composing a CDS sequence G C . Note that |E(G C )| = |C|. For the example of gene G with a CDS C 1 given above, |G C1 | = 30 and E(G C1 ) = {(1, 10), (11,30)}. Gene family, orthology relationships: A gene family G is a set of homologous genes that have derived from the same original gene by duplication and speciation events. The genes of a gene family are supposed to have similar segments and biochemical functions. The evolution of a gene family is represented by a rooted binary tree T whose set of leaves is G and internal nodes are ancestral genes that precede an ancestral event labeled as duplications or speciations. A node is an ancestor of a gene (leaf) if it is on the path between this leaf and the root of the tree. The lowest common ancestor (LCA) of two genes is the common ancestor to both genes that is the most distant from the root. The orthology relationship between genes can be defined from two different points of view. From a gene tree point of view, two genes of G are called orthologs if their lowest common ancestor in the tree T is a speciation, otherwise they are called paralogs. From a similarity point of view, orthologous genes have similar sequences, structures and functions. For instance, reciprocal best hits are a common approach for the definition of orthology relations in comparative genomics. Splicing and spliced alignment: In molecular biology, a eukaryotic mature RNA sequence is obtained from a RNA transcript by a phenomenon called splicing that consists in removing the sequences between any two successive exons of a CDS called introns and joining the resulting exon extremities. The CDS is the segment of the mature RNA sequence that is translated into a protein sequence. In practice, a full CDS sequence has a length that is multiple of 3 and it starts with a codon "ATG" and ends with a codon "TAA", "TAG" or "TGA".
A spliced alignment is an alignment between a CDS and a gene sequence that allows to identify conserved exons sequences. See Figure 1 for example. Formally, a spliced alignment of a CDS sequence G C of length m on a gene H of length n is a chain of quadruplets A = {(k 1 , l 1 , a 1 , b 1 ), . . . , (k j , l j , a j , b j )} called blocks such that for any block (k, l, a, b) of A, 1 ≤ k ≤ l ≤ m, a = b = 0 or 1 ≤ a ≤ b ≤ n, and for any two successive blocks (k The blocks (k, l, a, b) of A such that 1 ≤ a ≤ b ≤ n correspond to conserved exon segments G C [k, l] and H[a, b] between the CDS and the gene sequence. We call them conserved blocks and we denote the set of conserved blocks of A by The blocks (k, l, a, b) such that a = b = 0 correspond to exon segments G C [k, l] in the CDS sequence that are absent in the gene sequence. We call them deleted blocks. In other terms, the blocks composing A correspond to a chain of non-overlapping segments of the CDS sequence G C that are increasingly located on G C and cover entirely G C . Moreover the conserved blocks correspond to a chain of non-overlapping segments of the gene sequence H that are also increasingly located on H. Thus, there is no crossing in the alignment. We denote by A[i] the i th block of A.
The spliced alignment A also induces a set of putative gene intron segments. These intron segments are the gene segments that lie between two successive blocks of A that are both conserved in the gene sequence. We denote the set of introns induced by the alignment A = {(k 1 , l 1 , a 1 , b 1 ), . . . , (k j , l j , a j , b j )} by Note that if all blocks composing A are conserved, then |Intron(A)| = |A| − 1.

Computation of spliced alignment
In this section, we first re-call two well-known versions of the spliced alignment problem and we introduce a third version that we study in this paper. Next, we describe a heuristic algorithm for the last problem. denote the score of an optimal global alignment between them. The following is a reformulation of the less constrained version of the spliced alignment problem studied in [9]. This formulation allows a unified framework to compare the problem with some more constrained versions of the problem formulated thereafter.

Spliced Alignment Problem I (SAP I):
The (SAP I) problem only accounts for the alignment scores between the segments composing the blocks of the spliced alignment. In practice, in more than 99% of real cases of splicing, the removed intron sequences start with a dinucleotide sequence "GT" and ends with a dinucleotide sequence "AG" respectively called canonical donor and acceptor splice sites [3,15]. Thus, in order to improve the accuracy of the spliced alignment, a more constrained version of the problem allows to account for the intron segments induced by an alignment. Given a nucleotide sequence S, let intr(S) denote the intron score of S accounting for the presence or absence of canonical splice sites at the extremities of S. A sequence with two canonical splice sites at its extremity has a higher score than a sequence with only one which has a higher score than a sequence without canonical splice sites. A more constrained version of the spliced alignment problem studied in [9] is the following.

Spliced Alignment Problem II (SAP II):
Input: A CDS sequence G C from a gene G ; a gene sequence H.
The SAP I and SAP II problems do not account for the exon structure of the CDS sequence and the exon-intron structure of the gene sequence. In order to further improve the accuracy of the spliced alignment, we consider a more constrained version of the problem that accounts for the actual exons in the CDS and the gene sequence. Given a conserved block (k, l, a, b) of a spliced alignment, let exon E(G C ),E(H) (k, l, a, b) denote the score of a conserved block accounting for the correspondence of the block with an actual exon in the CDS or in the gene sequence. A block that corresponds to an exon in both sequences has a higher score than a block with an exon correspondence in only one of the sequences which has a higher score than a block without any exon correspondence. The more constrained version of the problem is defined as follows.

Spliced alignment Problem III (SAP III):
Heuristic algorithm for the (SAP III) problem: The general approach followed by most heuristic spliced alignment methods for the SAP I and SAP II problems consists in first computing high-identity local alignments between the CDS and the gene sequence. In a second step, these local alignments are used as anchors and a global alignment algorithm is applied in the regions between the anchor alignments in order to complete the spliced alignment. Our heuristic spliced alignment method for the SAP III problem also starts with the computation of highly conserved local alignments used as anchors. However, in a second step, we use the exon structure of the CDS sequence to extend the anchors towards the extremities of the CDS exons containing them. This steps drastically reduces the length of the CDS sequences remaining between the extended anchor in the deleted blocks. In a third step, we apply a global alignment in the remaining regions. In a fourth step, we correct the splicing junctions in order to recover GT/AG canonical splice sites that were missed because of the insertion or deletion of nucleotides in the exon sequences, resulting in gaps in the block alignments. Finally, in a fifth step, we use the exon-intron structure of the gene sequence to further correct the block extremities. The details of the steps composing the method are given below.
Step 1. Local alignment. This step is achieved using Translated Blast (tblastx) with an E-value threshold of 10 −2 in order to obtain local alignments. Tblastx is used in order to account for the translation of the sequences into amino acid sequences. This allows to detect amino acid sequence conservation even in the presence of translational frameshifts.
Step 2. Block extension based on CDS exons. This step is repeated iteratively in order to prioritize highly conserved local alignments. Local alignments are classified into four groups according to their E-value scores, less than 10 −7 , between 10 −7 and 10 −5 , between 10 −5 and 10 −3 and greater than 10 −3 . For each group from the more to the less conservative, a maximum size set of pairwise compatible anchors is kept as conserved blocks and the blocks are extended as follows. Given a block extremity that does not correspond to a CDS exon extremity, the block extremity is extended until the closest CDS exon extremity if this extension does not decrease the percentage of nucleotide identity in the block alignment. If the extension induces an overlapping with another conserved block on the CDS sequence, the other block is trimmed in order to free space for the extension. In this case, the extension and the trimming are applied only if they strictly increases the nucleotide identity of the two block alignments. See Figure 2 for an illustration. For each comparison, we make use of an exact semi-global alignment that penalizes no end gaps in the two compared sequences.
Step 4. Correction of exon junctions. In this step, the introns induced by the spliced alignment are refined. Given an induced intron (b, a) between two conserved blocks of the alignment, if the donor and acceptor splice sites of the intron are not both canonical, we look for a shift of the block junction that corresponds to a pair of canonical splice sites and does not decrease the percentage of nucleotide identity in the flanking block alignments. The space search is limited to a maximum of 30-nucleotide shift and 3-codon gaps.
Step 5. Correction of block extremities based on gene exons. This step is composed of two phases. First, each extremity of a conserved block that does not correspond to a gene exon extremity or a CDS exon extremity is trimmed until the closest gene exon extremity if there exists such a gene exon extremity. This allows to free space in order to extend neighboring conserved block extremities in a second phase. In the second phase, each conserved block extremity that still does not correspond to a gene exon extremity or a CDS exon extremity is extended until the closest available gene exon extremity if this extension does not decrease the nucleotide identity of the block alignment. See Figure 2 for an illustration.
Note that, the algorithm used for each step of the method can be replaced by any other algorithm solving the same problem. For instance, Step 1 can use a faster local alignment algorithm than tblastx. Steps 3 and 4 can be merged H" Gene"exon" Gene"exon" Trimming" Extension" CDS"exon" CDS"exon" Trimming" by using a global alignment algorithm that accounts simultaneously for the sequence similarities and the intron scores such as the global alignment algorithm developed [9]. The output of the spliced alignment method is a spliced alignment whose blocks extremities maximize the correspondence with the exon extremities of the input CDS and gene sequences.

Identification of CDS ortholog groups
In this section, we give a definition of CDS ortholog groups based on pairwise spliced alignments. We first start with a definition of orthologous CDS.
In [17], an extension of the concept of gene orthology to spliced transcript orthology was introduced. They defined orthologous transcripts as two structurally similar transcripts from two orthologous genes. In [12], we have introduced a new protein tree model of transcript evolution along gene trees, that is similar to the gene tree model of gene evolution along species tree with duplication and speciation events. In the model of transcript evolution along gene trees, we have introduced a third type of event called creation representing the separation of two lineages of structurally different transcripts. This model have led to an extended definition of protein and transcript orthology, based on the protein tree model and relaxing the constraint that two orthologous transcripts should come from orthologous genes. We call two transcripts ortholog if their LCA in the protein tree is a speciation or a duplication (not a creation), otherwise they are called paralogs. From a similarity point of view, the corresponding definition of transcript orthology is that two transcripts are orthologs if they are structurally similar and come from two homologous genes, not necessarily orthologous genes. Next, two orthologous transcripts are ortho-orthologs if they come from orthologous genes, and para-orthologs otherwise. Note, that our definition of ortho-orthologs then corresponds to the definition of orthologous transcripts from [17]. Thus, the orthology relationship between two transcripts from two homologous genes relies on the evaluation of the structural similarity between the transcripts. Here, we evaluate this structural similarity using the CDS associated to the transcripts.
CDS orthology: Let C 1 and C 2 be two CDS from two homologous genes G and H respectively. Let A 1 be a spliced alignment of the CDS sequence G C1 on the gene sequence H, and A 2 a spliced alignment of H C2 on G. Using A 1 and A 2 , we define C 1 and C 2 as orthologs if: In other terms, C 1 are C 2 as orthologs if (1) they have the same number of exons, (2) the spliced alignment of C 1 on H induced the same introns as C 2 in H or the spliced alignment of C 2 on G induces the same introns as C 1 in G and (3) the lengths of each pair of corresponding exons in C 1 and C 2 should be congruent modulo 3. The conditions (1) and (2) ensure that the two CDS have the same exon structure. The condition (3) ensures that the two CDS are translated in the same codon phase in each pair of corresponding exons.
Note that this definition only requires that one of the spliced alignments A 1 and A 2 supports the orthology relation. An alternative more stringent definition of CDS orthology consists in requiring the reciprocity, i.e. that A 1 and A 2 both support the orthology relation, by using the 'and' statement instead of the 'or' statement in the condition (2). CDS ortholog groups: Given a set of CDS C from a set of homologous genes G, the transitivity of the CDS orthology relation is used to identify distant orthologs and co-orthologs in C that cannot be identified by means of the CDS structural similarity. Such orthologs are typically missed because of partial spliced alignments due to low sequence similarity.
The CDS orthology relation on C is then extended into an equivalence relation such that for any three CDS C 1 , C 2 , C 3 in C, if C 1 and C 2 are orthologs and C 2 and C 3 are orthologs, then C 1 and C 3 are also orthologs. The CDS ortholog groups are defined as the equivalence classes of the resulting equivalence relation.

Computation of multiple spliced alignment
The concept of spliced alignment was extended in order to define multiple spliced alignments [2,8]. A multiple spliced alignment was then defined as a simultaneous spliced alignment of several RNA sequences on a single gene sequence. Multiple spliced alignment allows to improve the accuracy of gene structure prediction by making use of several RNA targets simultaneously. Here, we further extend the concept in order to consider the simultaneous spliced alignment of several CDS sequences on several genes. See Figure 3 in Appendix for example.
Based on the formalism introduced in Section 2, we defined a multiple spliced alignment of a set of CDS sequences C on a set of genes G as a chain A = {A [1], . . . , A[j]} of sets called multi-blocks such that . Given a spliced alignment A x,y of a CDS sequence x on a gene sequence y, let S(A x,y ) denote the score of the spliced alignment A x,y as defined for the SAP I, SAP II or SAP III problem. We defined the corresponding multiple spliced alignment problem as follows: Multiple spliced alignment Problem (MSAP S): Input: a set of CDS sequence C ; a set of genes G. Output: A multiple spliced alignment A of C on G that maximizes the sum of the scores of induced pairwise spliced alignments: (x,y)∈C×G S(A x,y ).
Heuristic algorithm for the MSAP III problem: We now describe a heuristic algorithm for building a multiple spliced alignment of a set of CDS sequence C on a set of genes G, given the spliced alignments for each pair of sequences (x, y) ∈ C × G computed using a pairwise spliced alignment method.
The idea behind the algorithm is to progressively merge the blocks of the input pairwise spliced alignments into multi-blocks of the target multiple spliced alignment. In order to merge a block into an existing multi-block, we rely on the compatibility between the segments composing the blocks. Given a sequence x, and two segments of x, x[s 1 , e 1 ] and x[s 2 , e 2 ], we say that the segments are compatible if s 1 = s 2 ± and e 1 = e 2 or s 1 = s 2 and e 1 = e 2 ± . In other terms, the segments are compatible if they are nested with an equality of at least one of the two extremities of the segments, and the other extremities have a difference of at most . If x[s 1 , e 1 ] and x[s 2 , e 2 ] are compatible, we denote by max((s 1 , e 1 ), (s 2 , e 2 )) the largest of the two intervals. In practice is set to 50.
Next, we define the compatibility between a pairwise block and multi-block, and the compatibility between two multi-blocks as follows. Let A be a multiple spliced alignment of C on G obtained at some step of the algorithm. Let (k, l, a, b) be a block in a spliced alignment of a pair of sequence (x, y) ∈ C × G. The block  (k, l, a, b) is compatible with a multi-block (1) Let (k, l, a, b) be a block in the spliced alignment of x ∈ C on y ∈ G considered at some iteration of the algorithm. By construction, the block (k, l, a, b) can be compatible with zero, one or two multi-blocks of A. Depending on the case, the multiple spliced alignment A is refined as follows.  ((a, b), (s y i , e y i )). Case 3. If (k, l, a,  is erroneous and it should be discarded. One should be very careful in choosing one of these cases in order to avoid subsequent conflicts caused by a wrong decision at this step. In order to decide of the correctness of the block, the percentage of nucleotide identity in the block alignment (k, l, a, b) is one criteria but it is not sufficient. Indeed, when an exon is duplicated in a gene (as we will see in one example in Section 6), a block with a high percentage of nucleotide identity in the alignment can be erroneous. Consequently, we give a very conservative definition of a correct block in the case of conflict. We say that the block is correct if the percentage of nucleotide identity in its alignment is greater or equal to a threshold τ and the block alignment contains no gaps. In practice τ is set to 60%. Case 3.b.i. If the block (k, l, a, b) is not correct then it is not added to the multiple spliced alignment. Case 3.b.ii. Otherwise, we define some scores for the occurrence of the segment (s x i , e x i ) in A[i] and the occurrence of the segment (s y i , e y i ) in A[j]. The occurrence whose score is the higher is kept and the block (k, l, a, b) is added in the corresponding multi-block, while the other occurrence is discarded from its multi-block. The output of this algorithm is a multiple spliced alignment whose multiblocks extremities maximize the correspondence with the exon extremities of the input CDS and gene sequences. We now describe how to use a multiple spliced alignment in order to cluster the CDS of set of homologous genes into groups of orthologs and close paralogs. Otherwise, x and y are co-orthologs or orthologs. Multiple CDS alignment: The multiple spliced alignment A of the set of CDS sequences C on the set of genes G is also used to define a Multiple Sequence Alignment (MSA) of the CDS sequences and the genes exon segments. For each multi-block A[i] ∈ A, the set of segments {x[s x i , e x i ] |x ∈ C ∪ G } is aligned using a MSA tools, and the resulting alignments are concatenated in order to obtain a global multiple alignment of all the CDS sequences with the concatenation of the gene exon segments. In practice, we use the sequence aligner Muscle [6] for the multiple alignment of the segments composing a multi-block.

Application
We applied our algorithms on two sets of homologous genes from the Ensembl-Compara database release 89 [5]. We wanted to evaluate the aptitude of the methods (1) to compute pairwise high CDS-coverage spliced alignments maximizing the correspondence with actual exon extremities in the CDS and gene sequences, (2) to identify accurate CDS ortholog groups covering phylogenetically distant genes based on pairwise spliced alignments, and (3) to compute accurate multiple CDS alignments, and ortholog and close paralog groups based on multiple spliced alignments. In the remaining of the section, our methods are named SpliceFAmAlign (SFA).
Dataset: The dataset contains 16 genes with their CDS sequences from two gene families, FAM86 and MAG: 8 genes per family, 14 CDS for FAM86 and 26 for MAG. For each family, the genes are from seven different amniote species which are human, chimpanzee, mouse, rat, cow and chicken (for FAM86) or lizard (for MAG). In particular, the 3 MAG human genes contain a pair of duplicated exons separated by a third smaller exon. Table 1 in Appendix gives more details about the dataset.

Evaluation of the pairwise spliced alignments:
We compared the results of our pairwise spliced alignment algorithm SFA for the SPA III problem with the results of the current best tool SPlign [9] developed for the SPA II problem. We wanted to compare the ability of the methods to correctly identify actual exon extremities in the CDS and the gene sequence, canonical splice sites and actual introns in the gene sequence. We also evaluated the CDS coverage of the spliced alignments computed using the two methods. For the two methods, we considered different values for a minimal block nucleotide identity parameter called min idty such that the blocks with lower identity ratio than this threshold in their alignment are discarded from the spliced alignment. For SPlign, we tested the values 0.0 and 0.75 (default parameter), and for SFA the values 0.0, 0.6 (default parameter), 0.7, and 0.75 for min idty. The full results are shown in Appendix in Table 2. As expected, the spliced alignments computed by SFA cover a larger percentage of the CDS sequences. For the ability to correctly recover actual exon extremities in the CDS and the gene sequences, SFA recovers almost twice more actual exon extremities than SPlign, but its number of blocks is also doubled. Thus, the ratios of actual CDS or gene exons extremities, and the ratio of canonical splice sites is lower for SFA.
2. Evaluation of the CDS ortholog groups computed based on pairwise spliced alignments: We applied our ortholog clustering method using pairwise spliced alignment in order to obtain a set of ortholog groups based on structural similarity. We also computed a set of ortholog clusters using the OrthoMCL clustering tool [13] based on sequence similarity. We wanted to evaluate the ability of the SFA method to identify accurate CDS ortholog groups covering phylogenetically distant genes, and to compare the structure-based ortholog clusters obtained with the SFA method to the clusters obtained using a sequence-based clustering method.
The sets of ortholog clusters computed by the two methods are depicted in Figures 7 and 5 in Appendix. For FAM86, SFA allows to recover a 5-CDS ortholog group composed of CDS conserved across 2 human, 1 mouse, 1 rat and 1 cow genes. It also returns a pair of orthologous CDS from 2 orthologous human and chimpanzee genes. For the MAG family, SFA allows to recover a 3-CDS ortholog group with CDS from 1 human, 1 rat and 1 cow genes, and a pair of orthologous CDS from 2 orthologous human and rat genes. The multiple alignment of the corresponding protein sequences confirms the accuracy of the ortholog groups (data not shown). The structural similarity between the orthologous CDS is also highlighted in Figures 7 and 5 (Appendix). Concerning the comparison between the OrthoMCL and SFA clusters, for FAM86, SFA returns 9 clusters and OrthoMCL returns 5 clusters. For MAG, SFA returns 23 clusters while OrthoMCL returns 3 clusters. In both cases, the two sets of clusters are compatible as the clusters computed by the SFA methods are included in larger clusters computed by OrthoMCL. Figures 7 and 5 (Appendix) reveals that the SFA clusters merged by OrthoMCL share some sequence similarities indeed, but they also have different exon structures.
3. Evaluation of the multiple spliced alignment: We applied our multiple spliced alignment algorithm using the pairwise spliced alignments as input. Next, we computed a set of ortholog and close paralog groups, and a multiple alignment of the CDS sequences using the methods described in Section 5. First, we wanted to evaluate the accuracy of the clusters reconstructed based on the multiple spliced alignment and their compatibility with the clusters reconstructed using the more stringent method based on pairwise spliced alignments. For both families, all pairwise SFA clusters are included in larger clusters defined using the multiple SFA clustering method (data not shown). Second, we wanted to evaluate the accuracy of the reconstructed multiple CDS alignments. To do so, we compare them with the multiple CDS alignments obtained using the only existing tools MACSE [14] that allows to align multiple CDS sequences while accounting for translational frameshifts. MACSE [14] was ran with its default parameter values to compute a multiple CDS alignment for each of the gene families. For each gene family, we compare the two CDS multiple alignments, the one from our method SFA and the one from MACSE, based on two criteria which are the percentage of column identity in the alignment and the ratio of long gaps in the alignment that correspond to real exon-exon junctions in the CDS sequences. We define a long gap as a gap of length greater or equal to 20 between two nucleotides x[i] and x[i + 1] of a CDS sequence x in the multiple CDS alignment. This gap corresponds to a real exon-exon junction in the CDS if x[i] and x[i + 1] belong to two different (consecutive) exons in x.
The results are shown in Table 3 in Appendix. In all cases, SFA recovers a higher number and a higher ratio of long gaps corresponding to real exon junctions in the CDS sequences.

Conclusion
This paper introduces a new constrained version and an extended version of the Spliced Alignment Problem for the purpose of identifying CDS ortholog groups within a set of CDS from homologous genes, and for computing multiple CDS alignments. We have developed a variety of algorithmic solutions for the computation of pairwise and multiple spliced alignments. We show how this framework can be used to improve the definition of CDS orthology clusters and multiple CDS alignments. The application of the algorithms to real datasets shows that the framework is particularly useful for identifying CDS ortholog groups that are conserved accross phylogenetically distant genes.
On the algorithmic side, the new problems presented require a more in-depth investigation. For instance, the heuristic algorithm described for the constrained spliced alignment problem could make use of a exact global alignment method with scoring schemes defined to account for all the constaints of the alignments. The heuristic algorithm described for the multiple spliced alignment problem follows a classical greedy approach used by most multiple sequence aligner based on pairwise alignment. What is the complexity of the problem ? It can be seen as a problem of reconstructing a macroscopic multiple alignment and then it is likely to be NP-hard. If so, it is possible that the problem is fixed-parameter tractable with respect to a parameter such as the maximum number of blocks in a pairwise alignment or the maximum size of a conflicting set of blocks. Future work will also make use of benchmark datasets in order to confirm the experimental results of this study.  [3], A [5] have occurrences in all sequences, while A [2] is conserved in G, H and HC 4 and A [4] is only present in G and GC 2 . The pairwise spliced alignment induced for the CDS sequence GC 2 on the gene sequence H is depicted in gray color. For each gene family, the family identifier, the Ensembl identifier of the genes, the number of complete CDS for each gene and the species of genes are given. Fig. 4. The sets of ortholog clusters obtained using the SFA and the OrthoMCL methods for the FAM86 family. The clusters containing more than one CDS are depicted in red color for the SFA method, and in dotted blue colored line for OrthoMCL. The visualization of the alignment was generated using CLC Sequence Viewer [11] . Table 2. Results for the comparison of SFA with the pairwise spliced alignment method SPlign [9].      6. The multiple CDS alignment obtained using the SFA method for the FAM86 family. The visualization of the alignment was generated using CLC Sequence Viewer [11]. Fig. 7. The multiple CDS alignment obtained using the MACSE method [14] for the FAM86 family.