• 沒有找到結果。

Confidence intervals for the substitution number in the nucleotide substitution models

N/A
N/A
Protected

Academic year: 2021

Share "Confidence intervals for the substitution number in the nucleotide substitution models"

Copied!
8
0
0

加載中.... (立即查看全文)

全文

(1)

Confidence intervals for the substitution number in the nucleotide

substitution models

Hsiuying Wang

Institute of Statistics, National Chiao Tung University, Hsinchu, Taiwan

a r t i c l e

i n f o

Article history:

Received 5 February 2011 Revised 5 May 2011 Accepted 18 May 2011 Available online 26 May 2011 Keywords: One-parameter model Two-parameter model Confidence interval Coverage probability Binomial distribution Substitution rate

a b s t r a c t

In the nucleotide substitution model for molecular evolution, a major task in the exploration of an evolutionary process is to estimate the substitution number per site of a protein or DNA sequence. The usual estimators are based on the observation of the difference proportion of the two nucleotide sequences. However, a more objective approach is to report a confidence interval with precision rather than only providing point estimators. The conventional confidence intervals used in the literature for the substitution number are constructed by the normal approximation. The performance and construc-tion of confidence intervals for evoluconstruc-tionary models have not been much investigated in the literature. In this article, the performance of these conventional confidence intervals for one-parameter and two-parameter models are explored. Results show that the coverage probabilities of these intervals are unsatisfactory when the true substitution number is small. Since the substitution number may be small in many situations for an evolutionary process, the conventional confidence interval cannot provide accurate information for these cases. Improved confidence intervals for the one-parameter model with desirable coverage probability are proposed in this article. A numerical calculation shows the substantial improvement of the new confidence intervals over the conventional confidence intervals.

Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction

A basic process in the evolution of DNA sequences is the substitution of one nucleotide for another during evolution. But since the substitution of one allele for another in a population gen-erally takes thousands of years or longer to complete, the process cannot be directly observed. Thus, to detect evolutionary changes in a DNA sequence, we need to compare two sequences that have descended from a common ancestral sequence. If two sequences of length L differ from each other at X sites, the proportion of differ-ences, X/L, is referred to as the observed or uncorrected divergence. When the degree of divergence between the two sequences com-pared is small, the chance for more than one substitution to have occurred at a site is negligible, and the number of observed differ-ences between the two sequdiffer-ences is close to the actual number of substitutions. However, if the degree of divergence is substantial, the observed number of differences is likely to be smaller than the actual number of substitutions due to multiple hits at the same site. Many methods have been proposed to correct for multiple hits (Holmquist, 1971; Jukes and Cantor, 1969; Kaplan and Risko, 1982; Kimura, 1980, 1981; Lanave et al., 1984). The simplest and most frequently used models are the Jukes and Cantor (1969) one-parameter model and the Kimura (1980) two-parameter model

(Graur and Li, 1999). For a DNA sequence, the Jukes and Cantor one-parameter model assumes that substitutions occur with equal probability, say

a, among the four nucleotide types, A, T, C, G. Since

the time of divergence between two sequences is usually un-known, we cannot estimate

a

directly. Instead, we compute K, the number of substitutions per site since the time of divergence between the two sequences. In the one-parameter model case, K = 2(3at), where 3at is the expected number of substitutions per site in a single lineage.Jukes and Cantor (1969)derived the follow-ing formula: K ¼ 3 4ln 1  4 3p   ð1Þ where p is the probability that the two sequences are different at a site at time t. They proposed the estimator

K1¼  3 4ln 1  4 3^p   ð2Þ to estimate K, where ^p ¼ X=L is the observed proportion of different nucleotides between the two sequences.

The variance of K can be approximated by

VðKÞ ¼ p  p

2

L 1 4

3p

 2:

1055-7903/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.ympev.2011.05.013

E-mail address:wang@stat.nctu.edu.tw

Contents lists available atScienceDirect

Molecular Phylogenetics and Evolution

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / y m p e v

(2)

ByKimura and Ohta (1972), an estimator for the variance of K is VðK1Þ ¼ ^ p  ^p2 L 1 4 3^p  2:

Although the Jukes and Cantor model is a simple model and many substitution models have been constructed in the literature to compete with it, it is still widely-used due to its simplicity and adaptability for many applications (Fu, 1995; Wirgart et al., 1998; Rosenberg, 2005; Chor et al., 2006, etc.).

In the case of the two-parameter model, the differences be-tween two sequences are classified into transitions and transver-sions. Transitions are changes between A and G (purines) or between C and T (pyrimidines). Transversions are changes between a purine and a pyrimidine. The substitute probabilities of transition and transversion are assumed to be different. Let ^p ¼ X1=L and

^

Q ¼ X2=L be the observed proportions of transitional and

transver-sional differences between the two sequences, respectively, where X1and X2are the numbers of transitional and transversional

differ-ences between the two sequdiffer-ences. ByKimura (1980), the number of nucleotide substitutions per site between the two sequences, K2,

is estimated by K2¼ 1 2ln 1 1  2bP  bQ ! þ1 4ln 1 1  2 bQ !

and the sampling variance is approximately given by

VðK2Þ ¼ 1 L bP 1 1  2bP  bQ !2 þ bQ 1 2  4bP  2 bQþ 1 2  4 bQ !2 0 @  bP 1  2bP  bQþ b Q 2  4bP  2 bQþ b Q 2  4 bQ !21 A: ð3Þ

Since the above two variance estimators underestimate the true variances in most mircumstances,Wang et al. (2008)derive im-proved variance estimators using a higher-order Taylor expansion and empirical methods.

The above illustration is under the assumption that the rate of nucleotide substitution is the same for all nucleotide sites. How-ever, this assumption may not hold in some situations because the nucleotide sequences have functional constraints and usually form a secondary structure consisting of loops and stems that have different substitution rates. Kocher and Wilson (1991), Tamura and Nei (1993) and Wakeley (1993, 1994) suggest that the substitution rate varies from site to site according to the gamma distribution for this case.

When the nucleotide substitution at each site follows the Jukes and Cantor model but the substitution rate 3avaries with the gam-ma distributionC(a, b), by Golding (1983) and Nei and Gojobori (1986), the expected number of substitutions per site becomes

H ¼3 4a 1  4 3p  1=a  1 " #

and the variance for the number of the substitutions per site is

VðHÞ ¼pð1  pÞ n 1  4 3p  2ð1=aþ1Þ " #

where a is the shape parameter of the gamma distribution with the density function f(x) = [ba/

C(a)]ebxxa1. Note that H and V(H)

de-pend on only one parameter of the gamma distribution, but not the two parameters a and b because a/b is the mean of the substi-tution rate 3a, and b is a function of a and

a

(Nei and Gojobori, 1986). The estimators H1¼ 3 4a 1  4 3^p  1=a  1 " # and VðH1Þ ¼ ^ pð1  ^pÞ L 1  4 3^p  2ð1=aþ1Þ " #

are used to estimate H and V(H).

The true number of substitutions per site is usually approxi-mated by point estimators. However, from a statistical point of view, a better approach is to report a confidence interval of the substitution number instead of a point estimator for the number of substitutions because the point estimation can only provide a rough estimate without any information about its precision. The confidence interval estimation can quantify the uncertainty associ-ated with the estimate such that we can have the confidence de-gree if the true K or H is belonged to the intervals. In many studies, the confidence intervals are reported associated with the point estimation, e.g. Yang (2007). The conventional confidence interval is constructed using the normal approximation, which can achieve the desirable coverage probability when the sample size is large enough. Here the sample size is the length of a se-quence. However, even when the length is large, from the study shown in Section2, the conventional confidence interval suffers from the serious drawback of unsatisfactory coverage probability when the true substitution number is small. Since in the evolution-ary process of DNA sequences, the true substitution number per site may be very small, the behavior of a confidence interval for the small substitution number case is especially important. Accordingly, the information provided by the conventional confi-dence interval is not very accurate.

In this article, modified confidence intervals for the one-param-eter model with more satisfactory coverage probability are pro-posed in both constant substitution rate and variable substitution rate models. These modified confidence intervals are constructed from a modification approach used in the literature for construct-ing confidence intervals of a binomial proportion.

This article is organized as follows. The coverage probability and expected length of the conventional confidence interval for substitution models with constant substitution rate and variable substitution rate are shown in Section2. Section3gives the pro-posed confidence intervals as well as their performances. Section 4provides an algorithm for selecting a factor such that the cover-age probability of the confidence interval is close to a desirable le-vel. The proposed methods are illustrated by a real data example analyzing the substitution number of owlet-nightjars species in Section5. The article concludes in Section6with a summary. 2. The existing methods

We first consider the constant substitution rate case. A statisti-cal interval ðLð^pÞ; Uð^pÞÞ is said to be a level 1 

c

confidence interval of K if it can cover the true K with at least probability 1 

c, which

is defined as

PKðLð^pÞ < K < Uð^pÞÞ P 1 

c

:

The probability PKðLð^pÞ < K < Uð^pÞÞ is the coverage probability of

the confidence interval, and EKðjUð^pÞ  Lð^pÞjÞ is the expected length

of the confidence interval under the true substitution number K. Usually the coverage probability of a level 1 

c

confidence interval we constructed based on the normal approximation may not be equal to 1 

c

and may be close to the nominal level 1 

c

only when the sample size is large. Furthermore, its coverage probability

(3)

could be far away from the nominal level 1 

c

for a fixed sample size. We usually evaluate confidence intervals in terms of their cov-erage probabilities. A confidence interval with covcov-erage probability close to the nominal level is regarded as better than a confidence with coverage probability far away from the nominal level.

In the Jukes and Cantor one-parameter model, when the substi-tution rate is assumed to be the same for all sites, the conventional confidence interval is constructed by the normal approximation that the statistic

K1 K

ffiffiffiffiffiffiffiffiffiffiffiffiffi

VðK1Þ

p ð4Þ

has an asymptotic standard normal distribution. By inverting from the probability P jK1ffiffiffiffiffiffiffiffiffiffiffiffiffi Kj VðK1Þ p <z1c 2 ! ¼ 1 

c

; ð5Þ

we have the level 1 

c

conventional confidence interval

ðK1 z1c2 ffiffiffiffiffiffiffiffiffiffiffiffiffi VðK1Þ p ;K1þ z12c ffiffiffiffiffiffiffiffiffiffiffiffiffi VðK1Þ p Þ; ð6Þ

where z1c/2is the upper 1 

c/2 cutoff point of the standard

nor-mal distribution.

Note that the equality in(5)only holds when L is large enough. Consequently, the coverage probability of(6)is not exactly equal to 1 

c

for a fixed sample size. We explore its coverage probability by a simulation study.

When the substitution rate is assumed to follow the gamma dis-tribution with shape parameter a, a confidence interval con-structed by the usual normal approximation uses the fact that the statistic

H1 H

ffiffiffiffiffiffiffiffiffiffiffiffiffi

VðH1Þ

p ð7Þ

has an asymptotic standard normal distribution. Consequently, we have the level 1 

c

conventional confidence interval

ðH1 z12c ffiffiffiffiffiffiffiffiffiffiffiffiffi VðH1Þ p ;H1þ z12c ffiffiffiffiffiffiffiffiffiffiffiffiffi VðH1Þ p Þ ð8Þ

We conduct a simulation to explore their coverage probability. For the constant substitution rate case, the simulation method is to generate two descendant sequences from an ancestral sequence. First we set a value for

at, say

v

0, then generate descendant

se-quences with probability 1/4 + 3/4e4atthat the nucleotide at a site

in a descendant sequence is the same as that in an ancestral se-quence and with probability 1/4  1/4e4atthat the nucleotide at

a site in a descendant sequence is equal to one of the three other bases from the ancestral sequence. Then we compute the propor-tion of the different nucleotide in these two descendant sequences and use the proportion as ^p to derive K1and V(K1). To derive the

coverage probability of the confidence interval, it is necessary to calculate the proportion if the true K = 2(3at) belongs to the inter-val from a simulation study. We replicate the above process 1000 times in generating the sequences and deriving K1, and then

calcu-late the proportion that the interval based on K1covers the true K.

The proportion is the coverage probability of the confidence interval approximated by the simulation. The simulation for the Kimura two-parameter model is to generate the descendant se-quences from a common origin using the probability setup for the two-parameter model. By using a method similar to the one-parameter model, we can derive the coverage probability for the two-parameter model. The model and simulation approach are re-ferred to Graur and Li (1999), Nei and Kumar (2000) and Yang (2007).

For the variable substitution rate case, the simulation method is that for each site, we generate

a

0= 3a value from a gamma

distribution C(a, b) for each site, then use these

a

0/3 values as

the substitution rate to generate two descendant sequences from an ancestral sequence. Then by a similar argument as the constant substitution rate case, the coverage probability of a confidence interval for H = 2(a/b) can be calculated.

Figs. 1 and 2plot the coverage probability and expected length of the confidence interval(6)corresponding to different K values, which shows that the coverage probability of the conventional confidence interval is much lower than the nominal level when L = 100 and 500 for K 6 0:06.

In the Kimura two-parameter model, by an argument similar to that in the one-parameter model, the 1 

c

level confidence inter-val of K2is ðK2 z1c=2 ffiffiffiffiffiffiffiffiffiffiffiffiffi VðK2Þ p ;K2þ z1c=2 ffiffiffiffiffiffiffiffiffiffiffiffiffi VðK2Þ p Þ: ð9Þ

Fig. 3shows the coverage probabilities of the confidence inter-vals for L = 100 and 500 when the expected substitution number per site is less than 0.11.

The coverage probability, increasing as K increases, of the inter-val(6)for the Jukes Cantor model is much lower than the nominal level 0.95 when the true K is small. Note that the coverage proba-bility is not a very smooth curve of K shown inFig. 1because X is a discrete random variable which leads to the oscillation of the cov-erage probability. We can see that the performance of the conven-tional confidence interval is not satisfactory because its coverage probability cannot reach the nominal level. This is the same as the Kimura two-parameter model, where the coverage probability is much lower than the nominal level 0.95.

In fact, the performance of the conventional confidence interval for the two-parameter model is even worse. For the one-parameter model,Fig. 1shows that the coverage probability increases to 0.95 when L increases. However, fromFig. 3, the deviation of the cover-age probability to the nominal level 0.95 multiplies when L in-creases. This indicates that the confidence interval constructed for the two-parameter model is less satisfactory in estimating the

One−parameter model, L=100 K coverage probability 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.2 0.4 0.6 0.8 1.0 One−parameter model, L=500 K coverage probability

Fig. 1. The coverage probabilities of the level 0.95 confidence intervals(6)when L = 100 and 500 for K 6 0:06.

(4)

true K. This drawback may be due to the bias of the estimator K2.

Accordingly, we mainly focus on constructing improved confidence intervals based on the one-parameter model in this article.

It is worth noting that there are other approaches such as the bootstrap method and jackknife method for constructing

confidence intervals in phylogenetic study (Dopazo, 1994, etc.). However, these methods cannot provide a closed form. In this study, we aim to establish improved closed form inter-vals such that they can be easily implemented in real applications.

3. New confidence interval

It is not surprising that the coverage probability of the conven-tional confidence interval goes to zero as K goes to zero. This phe-nomenon also commonly occurs in a simple model, the binomial distribution (see Agresti and Coull, 1998; Wang, 2007). In the one-parameter and two-parameter models, the parameters of interest are functions of the binomial proportion p. Since there are several approaches in the literature to modify the conventional confidence interval such that the coverage probability for small p can be improved, we extend these approaches to the one-parame-ter model.

Both models for the constant substitution rate and variable sub-stitution rate are considered.

3.1. Constant substitution rate

First we deal with the model that the substitution rate is as-sumed to be the same for each site.

The first method is the extension from the score approach. The score 1 

c

confidence interval for the binomial proportion of a binomial distribution is constructed by inverting from the set fp : jx=n  pj=pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipð1  pÞ=n6z1c=2g. The endpoints of the score confidence interval are the two solutions of p found by solving the equation ðx=n  pÞ2=ðpð1  pÞ=nÞ ¼ z2

1c=2. By a similar

argument, the score interval for K is obtained by inverting from the set fjK1 Kj= ffiffiffiffiffiffiffiffiffiffiffi VðKÞ p 6z1c=2g: ð10Þ

Note that by(1), we have

p ¼ 3=4ð1  e4=3KÞ: ð11Þ

Replacing it in V(K), we have

VðKÞ ¼3ð3 þ 2e

4K=3þ e8K=3Þ

16L : ð12Þ

The score confidence interval can be derived by replacing(12)in (10)and then solving the equation

ðK1 KÞ

2

=VðKÞ ¼ z2

1c=2 ð13Þ

in K. However, the above equation does not have a closed form, which would need to be solved by a numerical calcula-tion. To prevent this disadvantage, which may cause the inconvenience of usage of the interval, we can use the Taylor expansion to approximate the terms e4K/3 and e8K/3 in (12) by

1 + 4K/3 and 1 + 8K/3 because the substitute number should be small. Consequently, the variance can be approximated by K/L. Replacing the approximated variance in (13) and solv-ing the equation, we have the approximated 1 

c

score con-fidence interval 2K1L þ z21c=2 z1c=2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4K1L þ z21c=2 q 2L ; 2K1L þ z21c=2þ z1c=2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4K1L þ z21c=2 q 2L 0 @ 1 A: ð14Þ The second approach is an extension of the Agresti and Coull approach (1998)to modify the confidence interval(6)by replacing X and L by X þ z2

1c=2=2 and L þ z21c=2=2, respectively. Let

One−parameter model, L=100 K expected length 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.000 0.010 0.020 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.000 0 .004 0.008 0.012 One−parameter model, L=500 K expected length

Fig. 2. The expected lengths of the level 0.95 confidence intervals(6)when L = 100 and 500 for K 6 0:06. Two−parameter model, L=100 K coverage probability 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.4 0.5 0.6 0.7 0.8 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.00 0.02 0.04 0.06 0.08 Two−parameter model, L=500 K coverage probability

Fig. 3. The coverage probabilities of the level 0.95 confidence intervals(9)for L = 100 and 500 for K 6 0:11. Here the substitution probability of transition is fixed at 0.05, and of the substitution probability of transversion ranges from 0 to 0.03.

(5)

^ ^ p ¼X þ z 2 1c=2=2 L þ z2 1c=2 ;Kc1¼  3 4ln 1  4 3 ^ ^ p   and VcðK1Þ ¼ ^ ^ p  ^^p2 L 1 4 3 ^ ^ p  2:

The interval is proposed to have the form

ðKc1 z1c=2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi VcðK1Þ q ;Kc1þ z1c=2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi VcðK1Þ q Þ: ð15Þ

This approach can successfully increase the coverage probability when K goes to zero. We call this interval an adjusted confidence interval.

To compare the intervals, we conduct a simulation to explore their coverage probabilities and expected lengths, as shown inFigs. 4 and 5.

The comparison of the two proposed level 0.95 confidence intervals of K for L = 100 and 500 is shown inFigs. 4 and 5for small K. The solid and dashed lines represent the coverage probabilities of level 0.95 score and adjusted confidence intervals in Fig. 4. The coverage probability of a good confidence interval should be very close to 0.95 for any K. The simulation results show that the new confidence intervals substantially improve the coverage prob-ability when K is small, compared with the conventional confi-dence interval. We also conduct a simulation for larger, L such as L = 2000. They have similar performance as the above cases that the coverage probability of the adjusted interval is higher than the score interval.

The comparison of the expected lengths of the three confidence intervals shows the conventional confidence interval has the short-est expected length. In the criterion of evaluating a confidence interval, although we prefer a confidence interval with shorter ex-pected length, the most important thing is to evaluate its coverage probability performance. A conventional interval, with a smaller

expected length, has poor coverage probability because it cannot cover the true K with desirable frequency.

For the performance of the two new intervals, the simulation study shows that the adjusted interval always has higher coverage probability and longer expected length than the score confidence interval. In this case, we would prefer the score confidence interval because its average coverage probability is closer to the nominal level when the true K is small.

When the true K is not small, the performance of the two new intervals is shown inFig. 6. The coverage probability of the ad-justed interval is close to the nominal level 0.95, but the coverage probability of the score interval is lower than 0.95. The perfor-mance of the score interval is worse than the adjusted interval when the true K is not small. We recall the method used to con-struct (14) is to approximate e4K/3 and e8K/3 by 1 + 4K/3 and

1 + 8K/3, respectively, in(12) and (13). The approximation is only appropriate when K is small. When K is not small, to derive a more accurate score interval, a better approximation for e4K/3 is

using a Taylor expansion up to the second order term 1 + 4K/ 3 + (4K/3)2/2 or up to a higher order term. We also can derive the score interval by a numerical method to solve (13). Thus, when we do not have information about the range of the true K, we recommend the adjusted interval, or using a more accurate score interval.

Note that according toWang et al. (2008), there exist better var-iance estimators for the varvar-iance of K in the one-parameter and two-parameter models instead of the variance estimators V(K1)

and V(K2). Consequently, we can replace the usual variances in

(6), (9) by the improved variances proposed in Wang et al. (2008) to construct another confidence interval. However, this method cannot substantially improve the coverage probability when the true K is small from a simulation study. Thus, to improve the coverage probability, we propose the score approach and the adjusted approach to construct new confidence intervals. New confidence interval, L=100

K coverage probability 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.88 0.92 0.96 1.00 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.90 0.94 0.98

New confidence interval, L=500

K

coverage probability

Fig. 4. Solid and dashed lines present the coverage probabilities for the level 0.95 score and adjusted confidence intervals, respectively when L = 100 and L = 500 for K 6 0:06 for the constant substitution rate case.

New confidence interval, L=100

K expected length 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.04 0.06 0.08 0.10 0.12 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.00 0.02 0.04 0.06

New confidence interval, L=500

K

expected length

Fig. 5. Solid and dashed lines present the expected lengths for the level 0.95 score and adjusted confidence intervals, respectively when L = 100 and L = 500 for K 6 0:06 for the constant substitution rate case.

(6)

3.2. Variable substitution rate

We proceed to consider the Juke and Cantor model when the substitution rate is assumed to follow a gamma distribution C(a, b) with a shape parameter a.

For this case, applying the score approach leads to a messy for-mula for the confidence interval of H because the endpoints of the interval are approximated by solving a polynomial equation with degree 3 after simplification. Although it is hard to provide a closed form for the score interval of H, it is feasible to derive the confi-dence interval by numerical calculation.

In this case, the second approach, the adjusted approach, may be simpler and more useful here. The adjusted confidence interval is constructed by replacing H1and V(H1) in(8)as

Hc 1¼ 3 4a 1  4 3 ^ ^ p  1=a  1 " # and VðHc1Þ ¼ ^ ^ p 1  ^^ p L 1  4 3 ^ ^ p  2ð1=aþ1Þ " # ; respectively.

This leads to the confidence interval

ðHc1 z1c2 ffiffiffiffiffiffiffiffiffiffiffiffiffi VðHc1Þ q ;Hc1þ z1c2 ffiffiffiffiffiffiffiffiffiffiffiffiffi VðHc1Þ q Þ: ð16Þ

To compare the performance of the intervals, we conduct a sim-ulation study to compare their coverage probabilities and expected lengths, which are shown inFigs. 7 and 8.

The comparisons of the standard and proposed level 0.95 confi-dence intervals of H for L = 500 and 1000 are shown inFigs. 7 and 8. Here we consider the longer sequence lengths L = 500 and 1000 instead of L = 100 and 500 in the constant substitution rate case

because the estimation for the variable substitution rate case is more unstable than that for the constant substitution rate case. The solid and dashed lines represent the coverage probabilities of level 0.95 standard and adjusted confidence intervals inFig. 7for New confidence interval, L=100

K coverage probability 0.10 0.15 0.20 0.25 0.30 0.80 0.85 0.90 0.95 1.00 0.10 0.15 0.20 0.25 0.30 0.80 0.85 0.90 0.95 1.00

New confidence interval, L=500

K

coverage probability

Fig. 6. Solid and dashed lines present the expected lengths for the level 0.95 score and adjusted confidence intervals, respectively when L = 100 and L = 500 for K 6 0:3 for the constant substitution rate case.

0.70

0.80

0.90

1.00

variable substitution rate, L=500

H coverage probability 0.01 0.02 0.03 0.04 0.05 0.01 0.02 0.03 0.04 0.05 0.80 0.85 0.90 0.95 1.00

variable substitution rate, L=1000

H

coverage probability

Fig. 7. Solid and dashed lines present the coverage probabilities for the level 0.95 standard and adjusted confidence intervals, respectively when L = 500 and L = 1000 for H 6 0:05 for the variable substitution case.

variable substitution rate, L=500

H expected length 0.01 0.02 0.03 0.04 0.05 0.0 0.2 0.4 0.6 0.8 1.0 0.01 0.02 0.03 0.04 0.05 0.00 0.05 0.10 0.15 0.20

variable substitution rate, L=1000

H

expected length

Fig. 8. Solid and dashed lines present the expected lengths for the level 0.95 standard and adjusted confidence intervals, respectively when L = 500 and L = 1000 for H 6 0:05 for the variable substitution case.

(7)

H 6 0:05. The coverage probability of the standard confidence interval is clearly lower than the nominal level, which reveals the same disadvantage as in the constant substitution rate case. The adjusted confidence interval has higher coverage probability. Although it seems higher than the nominal level, we still prefer the adjusted confidence interval because its minimum coverage probability is near the nominal level 0.95. As we expect, the ex-pected length of the adjusted confidence interval is longer than that of the standard confidence interval because the adjusted con-fidence interval has higher coverage probability. The simulation re-sults of L = 500 and L = 1000 reveal that the intervals has better performance when the length L increases.

The average coverage probabilities of the proposed interval and the conventional interval are included inTable 1 to support the advantages of the proposed confidence intervals.Table 1 shows the coverage probabilities of the level 0.95 conventional interval and the level 0.95 adjusted interval for K in (0, 0.06) in the variable substitution model for different lengths. It reveals that the average coverage probability of the adjusted interval is closer to the nom-inal level 0.95 than the conventional interval.

4. Selection of

c

value

Although the coverage probabilities for the new confidence intervals can be substantially improved over the conventional interval when the true K or H is small,Figs. 4–7show that they still cannot reach the nominal level 0.95 for all p. Since the conven-tional and new confidence intervals are constructed based on the normal approximation, when the sample size is not large enough, their performance may not be good enough. To prevent this prob-lem from occurring at the small sample size case, we propose an algorithm that is modified from algorithms inWang (2007, 2009) to select an appropriate

c

values in(14)–(16)such that the mini-mum coverage probability or the average coverage probability can reach a desirable level. The minimum coverage probability de-notes the minimum value of the coverage probability of a confi-dence interval when K ranges over the domain. The average coverage probability denotes the average value of the coverage probability with respect to a prior for K under the Bayesian setup. The details are referred toWang (2009).

For the model with constant substitution rate, we can first con-sider the minimum coverage probability criterion. To construct a modified confidence interval by applying the algorithms inWang (2007, 2009), first we transform the confidence interval of K to the confidence interval of p. By(1), K is a one-to-one function of p, where p can be expressed as(11)in terms of K. Thus, for a level 1 

c

confidence interval (L(X), U(X)) of K, we can make a transfor-mation such that the interval (3/4(1  e4/3L(X)), 3/4(1  e4/3U(X)))

is a level 1 

c

confidence interval of p. If we hope to obtain a confidence interval with a minimum coverage probability 0.95, then we can select an appropriate z1-c/2 by using a procedure

modified fromWang (2007).

The algorithm of the modified procedure inWang (2007)is as follows. Since the length of the sequence is L, the possible values of observed X are {0, 1, . . . , L}. To adopt the procedure, the confidence interval needs to satisfy the assumption that 3/4(1  e4/3L(X)) and 3/4(1  e4/3U(X)) are increasing functions of

X (Wang, 2010). In addition to this increasing condition, from(2), since the log function must define on a positive domain, it needs to require that X/L is less than 3/4. Thus, we have to check if 3/4(1  e4/3L(X)) and 3/4(1  e4/3U(X)) are increasing functions of

X for X 6 3=4L. By checking the three confidence intervals(6, 14, and 15)in the one-parameter model, we found that only the score interval(14)satisfies the assumption. Therefore, we modified the score interval using the following procedure to select a

c

value in the confidence interval(14)such that the confidence interval has minimum coverage probability 0.95.

Algorithm 1. Procedure to derive a factor in a score confidence interval (L(X), U(X)) with the form(14) such that the confidence interval has a minimum coverage probability 0.95.

Step 1. Let (L⁄(X), U(X)) = (3/4(1  e4/3L(X)), 3/4(1  e4/3U(X))) .

Step 2. Calculate the endpoints L(X) and U(X) for x

e

{0, 1, . . . , [3/4L]}, where [w] denotes the largest integer less than w. Then list the endpoints L(X) and U(X) that are greater

than zero and smaller than 1.

Step 3. Calculate the coverage probabilities for p in the set of endpoints of Step 2 which are greater than zero and smaller than 1. The minimum value of these coverage probabilities is the minimum coverage probability of the confidence interval. The details refer toWang (2007, 2009).

Step 4. Find the

c

value such that the minimum coverage prob-ability derived in Step 3 is equal to 0.95. Then the confidence interval (L(X), U(X)) with this

c

value is the confidence interval with a minimum coverage probability 0.95.

The confidence level 0.95 in the algorithm can be replaced by other values depending on the precision we prefer.

For a model with variable substitution rate, we can consider the average coverage probability criterion. An argument similar to that inWang (2009)can be used to construct a confidence interval of H such that it has the desirable average coverage probability with re-spect to a prior.

5. Illustrative example

We use the DNA sequences of the avian family Aegothelidae (commonly known as owlet-nightjars) to illustrate the proposed approaches. Owlet-nightjars are small nocturnal birds related to the nightjars and frogmouths. Most are native to New Guinea, but some species extend to Australia, Moluccas, and New Caledo-nia. There is a single monotypic family Aegothelidae with the genus Aegotheles. The family Aegothelidae comprises only nine ex-tant species, all in a single genus, Aegotheles.

Dumbacher et al. (2003)used mitochondrial DNA sequence to construct a phylogeny of the owlet-nightjars. They analyzing mtDNA sequences cytochrome b and ATPase subunit 8 suggests that there are 11 living species of owlet-nightjar and one that went extinct early in the second millennium AD. The taxon listed in Ta-ble 2ofDumbacher et al. (2003)includes albertisi albertisi, wallacii wallacii, wallacii gigas, etc. The Genbank numbers for the se-quences are AY090664-AY090698 (for cytochrome b) and AY090699-AY090736 (for ATPase 8).

Table 1

The average coverage probabilities of the level 95% conventional interval and adjusted interval for 0 < K < 0.06 in the variable substitution model.

L Conventional interval Adjusted interval

500 0.858 0.985

1000 0.891 0.975

2000 0.913 0.97

Table 2

The level 95% conventional, score and adjusted confidence intervals for sequences AY090696 and AY09069.

Conventional confidence interval (0.0001, 0.0145)

Score confidence interval (0.0028, 0.01875)

(8)

We use the DNA sequences AY090696 and AY090697 of wallacii wallacii and wallacii gigas to illustrate the new confidence inter-vals for the one-parameter model. First we apply the multiple sequences alignment procedure of MEGA software to these sequences, and then count the number of different nucleotides between the sequences. The number of different nucleotides between the two sequences is four among 550 which is the length size. The three intervals for the number of nucleotide substitutions per site between the two sequences are listed inTable 2.

Table 2shows that the lower confidence bounds of the score and adjusted confidence intervals are much larger than those of the conventional confidence interval. This reveals if we adopt the conventional confidence interval, then we may think that it is pos-sible that the number of nucleotide substitutions per site can be near 0.0001. But from the two other intervals, we have 0.95 confi-dence that the number of nucleotide substitutions per site is great-er than 0.002. By the analysis from Sections2 and 3, when the proportion of the different nucleotide between two sequences is low, the score and adjusted confidence intervals are more reliable than the conventional confidence interval.

6. Conclusion

The performance and construction of confidence intervals for evolutionary models have not been much investigated in the liter-ature. In this article, we explore the conventional confidence inter-val performance and provide new confidence interinter-vals, which are shown to be better than the conventional confidence interval for estimating the nucleotide substitution number per site when the true number of substitutions is small. Both the constant substitu-tion rate and the variable substitusubstitu-tion rate cases are considered. One of the proposed approaches is to extend the score confidence interval for the binomial proportion to construct the improved confidence interval. The other approach is based on an adjusted ap-proach used for the binomial distribution. The simulation results show that the proposed confidence intervals are substantially im-proved over the standard confidence intervals. For the constant substitution rate case, the two new intervals outperform the con-ventional interval. For the variable substitution rate case, since the score interval does not have a closed form, the adjusted inter-val is more feasible for the real applications.

Since the proposed methodologies as well as the conventional approach are based on the normal approximation, when the length of sequences is long, the approaches can perform well. But when the length is not long enough, which can be viewed as a case with small sample size case, the methods are not satisfactory. In this cir-cumstance, a more accurate method is to select an appropriate fac-tor value, which can deal with the shorter length case with high precision.

The confidence interval approach can provide more useful information for estimating the nucleotide substitution number than the point estimation. The approach proposed in this article can provide a more efficient and accurate way in the nucleotide

substitution number estimation than the conventional confidence interval.

Acknowledgments

This study was supported by the National Science Council and National Center for Theoretical Sciences in Taiwan.

References

Agresti, A., Coull, B.A., 1998. Approximate is better than ‘‘exact’’ for interval estimation of binomial proportions. Am. Stat. 52, 119–126.

Chor, B., Hendy, M.D., Snir, S., 2006. Maximum likelihood Jukes–Cantor triplets: analytic solutions. Mol. Biol. Evol. 23, 626–632.

Dopazo, J., 1994. Estimating errors and confidence intervals for branch lengths in phylogenetic trees by a bootstrap approach. J. Mol. Evol. 38, 300–304. Dumbacher, J.P., Pratt, T.K., Fleischer, R.C., 2003. Phylogeny of the owlet-nightjars

(Aves: Aegothelidae) based on mitochondrial DNA sequence. Mol. Phylogenet. Evol. 29, 540–549.

Fu, Y., 1995. Linear invariants under Jukes’ and Cantor’s one-parameter model. J. Theor. Biol. 173, 339–352.

Golding, G.B., 1983. Estimates of sequence divergence: An examination of some assumptions. Mol. Biol. Evol. 1, 125–142.

Graur, D., Li, W.H., 1999 Fundamentals of Molecular Evolution. Sinauer Associates, Sunderland, MA.

Holmquist, R., 1971. Theoretical foundations for a quantitative approach to paleogenetics. Part I: DNA. J. Mol. Evol. 1, 115–133.

Jukes, T.H., Cantor, C.R., 1969. Evolution of protein molecules. In: Munro, H.N. (Ed.), Mammalian Protein Metabolism. Academic Press, New York, pp. 21–132. Kaplan, N., Risko, K., 1982. A method for estimating rates of nucleotide substitution

using DNA sequence data. Theor. Popul. Biol. 21, 318–328.

Kimura, M., 1980. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16, 111–120.

Kimura, M., 1981. Estimation of evolutionary distances between homologous nucleotide sequences. Proc. Natl. Acad. Sci. USA 78, 454–458.

Kimura, M., Ohta, T., 1972. On the stochastic model for estimation of mutational distance between homologous proteins. J. Mol. Evol. 2, 87–90.

Kocher, T.D., Wilson, A.C., 1991. Sequence evolution in mitochondrial DNA in humans and chimpanzees: Control region and a protein-coding region. In: Osawa, S., Honjo, T. (Eds.), Evolution of Life: Fossils, Molecules, and Culture. Springer-Verlag, New York, pp. 391–413.

Lanave, C., Preparata, G., Saccone, C., Serio, G., 1984. A new method for calculating evolutionary substitution rates. J. Mol. Evol. 20, 86–93.

Nei, M., Gojobori, T., 1986. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol. Biol. Evol. 3, 418–426.

Nei, M., Kumar, S., 2000. Molecular Evolution and Phylogenetics. Oxford University Press, Inc.

Rosenberg, M.S., 2005. MySSP: non-stationary evolutionary sequence simulation, including indels. Evol. Bioinf. Online 1, 81–83.

Tamura, K., Nei, M., 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol. Biol. Evol. 10, 512–526.

Wakeley, J., 1993. Substitution rate variation among sites in hypervariable region 1 of human mitochondrial DNA. J. Mol. Evol. 37, 613–623.

Wakeley, J., 1994. Substitution-rate variation among sites and the estimation of transition bias. Mol. Biol. Evol. 11, 436–442.

Wang, H., 2007. Exact confidence coefficients of confidence Intervals for a Binomial Proportion. Stat. Sinica 17, 361–368.

Wang, H., 2009. Exact average coverage probabilities and confidence coefficients of confidence intervals for discrete distributions. Stat. Comput. 19, 139–148. Wang, H., 2010. Monotone boundary property and the full coverage property of

confidence intervals for a binomial proportion. J. Stat. Plan. Infer. 140, 495–501. Wang, H., Tzeng, Y.H., Li, W.H., 2008. Improved variance estimators for one- and two-parameter models of nucleotide substitution. J. Theor. Biol. 254, 164–167. Yang, Z., 2007. Computational Molecular Evolution. Oxford University Press.

數據

Fig. 3 shows the coverage probabilities of the confidence inter- inter-vals for L = 100 and 500 when the expected substitution number per site is less than 0.11.
Fig. 2. The expected lengths of the level 0.95 confidence intervals (6) when L = 100 and 500 for K 6 0:06
Fig. 4. Solid and dashed lines present the coverage probabilities for the level 0.95 score and adjusted confidence intervals, respectively when L = 100 and L = 500 for K 6 0:06 for the constant substitution rate case.
Fig. 8. Solid and dashed lines present the expected lengths for the level 0.95 standard and adjusted confidence intervals, respectively when L = 500 and L = 1000 for H 6 0:05 for the variable substitution case.

參考文獻

相關文件

Robinson Crusoe is an Englishman from the 1) t_______ of York in the seventeenth century, the youngest son of a merchant of German origin. This trip is financially successful,

fostering independent application of reading strategies Strategy 7: Provide opportunities for students to track, reflect on, and share their learning progress (destination). •

Strategy 3: Offer descriptive feedback during the learning process (enabling strategy). Where the

How does drama help to develop English language skills.. In Forms 2-6, students develop their self-expression by participating in a wide range of activities

Now, nearly all of the current flows through wire S since it has a much lower resistance than the light bulb. The light bulb does not glow because the current flowing through it

Students are asked to collect information (including materials from books, pamphlet from Environmental Protection Department...etc.) of the possible effects of pollution on our

To this end, we introduce a new discrepancy measure for assessing the dimensionality assumptions applicable to multidimensional (as well as unidimensional) models in the context of

There are existing learning resources that cater for different learning abilities, styles and interests. Teachers can easily create differentiated learning resources/tasks for CLD and