以下内容都是基于Christopher P. Randle教授在研究组上交流时的课件整理而来。
Stationary models for Maximum Likelihood estimation from DNA data
基于DNA数据，用最大似然法进行系统发育重建的模型
在系统发育学中，模型就是为了确定给定分支长度的分枝上发生某种性状状态改变的概率。分子性状状态的改变就是碱基的替代（substitution）。
为了建立模型，我们必须给出一些前提条件。这些前提条件可能会和现实情况不符，目的只是为了让事情简单化，使得我们有可能建模。这些前提条件包括：
1碱基替代是一个马尔科夫过程（Markov Process）
Markov Process: Changes on one branch occur independently of changes on another branch. This means that the probability of change depends only on the state optimizations.
是指每个分支上发生的情况都是相互独立的，仅受到分支末端的taxa的形状状态（character state）的影响。
2 碱基的替代是时间可逆的。
Time Reversibility: The direction of change is immaterial for the calculation.
因为似然值的估计常常是在无根树上进行的，所以 性状状态的改变（这里就是碱基的转变或碱基替代）也是没有方向的，比如，由A到C的改变和由C到A的改变发生的概率是一样的：
p(A->C）=p(C>A）
3 均一过程（Homogeneity）:
The process is homogenous across characters and throughout the branching process.
整个演化的过程在不同的性状和分支过程中都是同质的、均一的。
我们从前面知道，建立模型最重要的就是找到参数，来表示给出的假设。
Models of DNA sequence evolution include parameters for estimating the probability of transformations between nucleotides along branches with associated branch length.
建立DNA序列的演化模型需要的参数就是用来估算指定分支长度的分枝上发生某种碱基替代的概率。
那么我们怎么表示这个参数那？
Parameters take the form of elements in an instantaneous rate matrix, or Q-matrix, where each parameter specifies the probability of change over an infinitesimal increment of time. Ultimately, the Q-matrix will be used to estimate a matrix or transformation probabilities for a branch of a defined length.
这个参数是以一个矩阵中的元素形式表现出来的，被称为Q-矩阵。这个矩阵中的每一个元素就是一个参数，表示发生这种碱基替代的概率。
下面我们来介绍最大似然法中最重要的一个模型：General Time Reversible Model,，简称GTR。
它的元素就是将上图矩阵中每个元素对应的概率的的自然对数：
每一个元素都可以拆分为这种表达式：
=μ×
×
μ代表的是平均瞬时替代速率参数（mean instantaneous substitution rate），对矩阵的所有元素，该值恒定不变。
代表的是相对瞬时速率（relative rate parameter），在矩阵的不同元素中是变化的，用希腊小写字母表示。
代表的是j性状状态在数据矩阵中的出现频率的期望值，被称为期望参数（frequency parameter）。
GTR
我们说过μ在所有元素中都是数值恒定的，我们可以不去管。
根据碱基的替代是时间可逆的前提条件（比如，碱基A、C的相互转变） ，A到C和C到A两种碱基替代的相对瞬时速率（relative rate parameter）是相同的，因此我们现在有6种相对瞬时速率。
在一瞬间，碱基发生转变或者不发生转变，转变发生或不发生的概率是1， 即p（A->A）+p（A->C）+p（A->G）+p（A->T）=1。因此当碱基没有改变的情况下，如A->A,p（A->A）的值为1减去改行其他元素的值。
需要注意的是，我们在前面说过，因为似然值通常很小，是小数点后好几位的数，为了避免麻烦，我们通常用自然对数（lnN）来表示，在GTR中，每个元素也用概率值的自然对数来表示，因此每一行的加和就是ln(1)，为0，p（A->A）位置上的就是0减去其他三个元素的加和，也就是其他三个元素加和的负值。
By restricting the kinds of changes than can happen, we can make specialized models from GTR.
其他的模型就是在GTR模型的基础上，增加一些前提条件，一些参数假定是不变的，使得参数表达变得更简单一些。
F846（也称为HKY85，Hashagano Kishino Yano 1985）模型不再考虑每个碱基对之间转变的相对瞬时速率，而是把碱基转变分成两大类：颠换（transversion，嘌呤和嘧啶之间相互转变）和转换（transition，嘌呤或嘧啶之间的转变），并且假设颠换和转换各有一个相对瞬时速率a、b。
HKY85
更进一步，我们可以假设期望参数在元素中也是恒定的，π=0.25，这就是Kimura 2-Parameter（K2P）模型：
K2P
最简单的模型是假设相对瞬时速率和期望参数都恒定，这就是Jukes-Cantor（JC69）模型：
JC69
我们可以总结一下所有这些模型都是homogeneous和nested。
homogeneous（同质的）指的还是演化过程在不同的性状和树的分支上都是同质的、均一的。
nested（嵌套的）指的是所有的模型都是在GTR这一模型的基础上增加前提条件，变化来的。
Estimating likelihood given a model
用已有模型估算似然值
我们前门提到的Q-矩阵表示的只是瞬时替代速率，而我们要计算的碱基替代概率，因此，要把Q-矩阵过积分的数学方法，转化为代表碱基替代概率的P-matrix。它和Q之间可以用下面的算式推导：
要把我们之前提到的模型由给出的Q-矩阵转化为P-矩阵需要一些积分的数学方法，在这里只给出转化完之后的例子。
我们以K2P模型（转换和颠换的速率不同，期望参数相同）为例，转换和颠换的速率比率（α/β）为κ，则P-矩阵中的元素不外乎下面三种情况，其中i、j表示的是两种性状状态：
第一种表示当没有发生碱基替代时；
第二种表示发生的碱基替代是转换；
第三种表示发生的碱基替代是颠换。
我们假设μ= 0.3，κ=4.0，计算t=0.1的时长内的概率：
这时的P-矩阵为：
目前为止，我们还没有计算过一个树的似然值，那么下面我们用一个最简单的例子来试一下：
对这个给定的树，让我们估算一些它的似然值：
首先把各个分支的碱基替代概率标注在分枝上：
然后我们计算整棵树的似然值，就是各个分支的概率的乘积，而它们的最近祖先是A的概率为0.25，因为共有四种碱基（A、C、T、G）：
这棵树的似然值就是给定这棵树的分支长度和拓扑结构后，出现这样一棵树的概率，用概率的自然对数表示。
Rate heterogeneity 速率的不均一性/异质性
根据分子钟演化假说（clock-like evolution），在整个系统发育树中，所有位点的碱基替代速率应该是均一的，那么分支长度（代表发生了多少演化，这里也可以理解为发生了多少碱基替代）就只和时间有关。但是，实际情况往往不是这样的，碱基的替代速率在不同的分支上是不同的。就像龟兔赛跑，跑过同样长度的距离，一个是因为用的时间比较长，一个是因为跑的快，放在我们这里就是：同样的分支长度，可能是因为演化的时间长，也可能是碱基替代速率快。
这里我们给出一个名词：rate heterogeneity（速率的不均一性）
是什么导致了性状演化速率的不均一性？
原因可能是一些位点对生物的生存至关重要，比如指导一些重要大分子物质的合成，因此承受更大的选择压力，碱基的替代不太可能发生。或者在同一个密码子中，第一、二个位点的变化会造成编码氨基酸的改变（非同义突变），第三个位点的变化容易造成编码氨基酸的改变（同义突变），理论上第三个位点的演化速率就会比前两个位点快。
注：最新研究又认为同义突变也会影响生物功能。这里只是举个例子帮助理解，毕竟科学总是在进步的。
现在我们要把性状演化的不均一性考虑打模型中该怎么办？两种办法：
一种方法是假设一些位点是不变的，另一些可以自由转变。
对每个位点，我们都计算两个似然值：一个是整个树都没有发生变化的概率，如果这个位点的碱基在所有分类群中都不变，那么
的值为1，如果在分类群中不一样，
的值为0；一个是如果发生了碱基替代时每个位点的似然值
，该值的计算就像我们刚才举例的一样。
此外，还要增加一个参数𝝷，代表的是系统发育树的所有位点中未发生改变的位点的比例。
对一个性状j来说，树的似然值的计算式就是：
=𝝷
+(1-𝝷)
第二种方法是对每个位点都增加一个额外的速率参数：𝚪。
𝚪 apportions rate heterogeneity among sites according to a gamma distribution.
𝚪参数基于一个伽马分布来分派各个位点的速率不均一性。
这里我们就要介绍一下伽马分布：
Gammadistribution: A representation of moderately-very asymmetrical probability densities.
我们从图中可以看到伽马分布大概就是一条钟形曲线，居中值的y轴数值高，极端值的y轴数值低，作为概率密度曲线能比较好的反应出我们遇到的情况：大多数位点的演化速率是比较接近的，但存在个别位点的演化速率很快或者很慢。
Gamma Correction: Usually, once a shape and scale parameter are determined, the gamma distribution is broken up into rate classes(rj) and each characterassigned a rate parameter value proportional to the probability density ofgamma.
确定了伽马曲线的形状和数值范围之后，我们会把伽马分布进行分割，根据所代表的伽马分布的概率密度划分出等级（），把这些等级分配给对应的位点，代表该位点的异质性的演化速率。
P的表达式就变为：
有的时候，同一个模型中会同时使用这两种方法，即同时使用参数𝝷和𝚪。但我们要明白它们都是用来表示一件事：among-site rate heterogeneity。
𝝷和α的值是相关的，有时候会造成计算的不准确性。如果不发生变化的位点数目不多的话，我们只需要用𝚪参数就能模拟出位点之间的演化速率差异。
Maximum Likelihood 最大似然值
到现在为止，我们只说了如何计算一个给定树的似然值，那么如何找到我们想要的最优树呢？
最优树就是似然值最大的树。在这一点上，最大似然法和最大简约法很像。但是最大简约法更直观一些，而最大似然法要同时考虑多个方面同时达到最优：拓扑结构、分支长度和模型参数值等。但是模型的参数值对于解释演化过程并不会提供太多有用的信息，但分析时又不能避免，我们恶毒地称这些参数为nuisance parameter。
因此在计算上，最大似然法要比最大简约法复杂很多，也慢很多。但是计算机技术的不断进步已经在改善这个问题。
在选取最优树上，最大似然法和最大简约法很像：
选择一个起始树，不断调整分支长度和模型参数直到似然值达到最大，然后在该树的基础上变换出另一种拓扑结构，再次调整分支长度分布和模型参数，直到似然值达到最大，就这样一直进行下去。
事实上，为了提高搜寻最优树的速度，Q-矩阵的参数通常是保持不变的，调整的只是拓扑结构和分支长度，参数的起始值基于一个最大简约树估算而来。
网友评论