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[@b1-aabc-2-037] 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,[@b2-aabc-2-037] Kellis and colleagues,[@b3-aabc-2-037] Moses and colleagues,[@b4-aabc-2-037] Xie and colleagues,[@b5-aabc-2-037] Wei and Jensen,[@b6-aabc-2-037] Sinha and He,[@b7-aabc-2-037] 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](#f1-aabc-2-037){ref-type="fig"}).

The joint parallel evolution process probability Pr(Ancestor, Species 1, 2, and 3) is

Pr

(

Ancestor

)

∏

i

=

1

3

Pr

(

Species 

i

∣

Ancestor

)

,

and the joint phylogenic evolution process probability Pr(Root, Ancestor, Species 1, 2, 3) is

Pr

(

Root

)

Pr

(

Species 

1

∣

Root

)

Pr

(

Ancestor

∣

Root

)

∏

i

=

2

3

Pr

(

Species 

i

∣

Ancestor

)

.

Jukes and Cantor[@b8-aabc-2-037] 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

p

i

j

(

t

)

=

(

1

/

4

)

\+

(

3

/

4

)

e

\-

4

α

t

for 

(

i

=

j

)

,

and

(

1

/

4

)

\-

(

1

/

4

)

e

\-

4

α

t

for 

(

i

≠

j

)

.

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](#f1-aabc-2-037){ref-type="fig"}). 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)](#FD3){ref-type="disp-formula"}), 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

p

j

=

∑

i

=

1

4

(

1

/

4

)

p

i

j

(

t

)

=

1

/

4

,

i

,

j

∈

{

1234

}

=

{

A

T

C

G

}

.

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

p

(

i

∣

j

)

=

p

i

j

(

t

)

/

∑

k

=

1

4

p

k

j

(

t

)

,

i

=

1

,

2

,

3

,

4\.

We now briefly investigate the ambiguity extent to which different sources of sequence are aligned. For simplicity, we use the Jukes and Cantor[@b8-aabc-2-037] model and assume the ancestor vs species nucleotide identity ("ancestor = species") probability is

1

\-

3

p

,

which equals *p~ii~* (*t*), in [Eq. (3)](#FD3){ref-type="disp-formula"}, the substitution probability is *p*, which equals *p~ij~* (*t*) in [Eq. (3)](#FD3){ref-type="disp-formula"} for *i* ≠ *j*. The identity probability between two species ("species = species") nucleotides with equal evolution duration is thus

(

1

\-

3

p

)

2

\+

3

p

2

.

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

Pr

(

reference species

=

another species

)

×

Pr

(

X

∣

reference species

=

another species

)

\+

Pr

(

reference species

≠

another species

)

×

Pr

(

X

∣

reference

species

≠

another species

)

=

(

(

1

\-

3

p

)

2

\+

3

p

2

)

(

1

\-

3

p

)

\+

(

1

\-

(

1

\-

3

p

)

2

\-

3

p

2

)

p

=

\-

48

p

3

\+

36

p

2

\-

9

p

\+

1\.

The three identity probabilities ([Eqs. (6)](#FD6){ref-type="disp-formula"}, [(7)](#FD7){ref-type="disp-formula"}, and [(8)](#FD8){ref-type="disp-formula"}) are plotted in [Figure 2](#f2-aabc-2-037){ref-type="fig"} 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

)

.

Thus, alignment between the ancestor and species sequences may be less ambiguous than between two species sequences. [Figure 2](#f2-aabc-2-037){ref-type="fig"} 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](#f2-aabc-2-037){ref-type="fig"}) 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](#f1-aabc-2-037){ref-type="fig"}). 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

P

(

B

1

,

B

2

,

B

3

)

=

∑

B

a

P

(

B

1

,

B

2

,

B

3

∣

B

a

)

π

(

B

a

)

=

∑

B

a

P

(

B

1

∣

B

a

)

P

(

B

2

∣

B

a

)

×

P

(

B

3

∣

B

a

)

π

(

B

a

)

,

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[@b9-aabc-2-037] 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[@b9-aabc-2-037] 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.[@b7-aabc-2-037] Thus, one can underestimate the sequence set probability which is induced in a similar way to [Eq. (9)](#FD11){ref-type="disp-formula"} 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.[@b10-aabc-2-037],[@b11-aabc-2-037] Under a 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[@b9-aabc-2-037] 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 non-informative 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)](#FD11){ref-type="disp-formula"}, we turn to calculating *P* (*B*~1~, *B*~2~, *B*~3~) in an alternative way under parallel evolution

P

(

B

1

,

B

2

,

B

3

)

=

P

(

B

2

,

B

3

∣

B

1

)

P

(

B

1

)

=

\[

∑

B

a

P

(

B

2

,

B

3

∣

B

a

)

P

(

B

a

∣

B

1

)

\]

×

\[

∑

B

a

P

(

B

1

∣

B

a

)

π

(

B

a

)

\]

=

\[

∑

B

a

P

(

B

2

∣

B

a

)

P

(

B

3

∣

B

a

)

P

(

B

a

∣

B

1

)

\]

\[

∑

B

a

P

(

B

1

∣

B

a

)

π

(

B

a

)

\]

=

\[

I

\]

×

\[

I

I

\]

,

where \[II\] is obtained after multiplicity over all nucleotide marginal probabilities for Species 1 (see [Eq. (4)](#FD4){ref-type="disp-formula"}). 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)](#FD12){ref-type="disp-formula"} 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](#f1-aabc-2-037){ref-type="fig"}), the sequence set probability can be written as

P

(

B

1

,

B

2

,

B

3

)

=

P

(

B

2

,

B

3

∣

B

1

)

P

(

B

1

)

=

\[

∑

B

r

,

B

a

P

(

B

2

,

B

3

∣

B

r

,

B

a

)

P

(

B

r

,

B

a

∣

B

1

)

\]

×

\[

∑

B

r

P

(

B

1

∣

B

r

)

π

(

B

r

)

\]

=

\[

∑

B

r

,

B

a

P

(

B

2

,

B

3

∣

B

r

,

B

a

)

P

(

B

a

∣

B

r

)

P

(

B

r

∣

B

1

)

\]

×

\[

∑

B

r

P

(

B

1

∣

B

r

)

π

(

B

r

)

\]

=

\[

∑

B

r

,

B

a

P

(

B

2

,

B

3

∣

B

a

)

P

(

B

a

∣

B

r

)

P

(

B

r

∣

B

1

)

\]

×

\[

∑

B

r

P

(

B

1

∣

B

r

)

π

(

B

r

)

\]

=

\[

∑

B

r

,

B

a

P

(

B

2

∣

B

a

)

P

(

B

3

∣

B

a

)

P

(

B

a

∣

B

r

)

P

(

B

r

∣

B

1

)

\]

×

\[

∑

B

r

P

(

B

1

∣

B

r

)

π

(

B

r

)

\]

=

\[

I

\]

×

\[

I

I

\]

,

where *B*~1~ is Species 1 (Human) sequence, *B~r~* is the root sequence, and *B~a~* is the ancestor sequence in [Figure 1](#f1-aabc-2-037){ref-type="fig"}. However, if we use Species 2 (Mouse) sequence *B*~2~ as the reference sequence, then we have the following decomposition

P

(

B

1

,

B

2

,

B

3

)

=

P

(

B

1

,

B

3

∣

B

2

)

P

(

B

2

)

=

\[

∑

B

r

,

B

a

P

(

B

1

,

B

3

∣

B

r

,

B

a

)

P

(

B

r

,

B

a

∣

B

2

)

\]

×

\[

∑

B

r

P

(

B

2

∣

B

r

)

π

(

B

r

)

\]

=

\[

∑

B

r

,

B

a

P

(

B

1

∣

B

r

)

P

(

B

3

∣

B

a

)

P

(

B

r

,

B

a

∣

B

2

)

\]

×

\[

∑

B

r

P

(

B

2

∣

B

r

)

π

(

B

r

)

\]

Note that,

P

(

B

r

,

B

a

∣

B

2

)

=

P

(

B

r

∣

B

a

)

(

B

a

∣

B

2

)

and

∑

B

r

P

(

B

2

∣

B

r

)

π

(

B

r

)

=

∑

B

r

\[

∑

B

a

P

(

B

2

∣

B

a

)

P

(

B

a

∣

B

r

)

\]

π

(

B

r

)

.

As in [Eq. (4)](#FD4){ref-type="disp-formula"}, we assign 1/4 to the probability for each nucleotide along the reference sequence *B*~2~ after applying [Eq. (13)](#FD16){ref-type="disp-formula"}. Only pairwise alignment between the posterior ancestor sequence and species sequence is used for Monte Carlo integration ([Eqs. (10)](#FD12){ref-type="disp-formula"}, [(11)](#FD13){ref-type="disp-formula"}, and [(12)](#FD14){ref-type="disp-formula"}). 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 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[@b9-aabc-2-037] 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

\[

\-

A

T

C

G

\-

σ

(

t

)

ξ

(

t

)

ξ

(

t

)

ξ

(

t

)

ξ

(

t

)

A

γ

(

t

)

r

(

t

)

s

(

t

)

s

(

t

)

s

(

t

)

T

γ

(

t

)

s

(

t

)

r

(

t

)

s

(

t

)

r

(

t

)

C

γ

(

t

)

s

(

t

)

s

(

t

)

r

(

t

)

s

(

t

)

G

γ

(

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,

{

r

(

t

)

=

(

1

/

4

)

e

\-

β

t

\+

(

3

/

4

)

e

\-

(

4

α

\+

β

)

t

;

S

(

t

)

=

(

1

/

4

)

e

\-

β

t

\-

(

1

/

4

)

e

\-

(

4

α

\+

β

)

t

;

γ

(

t

)

=

1

\-

e

\-

β

t

;

ξ

(

t

)

=

(

1

/

4

)

(

1

\-

q

0

)

e

\-

β

t

;

σ

(

t

)

=

1

\-

(

1

\-

q

0

)

e

\-

β

t

.

For these generalized transition probabilities, we refer to the notations from the substitution model ([Eq. (3)](#FD3){ref-type="disp-formula"}) and denote the element (*u*, *v*) in the matrix ([Eq. (14)](#FD17){ref-type="disp-formula"}) to be *p* ~−1,~ *~v~*~−1~(*t*), where *u* (row index) and *v* (column index) = 1, 2, 3, 4, 5. Parameter 0 \< *q*~0~ ≤ 1 controls the background (nongap) frequency at time *t*. Specifically, letting *β* = 0 leads to the original Jukes and Cantor model ([Eq. (3)](#FD3){ref-type="disp-formula"}) and *q*~0~ = 1 excludes nucleotide insertion. Since each pair of neighboring ancestor nucleotides holds a potential insertion site (gap, "--") with an overall "gap:nongap" ratio of one, we assume a pentanomial distribution for general ancestor nucleotide states with sample space {--\[0\]*A*\[1\]*T*\[2\]*C*\[3\]*G*\[4\]} and normalized probability set (*p*~0~ = 1/2, *p*~1~ = *p*~2~ = *p*~3~ = *p*~4~ = 1/8). This assumption is useful for sampling the posterior ancestor state among {--*ATCG*} given the reference species state (--, *A*, *T*, *C* or *G*). If we denote the general species nucleotide state to be *J* ε {--*ATCG*}, then the posterior distribution among {--*ATCG*} for the unknown ancestor state is

P

(

i

∣

j

)

=

p

i

p

i

j

(

t

)

/

∑

k

=

0

4

p

k

p

k

j

(

t

)

,

i

=

0

,

1

,

2

,

3

,

4\.

Parallel evolution model
------------------------

We refer to the left panel of [Figure 1](#f1-aabc-2-037){ref-type="fig"}.

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)](#FD19){ref-type="disp-formula"} 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)](#FD19){ref-type="disp-formula"} 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)](#FD19){ref-type="disp-formula"} 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](#f1-aabc-2-037){ref-type="fig"} and transition parameter (*β* and *q*~0~ in [Eq. (15)](#FD18){ref-type="disp-formula"}) configurations are given in [Table 1](#t1-aabc-2-037){ref-type="table"}, where transition parameter (*α* in [Eq. (15)](#FD18){ref-type="disp-formula"}) is standardized into one unit.

Phylogenic evolution model
==========================

We refer to the right panel of [Figure 1](#f1-aabc-2-037){ref-type="fig"}.

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)](#FD19){ref-type="disp-formula"} with corresponding divergence times;

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 root and ancestor sequences are sampled using [Eq. (16)](#FD19){ref-type="disp-formula"} with corresponding divergence times;

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 root and ancestor sequences are sampled using [Eq. (16)](#FD19){ref-type="disp-formula"} 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](#f1-aabc-2-037){ref-type="fig"} and transition parameter (*β* and *q*~0~ in [Eq. (15)](#FD18){ref-type="disp-formula"}) configurations are given in [Table 1](#t1-aabc-2-037){ref-type="table"}, where transition parameter (α in [Eq. (15)](#FD18){ref-type="disp-formula"}) is standardized into one unit. We use the same transition parameter (*β* and *q*~0~ in [Eq. (15)](#FD18){ref-type="disp-formula"}) 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](#f3-aabc-2-037){ref-type="fig"}, [4](#f4-aabc-2-037){ref-type="fig"}, [5](#f5-aabc-2-037){ref-type="fig"} and [6](#f6-aabc-2-037){ref-type="fig"}. 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](#t2-aabc-2-037){ref-type="table"} and [3](#t3-aabc-2-037){ref-type="table"} 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](#f1-aabc-2-037){ref-type="fig"}) seems to show more inconsistency than the parallel evolution model (the left panel of [Figure 1](#f1-aabc-2-037){ref-type="fig"}) 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[@b12-aabc-2-037]- induced probability is obtained by moving along the sequence set which holds nucleotides {*ATCG*} and possible gaps (covalent bonds) and applying Rivas[@b9-aabc-2-037] 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,[@b13-aabc-2-037] we extracted the promoter regions of transcription factor *CREB* for three mammals (human, mouse and rat). We used MEGA 4.1 package[@b14-aabc-2-037] to construct the phylogeny tree with corresponding divergence times under uniform transition rate 1 (see [Figure 7](#f7-aabc-2-037){ref-type="fig"}). These are used for sampling the posterior ancestor and root.

Since no current packages give us *β* and *q*~0~ maximum likelihood estimation for Rivas[@b9-aabc-2-037] model ([Eq. (14)](#FD17){ref-type="disp-formula"}), 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](#t4-aabc-2-037){ref-type="table"}). 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](#f8-aabc-2-037){ref-type="fig"}. 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](#f9-aabc-2-037){ref-type="fig"}). The conservation levels show that evolutionary transition rates 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)](#FD19){ref-type="disp-formula"} using the evolution model[@b9-aabc-2-037] 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,[@b15-aabc-2-037] Redelings and Suchard,[@b11-aabc-2-037] and MEGA package[@b14-aabc-2-037] are useful. Recently, Wong and colleagues[@b16-aabc-2-037] 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](#f9-aabc-2-037){ref-type="fig"}), 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 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.

We thank Terence P Speed for his directions on evolution models when he visited Yale Center for Statistical Genomics and Proteomics (YCSGP) in May 2004. We are also grateful to Stéphane Robin and many anonymous referees for their constructive and insightful comments which greatly improved our work.

**Disclosure**

The authors report no conflicts of interest in this work.

![Evolution models.](aabc-2-037f1){#f1-aabc-2-037}

![Identity probabilities for matched nucleotide pair ("ancestor vs species," "species vs species" and "posterior ancestor vs species").](aabc-2-037f2){#f2-aabc-2-037}

![Log(probability) densities from parallel model 1--6.\
**Notes:** Solid line, Species 1; dashed line, Species 2; dotted line, Species 3.](aabc-2-037f3){#f3-aabc-2-037}

![Log(probability) densities from parallel model 7--12.\
**Notes:** Solid line, Species 1; dashed line, Species 2; dotted line, Species 3.](aabc-2-037f4){#f4-aabc-2-037}

![Log(probability) densities from phylogenic model 1--6.\
**Notes:** Solid line, Species 1; dashed line, Species 2; dotted line, Species 3.](aabc-2-037f5){#f5-aabc-2-037}

![Log(probability) densities from phylogenic model 7--12.\
**Notes:** Solid line, Species 1; dashed line, Species 2; dotted line, Species 3.](aabc-2-037f6){#f6-aabc-2-037}

![Phylogeny tree for orthologous *CREB* promoters.](aabc-2-037f7){#f7-aabc-2-037}

![Distribution of LogPrs (*CREB* promoters).](aabc-2-037f8){#f8-aabc-2-037}

![Nucleotide identity proportion along the upstream promoter regions for transcription factor *CREB*.](aabc-2-037f9){#f9-aabc-2-037}

###### 

Simulation configurations (parallel \[PA\] and phylogenic \[PH\] evolution models, α = 1.0)

  No   *β*    *t*~0~      *t*~1~       *t*~2~       *t*~3~       *q*~0~   *L*~0~
  ---- ------ ----------- ------------ ------------ ------------ -------- --------
              (PA, PH)    (PA, PH)     (PA, PH)     (PA, PH)              
  1    0.05   N/A, 0.05   0.10, 0.10   0.10, 0.05   0.10, 0.05   0.95     100
  2    0.05   N/A, 0.05   0.10, 0.10   0.20, 0.15   0.10, 0.05   0.95     100
  3    0.05   N/A, 0.10   0.20, 0.20   0.20, 0.10   0.20, 0.10   0.95     100
  4    0.20   N/A, 0.15   0.30, 0.30   0.30, 0.15   0.30, 0.15   0.90     100
  5    0.20   N/A, 0.15   0.30, 0.30   0.60, 0.45   0.30, 0.15   0.90     100
  6    0.20   N/A, 0.01   0.05, 0.05   0.05, 0.04   0.05, 0.04   0.80     100
  7    0.05   N/A, 0.05   0.10, 0.10   0.10, 0.05   0.10, 0.05   0.95     500
  8    0.05   N/A, 0.05   0.10, 0.10   0.20, 0.15   0.10, 0.05   0.95     500
  9    0.05   N/A, 0.10   0.20, 0.20   0.20, 0.10   0.20, 0.10   0.95     500
  10   0.20   N/A, 0.15   0.30, 0.30   0.30, 0.15   0.30, 0.15   0.90     500
  11   0.20   N/A, 0.15   0.30, 0.30   0.60, 0.45   0.30, 0.15   0.90     500
  12   0.20   N/A, 0.01   0.05, 0.05   0.05, 0.04   0.05, 0.04   0.80     500

###### 

Computation results (parallel evolution models)

  No   Mean (Log(probability))          0.5 variance (Log(probability))   Log(probability) (MCI)   Length (MA)   Log(probability) (MA)
  ---- -------------------------------- --------------------------------- ------------------------ ------------- -----------------------
       Reference species: 1, 2, and 3                                                                            
  1    −544, −528, −533                 267, 257, 263                     −277, −270, −270         104           −1343
  2    −545, −567, −538                 229, 252, 231                     −317, −315, −307         104           −1343
  3    −572, −565, −568                 197, 205, 212                     −376, −361, −356         109           −1344
  4    −538, −543, −536                 92, 99, 87                        −447, −444, −450         112           −1143
  5    −539, −528, −537                 88, 100, 91                       −450, −429, −446         108           −1118
  6    −565, −563, −556                 152, 130, 119                     −413, −432, −437         123           −1056
  7    −2767, −2727, −2734              1135, 1130, 1121                  −1632, −1597, −1613      536           −6825
  8    −2775, −2868, −2766              1026, 970, 853                    −1748, −1898, −1913      535           −6823
  9    −2835, −2815, −2822              1008, 937, 853                    −1827, −1879, −1969      533           −6801
  10   −2684, −2683, −2673              434, 459, 412                     −2249, −2225, −2261      533           −5786
  11   −2700, −2620, −2685              455, 451, 429                     −2245, −2170, −2257      538           −5678
  12   −2940, −2878, −2860              763, 535, 580                     −2177, −2342, −2280      630           −5413

**Abbreviations:** MCI, Monte Carlo integration; MA, multiple alignment.

###### 

Computation results (phylogenic evolution model)

  No   Mean(Log(probability))           0.5 variance (Log(probability))   Log(probability) (MCI)
  ---- -------------------------------- --------------------------------- ------------------------
       Reference species: 1, 2, and 3                                     
  1    −596, −506, −496                 402, 219, 218                     −193, −287, −279
  2    −604, −554, −511                 378, 253, 180                     −227, −301, −331
  3    −607, −553, −543                 302, 215, 195                     −305, −338, −348
  4    −557, −545, −546                 140, 93, 99                       −417, −452, −447
  5    −562, −550, −556                 119, 100, 95                      −443, −450, −461
  6    −628, −597, −593                 229, 134, 136                     −399, −464, −456
  7    −3063, −2696, −2697              1425, 1029, 1016                  −1638, −1667, −1680
  8    −3049, −2923, −2763              1041, 1074, 888                   −2008, −1849, −1875
  9    −3002, −2904, −2884              1072, 992, 937                    −1930, −1912, −1948
  10   −2782, −2768, −2766              545, 420, 410                     −2237, −2348, −2356
  11   −2804, −2745, −2801              491, 437, 439                     −2314, −2308, −2362
  12   −3215, −3087, −3077              928, 585, 598                     −2287, −2502, −2479

**Abbreviation:** MCI, Monte Carlo integration.

###### 

Computation results for *CREB* promoter sequence set (phylogenic evolution models)

  *β, q~0~*                        Mean (Log(probability))   0.5 variance (Log(probability))   Log(probability) (MCI)   *p*-value
  -------------------------------- ------------------------- --------------------------------- ------------------------ ------------------------
  Reference species: 1, 2, and 3                                                                                        
  0.2, 0.7                         −1912, −1408, −1462       5724, 249, 656                    3812, −1159, −806        = 0.11, \<0.01, \<0.01
  0.3, 0.7                         −1964, −1403, −1467       5015, 259, 777                    3051, −1144, −690        \>0.15, \<0.01, \<0.01
  0.4, 0.7                         −1998, −1401, −1467       4506, 307, 958                    2508, −1094, −518        = 0.04, \<0.01, \<0.01
  0.5, 0.7                         −2017, −1400, −1487       4501, 345, 1049                   2484, −1055, −438        \<0.01, \<0.01, \>0.15
  0.2, 0.8                         −1789, −1361, −1424       5231, 742, 927                    3442, −619, −497         \>0.15, \<0.01, \<0.01
  0.3, 0.8                         −1833, −1358, −1432       4892, 781, 1056                   3059, −577, −376         \>0.15, \<0.01, = 0.04
  0.4, 0.8                         −1816, −1359, −1446       4240, 822, 1207                   2424, −537, −239         \>0.15, \<0.01, \>0.15
  0.5, 0.8                         −1880, −1361, −1460       4334, 827, 1251                   2454, −534, −209         \>0.15, \<0.01, \>0.15

**Abbreviation:** MCI, Monte Carlo integration.
