6.4 最大似然法

前面介绍的系统发生分析方法隐含地使用了各种概率模型，说明生物分子序列是如何进化的，通过系统发生树研究序列之间的进化关系。最大似然法（Maximum Likelihood，ML）明确地使用概率模型，其目标是寻找能够以较高概率产生观察数据的系统发生树。最大似然法是一类完全基于统计的系统发生树重建方法的代表。该方法在每组序列比对中考虑了每个核苷酸替换的概率。例如，转换出现的概率大约是颠换的三倍。在一个三条序列的比对中，如果发现其中有一列为一个C，一个T和一个G，我们有理由认为，C和T所在的序列之间的关系很有可能更接近。由于被研究序列的共同祖先序列是未知的，概率的计算变得复杂；又由于可能在一个位点或多个位点发生多次替换，并且不是所有的位点都是相互独立，概率计算的复杂度进一步加大。尽管如此，还是能用客观标准来计算每个位点的概率，计算表示序列关系的每棵可能的树的概率。然后，根据定义，概率总和最大的那棵树最有可能是反映真实情况的系统发生树。

在基于DNA或蛋白质序列的系统发生分析方面，与最大简约法相似，最大似然法首先依赖于一个合理可靠的多重序列的比对，然后检测每一列的变化。对于每一个可能的树，计算在每一列发现真实序列变化的可能性，将每个排列位置的概率相乘，其结果作为每棵树的可能性。具有最大似然值的树就是最可能的树。

对于一棵给定的树，我们希望有一种评价的方法。可以用可能性得分评估我们所作出的假设，即评价所得到的系统发生树T。对于给定的一组分类单元，假设它们的观察值为M（M为向量），可以选择一棵树，使得P(M|T)最大，即最大似然法。下面假定已经知道树T的拓扑结构，目标是计算树T的似然值，寻找T的最优的分支长度。

令t_vu代表节点v和u之间的分支长度，反映的是遗传距离或者进化时间。假设特征是两两独立的，分支过程是一个Markov过程，即节点处于一个给定状态的概率仅仅是其父节点状态及该节点到父节点分支长度的函数。以概率P_x®y(t_vu)表示在时间t_vu内，从状态x转换到状态y的概率。假设在进化过程中特征x的出现频率是恒定的，其值为P(x)。

假设有一个矩阵M，它是关于n个分类单元的实际观察值，M描述每个分类单元m个特征的具体取值。同时假设存在一棵树T，其叶节点（如v、u）对应于这些分类单元，而树中的分支代表分类单元之间的距离t_vu ，求该树的似然值L = P(M|T)。

首先处理最简单的情况，即每个分类单元仅有一个特征。由于树种内部节点的标记未知，需要考虑所有可能的树，并加和对应的结果。例如，对于图6.11所示的树，可以写出公式：

其中，r、v是内部节点可能的标识。

将上述计算公式推广到多个特征，只需分别对每个特征进行反复的计算，然后将结果相乘（因为假定特征两两独立）。

利用最大似然法得到的树一般是有根树，然而，如果特征替换是可逆的，即P_x®y(t)= P_y®x(t)，则树实际上是无根的，在不改变树的似然得分情况下，可以任意选择树根。

对于一个特征j，定义

                    C_j(x,v) = P(根为v的子树|v_j=x)

C_j(x,v) 是关于子树节点v的条件概率，即节点v的第j个特征具有标识x成为子树的可能性。下面利用动态规划算法计算树的可能性。

（1）首先进行初始化工作。对于所有的叶节点v和标识符号x，置

（2）递归计算：按后序方式遍历树（即先计算子树，然后再计算树根），对于每个内部节点v，设其子节点为u和w，对每个标识符x计算

（3）最终结果为

    可以根据给定的树的拓扑结构寻找最优分支长度。这实际上是按照一定方式给各个分支长度赋值，使得整棵树的似然值L 最大。可以采用许多不同的方法使似然值L最大，如Newton-Raphson方法，或者最大期望EM 算法。