On calculating the probability of a set of orthologous sequences.

Probabilistic DNA sequence models have been intensively applied to genome research. Within the evolutionary biology framework, this article investigates the feasibility for rigorously estimating the probability of a set of orthologous DNA sequences which evolve from a common progenitor. We propose Monte Carlo integration algorithms to sample the unknown ancestral and/or root sequences a posteriori conditional on a reference sequence and apply pairwise Needleman–Wunsch alignment between the sampled and nonreference species sequences to estimate the probability. We test our algorithms on both simulated and real sequences and compare calculated probabilities from Monte Carlo integration to those induced by single multiple alignment.


Introduction
Comparative genomics/proteomics research often focuses on a set of orthologous sequences arising from evolutionary speciation. For example, multiple related species (for example, human, mouse, and rat) can have a common gene as well as the corresponding promoters in the upstream region of such a gene, although these matched sequences may have minor difference across species. For simplicity the set of sequences studied in the sequel are assumed to have almost equal length in light of these examples. Sequence alignment algorithms 1 have substantially facilitated comparative genomics/proteomics research by showing conservation pattern along orthologous sequences, and biologically functional segments are likely to be those more conserved regions along the genome. For the vast body of related literature, we refer to Liu and colleagues, 2 Kellis and colleagues, 3 Moses and colleagues, 4 Xie and colleagues, 5 Wei and Jensen, 6 Sinha and He, 7 and many others. As another major tool, statistical modeling approaches are devoted to comprehensively describing the probabilistic uncertainties linked to those established biological evolution models which may include two topological structures: parallel and phylogenic models (see Figure 1).
The joint parallel evolution process probability Pr(Ancestor, Species 1, 2, and 3) is Pr(Ancestor) ∏ Jukes and Cantor 8 proposed the first probabilistic nucleotide evolution model which assumes substitution to take place randomly among four types of nucleotides "A[1]T [2]C [3]G [4]". The transition (from nucleotide i to nucleotide j) probability up to time t is derived as We assume that the substitution rate parameter (α) is constant for different species and the evolution duration (t) is represented by specific time period (t 0 , t 1 , t 2 , or t 3 ) for the associated divergence process (see Figure 1). Our question is how to effectively estimate the marginal probability for the given orthologous species sequence set without knowing the genotype of the ancestor and/or root.

Material and methods
For the unknown ancestor sequence, we simply assume that the nucleotide on any site follows a tetranomial distribution with categories {ATCG} and equal proportion, (1/4). We further assume that each nucleotide on the ancestor sequence evolves independently (under the probability law, Eq. (3)), so that each species sequence is a series of nucleotides which follow another tetranomial distribution identically and independently. The state space is {ATCG} and the state proportions are (P A , P T , P C , P G ) which can be calculated by Thus, each species' nucleotide follows the same tetranomial distribution as the ancestor nucleotide. Under the independence assumption, the probability for the species sequence is simply a product of all nucleotide marginal probabilities, (1/4). This formulation can also be used to sample the unknown ancestor state among {ATCG} given the reference species state j = (A, T, C or G), since the posterior distribution among {ATCG} for the unknown ancestor state can be easily derived to be We now briefly investigate the ambiguity extent to which different sources of sequence are aligned. For simplicity, we use the Jukes and Cantor 8 model and assume the ancestor vs species nucleotide identity ("ancestor = species") probability is which equals p t ii ( ), in Eq. (3), the substitution probability is p, which equals p t ij ( ) in Eq. (3) for i j ≠ . The identity probability between two species ("species = species") nucleotides with equal evolution duration is thus The statistical sequence evolution model works on probabilistic transition from the ancestor nucleotide to species nucleotide. Since the ancestor sequence is never known for a direct alignment, we may sample it a posteriori given the reference species nucleotide. The probability for event "X", nucleotide identity between such a posterior ancestor nucleotide and another species' nucleotide ("posterior ancestor = species") other than the reference species nucleotide, is derived as The three identity probabilities (Eqs. (6), (7), and (8)) are plotted in Figure 2 where the ancestor-to-species transition probability (p) varies. We find that, the identity probabilities for these three types of matched nucleotides follow the order Pr(ancestor = species)  Pr(species = species)  Pr(posterior ancestor = species), or (6)  (7)  (8). On calculating the probability of a set of orthologous sequences Dovepress submit your manuscript | www.dovepress.com

Dovepress
Thus, alignment between the ancestor and species sequences may be less ambiguous than between two species sequences. Figure 2 also indicates that, the difference among these three types of identity probability is most significant in the middle of interval (0, 1/4). The dominance of "Pr(ancestor = species) over Pr (reference species = another species)" is more significant than that of "Pr (reference species = another species) over Pr (posterior ancestor = another species)", and the difference between the lower two curves ( Figure 2) seems not to be relatively large. Thus pairwise alignment between the posterior ancestor and another species sequence may achieve similar unambiguity as alignment between species. Now we study a multispecies orthologous sequence set (say Human = Species 1, Mouse = Species 2, and Rat = Species 3 in Figure 1). We denote by B a , B 1 , B 2 , and B 3 the sequences of the unknown ancestor, Species 1, 2, and 3. Under the nucleotide substitution model and unambiguous matching, the probability for the set of sequences under parallel evolution is where π ( ) ⋅ is the identical and independent tetranomial distribution for the ancestor nucleotide with state space {ATCG} and equal (1/4) proportion. The result is obtained by integrating out four possible ancestor nucleotides on each site for a marginal nucleotide group (three members across species) probability and multiplying these individual marginal probabilities along the sequence. Similarly, the phylogenic evolution model requires integrating out both the ancestor and root nucleotides on each site to get the result. Note that multiple alignment is not needed under the substitution model since no gaps are allowed. For general nucleotide substitution-insertion-deletion, the probabilistic evolution model developed by Rivas 9 gives the overall "substitution, insertion, and deletion" probabilities from the ancestor to species given divergence time. Calculating the evolution probability from the ancestor (assumed to be known) to the observed species sequence using the Rivas 9 model may require multiple alignments up front in order to match those nucleotides between the ancestor and species. Aside from not knowing the ancestor sequence, unambiguous alignment may not exist due to moderate sequence divergence. 7 Thus, one can underestimate the sequence set probability which is induced in a similar way to Eq. (9) because it simply picks one alignment to calculate the sequence set probability without incorporating other possible alignments. Ignoring ambiguous alignment may also lead to incorrect phylogenic inference and/or misleading sequence taxa partition pattern. 10,11 Under a Dovepress moderate sequence length (∼100 nucleotides), a tetranomial distribution for each nucleotide along the ancestor sequence may be used to sample ancestor nucleotides independently to form a large number (n) of sequences which are further used to induce a probabilistically evolved set of species sequences based on the Rivas 9 model. The sequence set probability is simply estimated as the number of exact duplicates of the given sequence set divided by n. However, this is highly impractical under moderate sequence length due to the small chance of duplicate sequence sets. Another possible way is applying pairwise alignment between each sampled ancestor sequence and observed species-specific sequences, and the sequence set probability may be done by averaging these evolution probabilities over all sampled ancestor sequences. However, this may also be inefficient due to noninformative ancestor sampling and lack of reliable alignment between a random sequence and species sequences. Thus it becomes desirable to propose and investigate more efficient multiple-imputation-like approaches such as using posterior ancestor samples which may offer multiple representative alignment results conditional on a reference species sequences for sequence set probability elicitation. Instead of following the theme of Eq. (9), we turn to calculating P (B 1 , where [II] is obtained after multiplicity over all nucleotide marginal probabilities for Species 1 (see Eq. (4)). As for [I], since B a (posterior ancestor) is sampled from the reference sequence (offspring) B 1 and the integrand is the offspring (B 2 , B 3 ) probability derived from the representative ancestor B a which is already linked to offspring B 1 through posterior sampling. Monte Carlo integration introduced in Eq. (10) realistically implements the joint probability of multiple post-evolution sequences by working on pairwise alignments between the sampled ancestor sequence and observed species-specific sequences. Under phylogenetic (tree-structured) evolution (the right panel in Figure 1), the sequence set probability can be written as where B 1 is Species 1 (Human) sequence, B r is the root sequence, and B a is the ancestor sequence in Figure 1. However, if we use Species 2 (Mouse) sequence B 2 as the reference sequence, then we have the following decomposition Note that, As in Eq. (4), we assign 1/4 to the probability for each nucleotide along the reference sequence B 2 after applying Eq. (13). Only pairwise alignment between the posterior ancestor sequence and species sequence is used for Monte Carlo integration (Eqs. (10), (11), and (12)). Since the probability of a sequence evolving from an ancestor is obtained by multiplying over all individual nucleotide evolution probabilities along a sequence, a large sequence length (say 100) may result in an overly small probability and lead to numerical overflow. The log-probability (LogPr) for a species On calculating the probability of a set of orthologous sequences Dovepress submit your manuscript | www.dovepress.com Dovepress evolutionary sequence from the ancestor is the summation of individual nucleotide evolutionary log-probabilities, and the evolutionary probability expectation obtained from Monte Carlo integration (exp(LogPr) mean), can be implemented by using moment generating function with argument one. Normality of these randomly produced LogPrs leads to the simple result of exp (µ + σ 2 /2) where µ and σ 2 are the sample mean and variance for these LogPrs.

simulation and real data study
We first introduce in detail the extended Jukes and Cantor model by Rivas 9 which will be used for our simulation study. The transition probabilities among general states {-ATCG} ("-" is the gap or covalent bond between two nucleotides) until time t are

s t r t t s t s t r t s t t s t s t s t r t
where the {-ATCG} column to the left of the transition probability matrix represents the initial (ancestor) states and the {-ATCG} row on top of the matrix represents the final (species) states. Specifically,

Parallel evolution model
We refer to the left panel of Figure 1. 1. Simulate the ancestor sequence with length = L 0 ; 2. Simulate species "1, 2, 3" sequences from this simulated ancestor sequence; 3. Apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 2 and 3 conditional on Species 1 sequence, where the unknown ancestor sequences are sampled using Eq. (16) with corresponding divergence time; 4. As a numerical verification, we apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 1 and 3 conditional on Species 2 sequence, where the unknown ancestor sequences are sampled using Eq. (16) with corresponding divergence time; 5. As another numerical verification, we apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 1 and 2 conditional on Species 3 sequence, where the unknown ancestor sequences are sampled using Eq. (16) with corresponding divergence time; 6. We investigate the consistency among different references. 7. Various divergence time vector (t 1 , t 2 , t 3 ) in the left panel of Figure 1 and transition parameter (β and q 0 in Eq. (15)) configurations are given in Table 1, where transition parameter (α in Eq. (15)) is standardized into one unit.

Phylogenic evolution model
We refer to the right panel of Figure 1. 1. Simulate root sequence with length = L 0 and the evolved ancestor sequence for Species 2 and 3; 2. Simulate Species 1 sequence from this simulated root sequence, and simulate the Species 2 and 3 sequences from this simulated ancestor sequence; 3. Apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 2 and 3 conditional on Species 1 sequence, where the unknown root and ancestor sequences are sampled using Eq. (16) with corresponding divergence times; 4. As a numerical verification, we apply Monte Carlo integration to randomly produced log(evolution probabilities) for Carlo integration to randomly produced log(evolution probabilities) for Species 1 and 2 conditional on Species 3 sequence, where the unknown root and ancestor sequences are sampled using Eq. (16) with corresponding divergence times; 6. We investigate the consistency among different references. 7. Various divergence time vector (t 0 , t 1 , t 2 , t 3 ) in the right panel of Figure 1 and transition parameter (β and q 0 in Eq. (15)) configurations are given in Table 1, where transition parameter (α in Eq. (15)) is standardized into one unit. We use the same transition parameter (β and q 0 in Eq. (15)) from the parallel model simulation and make evolution divergence times comparable between the parallel and phylogenic models. We collect LogPrs from 1000 Monte Carlo simulations. The distribution of these LogPrs are plotted in Figures 3, 4, 5 and 6. A Kolmogorov-Smirnov normality test gives p-value (0.15) for all LogPr sets, which means that the difference between the produced LogPrs and a normally distributed random variable is not significant. The normality assumption for LogPrs holds and the probability approximation based on log-normal distribution is reasonable. For such an assumption, a heuristic justification without rigorous theoretical proof is as follows: Given each randomly produced ancestor sequence, each nucleotide (event) LogPr on the non-reference sequences acts as an independently and identically distributed random variable, and the summation of these LogPrs follows the central limit theorem for a large sample size (sequence length). By increasing the ancestor or root length from 100 to 500, we can see that the relationship between LogPr and the sequence length is approximately linear. Another observation from Tables 2 and 3 is that different reference species sequences may lead to inconsistent sequence set probabilities due to different evolution durations and/or topological locations within the phylogenic structure. The phylogenic evolution model (the right panel of Figure 1) seems to show more inconsistency than the parallel evolution model (the left panel of Figure 1) does due to the dual missing sequences (the root and ancestor) instead of ancestor only in the parallel evolution model. A reference sequence which is closer to the root and/or ancestor is preferable since the imputed multiple roots and/or ancestors tend to be more informative due to shorter divergence. The CLUSTAL W multiple alignment 12 -induced probability is obtained by moving along the sequence set which holds nucleotides {ATCG} and possible gaps (covalent bonds) and applying Rivas 9 model, where "one gap with two nucleotides" across three matched sequence sites stands for a deletion and "two gaps with one nucleotide" across three matched sequence sites stands for an insertion. For each simulated sequence set, the discrepancy between Monte Carlo integration (MCI) and single multiple alignment (MA) induced probabilities are clearly more significant than that among probabilities estimated from different reference species sequences.

CREB promoters study
From the ABS database, 13 we extracted the promoter regions of transcription factor CREB for three mammals (human,    On calculating the probability of a set of orthologous sequences Dovepress submit your manuscript | www.dovepress.com Dovepress mouse and rat). We used MEGA 4.1 package 14 to construct the phylogeny tree with corresponding divergence times under uniform transition rate 1 (see Figure 7). These are used for sampling the posterior ancestor and root.
Since no current packages give us β and q 0 maximum likelihood estimation for Rivas 9 model (Eq. (14)), we mainly investigate the sequence set probability sensitivity to β and q 0 input by trying different values. We report the means and variances as well as estimated Log(sequence set probabilities) under two parameter settings for β and q 0 under different reference species. We also report the p-values from Kolmogorov-Smirnov normality test ( Table 4). The normality test results are sensitive to parameter input and reference species selection, which may be due to the fact that conservation levels/transition probabilities are likely to be nonhomogeneous along the sequences. Two LogPr distributions under associated parameter inputs are plotted in Figure 8. As a verification, we apply multiple alignment to these three promoters and at each site we calculate the nucleotide identity proportion within the window (with size 23) starting from this site at the gene direction ( Figure 9). The conservation levels show that evolutionary transition rates  Dovepress are approximately constant piece-wisely, thus the central limit theorem discussed for the simulated data study may still apply to these LogPrs on each promoter segment with quasi-constant conservation level under certain nucleotide insertion-deletion parameter (β, q 0 ) values.

Discussion
We proposed and investigated some promising numerical algorithms for accurately estimating the probability of a set of orthologous sequences with equal length under certain assumptions. Our approach was to informatively shuffle the unknown ancestors and/or roots and to find the distributional characteristics of simulated log-probabilities in order to reasonably approximate the true probability. The merit of our approach depends on how well the ancestor and/or root is imputed based on certain pentanomial distribution proportions (p -, p A , p T , p C , p G ) in Eq. (16) using the evolution model 9 and how reliably the pairwise Needleman-Wunsch alignment is applied to cross-species matching of nucleotides which are supposed to come from the same ancestor entry {-, A, T, C or G}. The former depends on the divergence duration from the ancestor/root to the reference sequence and the latter may depend on the species-specific adjustment of pairwise alignments based on phylogenic information.
When this piece of information is not immediately available, the algorithms by Yang, 15 Redelings and Suchard, 11 and MEGA package 14 are useful. Recently, Wong and colleagues 16 demonstrated that various alignments may lead to quite inconsistent inference. Although distance estimation for multiple species from a common ancestor may lack some accuracy using only one sequence set (Figure 9), we used MEGA package for phylogenic structure information for real sequence set probability estimation. Note that we only use background sequences as examples to demonstrate our algorithms by assuming independent tetranomial distribution among {ATCG} along sequences. For the set of orthologous sequences involving many species (3), we follow the evolutionary process (described by a phylogenic tree) to sample the internal nodes within the phylogenic tree conditional on one selected reference sequence (a terminal node on the phylogenic tree) and apply Monte Carlo integration to these imputed internal nodes for obtaining LogPrs (we omit the details). As one referee points out, it may be unreliable to directly apply our algorithms to sequences with very irregular lengths, since the insertion-deletion events need to be identified by matching nucleotides across all involved species other than due to artificial sequence truncation. Thus a crude multiple alignment across such

Publish your work in this journal
Submit your manuscript here: http://www.dovepress.com/advances-and-applications-in-bioinformatics-and-chemistry-journal Advances and Applications in Bioinformatics and Chemistry is an international, peer-reviewed open-access journal that publishes articles in the following fields: Computational biomodelling; Bioinformatics; Computational genomics; Molecular modelling; Protein structure modelling and structural genomics; Systems Biology; Computational Biochemistry; Computational Biophysics; Chemoinformatics and Drug Design; In silico ADME/Tox prediction. The manuscript management system is completely online and includes a very quick and fair peerreview system, which is all easy to use. Visit http://www.dovepress.com/ testimonials.php to read real quotes from published authors.

48
Liu et al Dovepress submit your manuscript | www.dovepress.com Dovepress Dovepress sequences may overly produce insertion and/or deletions. A rough solution may involve first applying multiple alignment procedures to these sequences and then segmenting the aligned sequences into subsequences involving different numbers of species followed by segment-wise Monte Carlo integration. However, the internal edge-effects introduced by segmentation deserves further study. Lastly, we highlight that applying the proposed algorithms to real sequences is not so straightforward in view of heterogeneous conservation patterns along the orthologous sequences, which poses as an important future research topic.