6.2 基于距离的系统发生树构建方法
建立系统发生树的基本任务是：在给定的条件下（包括分类单元、分类单元的特征值或者序列），构造一棵最优的系统发生树。这里重点讨论针对DNA序列或者蛋白质序列构建系统发生树。
基于距离的系统发生树构建方法的基本思路是：给定一种序列之间距离的度量，在该距离度量下构建一棵系统发生树，使得该树能够最好地反映已知序列之间的距离。这种方法采用两两距离，建立一个距离矩阵，如表6.2所示，根据距离矩阵构造系统发生树。
表6.2 10条核酸序列的距离矩阵。
|
|
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
|
2 |
0.0516 |
|
|
|
|
|
|
|
|
|
3 |
0.0550 |
0.0031 |
|
|
|
|
|
|
|
|
4 |
0.0483 |
0.0221 |
0.0253 |
|
|
|
|
|
|
|
5 |
0.0582 |
0.0651 |
0.0685 |
0.0549 |
|
|
|
|
|
|
6 |
0.0094 |
0.0416 |
0.0450 |
0.0384 |
0.0549 |
|
|
|
|
|
7 |
0.0125 |
0.0584 |
0.0619 |
0.0551 |
0.0651 |
0.0157 |
|
|
|
|
8 |
0.0284 |
0.0687 |
0.0722 |
0.0654 |
0.0754 |
0.0317 |
0.0285 |
|
|
|
9 |
0.0925 |
0.1221 |
0.1259 |
0.1185 |
0.1370 |
0.0820 |
0.0786 |
0.0927 |
|
|
10 |
0.1921 |
0.2183 |
0.2228 |
0.2054 |
0.2309 |
0.1798 |
0.1795 |
0.1833 |
0.1860 |
|
|
这里的距离代表两条序列之间各位点核苷酸替换数的估计值。 |
6.2.1 最小二乘法
在很多情况下，这种距离矩阵传达了大部分进化信息。运用前面介绍的序列比较程序，可以计算出序列之间的距离。但是，在进行序列比较时，应根据比较的对象是DNA序列还是蛋白质序列，正确选用不同的打分矩阵。
为了便于分析，首先定义一种连续加和距离函数，在该函数下，两个分类单元之间的距离与系统发生树中连接这两个分类单元的分支总长度成正比。这样，如果分类单元a和分类单元b由经过中间节点v的两条边相连，两条边的长度分别为dav和dbv，则它们之间的距离为dav+ dbv。这样，可以在系统发生树中确定a和b的相对位置。进一步，假设三个分类单元之间的距离分别为dab、dac、dbc,如果分类单元a和分类单元b由经过中间节点v的两条边相连，再经过节点u与分类单元c相连,则可以通过求解线性方程计算出系统发生树的各种内部距离。
例，如果有三个分类单元，其两两距离如下：
dab = 0.5;
dac = 0.9;
dbc = 0.9
假设分类单元a和分类单元b的分歧起始时间是相同的，根据分子时钟假说， 
和 
的值应该是相等的，进一步假设节点u到其它节点的距离相同，则通过求解方程，得到如图6.2所示的一棵树。
图6.2是一个简单的例子，但是，在实际应用中，所要处理的分类单元可能很多，因而，需要求解的线性方程也很多，难以求解，或者方程组的求解过程存在着不确定性。因此，需要采用数学逼近的方法。
统计学上逼近的方法之一是最小二乘法。我们的目标是构造一棵树T，以该树的叶节点代表分类单元，用该树预测分类单元之间的距离。通过优化，使下式最小化：
这里，Dij为分类单元i和j的实际观察距离（或序列之间的计算距离），dij是分类单元i和j在系统发生树T 中的距离，Wij是与分类单元i和j相关的权值。SSQ（T）是树T所有预测值与实际观察值偏差的累加和。权值Wij一般为1，或
寻找一棵最小方差树是一个NP-完全问题，需要采用近似的算法。下面讨论几种计算时间复杂度为多项式的启发式方法，即连锁聚类方法（linkage clustering）、非加权组平均法（UPGMA）和邻近归并法（Neighbor Joining）。
6.2.2 连锁聚类方法及非加权分组平均法
连锁聚类属于一般的聚类分析方法，当用来构建系统发生树时，其假定的前提条件是：在进化过程中，核苷酸或氨基酸的替换速率是均等且恒定的，在每一次分歧发生后，从共同祖节点到两个分类单元间的分支长度一样。在构建系统发生树时，首先用n个叶节点表示n个分类单元（序列）， 每个分类单元自成一类，然后通过反复的聚类使所有的分类单元都聚为一类，并将进化过程中的祖先赋予树的内部节点，最终得到一个完整的系统发生树。假设若干条序列是从一个共同的祖先进化而来，则系统发生树将是一个有根树，并且从根节点出发到所有叶节点路径的长度相同。
对于给定的序列，通过序列之间的两两比对，计算序列之间的进化距离，然后根据距离矩阵构造系统发生树。算法的基本思路是首先从距离矩阵中选择距离最小的一对分类单元（序列），令它们分别为x和y，然后将这两个分类单元合二为一，形成一个新的对象（代表这两个分类单元的祖先，记为z），并重新计算这个新的对象与其它分类单元(或对象，以u表示)之间的距离d(z,u)。不同的实现方案采用不同的计算公式：
单连锁聚类: d(z,u)=min(d(x,u),d(y,u)) （6-8）
最大连锁聚类： d(z,u)=max(d(x,u),d(y,u)) （6-9）
平均连锁聚类： d(z,u)=(d(x,u)+d(y,u))/2 （6-10）
上述公式中z代表x和y的合并，u代表任意其它对象。每次合并所形成的新对象实际上是一个聚类，以一个内部节点表示，该节点到x、y 所在节点的距离相同，其值等于d（x，y）的一半，而到其它节点的距离按照上述公式计算。每次合并后，修改距离矩阵。重复上述过程，直到所有的分类单元都被合并到一类为止。
连锁聚类算法的执行过程如下：
（1） 初始化：使每个分类单元自成一类，如果有n个分类单元，则开始时共有n个类，每个类的大小为1，分别用n个叶节点代表每个类；
（2） 执行下列循环：
l 寻找具有最小距离dxy 的两个类x、y；
l 建立一个新的聚类z, 该聚类是类x和类y的合并；
l 在树中建立一个新的内部节点z，生长两个新的分支，将x 和y 连接到z，并使这两个分支到x、y 中各个叶节点的长度为d（x，y）/2；
l 按照公式（6-8）、（6-9）或（6-10）计算新的分类到其它类的距离；
l 在距离矩阵中删除与类x和类y 相应的行和列，为类z 加入新的行和列，其值按上述公式计算；
重复循环，直到仅剩一个类为止。
系统发生树的结构与聚类过程相对应。图6.3就是一个通过平均连锁聚类分析的实例。距离矩阵如图6.3（a）所示。初始化以后，共有5个分类，分别对应于5个分类单元A、B、C、D及E。第一次聚类时，选择距离最小的两个分类A和 C进行合并，形成一个新的分类（AC）。然后修改距离矩阵，（AC）到B的距离等于（8+8）/2=8，其它距离类推。最终结果见图6.3 (b)。在系统发生树中，两个分类单元之间的距离为连接这两个分类单元所在叶节点对应路径的长度。
运用连锁聚类算法，得到一棵有根的系统发生树，从树根到任何叶节点的分支长度全都一样，也就是说，所有物种的突变速率相同，存在一个固定节律的“分子钟”，各个物种从树根（代表某个历史起点）出发，踏着同样的节律，沿不同的路径，演化成为当前的形式。而任意两个物种之间的距离是连接这两个物种所在叶节点路径上各分支长度之和。
如果上述过程处理的分类单元是序列的话，则聚类的过程实际上也是一个进行序列多重比对的过程，将序列多重比对问题转化为序列两两比对问题。
非加权分组平均法（Unweighted Pair Group Method with Arithmetic mean, UPGMA）是一种较为常用的聚类分析方法，最早是用来解决分类问题的。使用该方法的前提条件也是稳定的进化速率。从分类过程来看，该方法与前面介绍的平均连锁聚类方法相似，但是关于各个分类（即若干个分类单元形成的集合）之间距离的计算公式不一样。在平均连锁聚类过程中，一个新类到其它类之间的距离就是简单的原距离平均值。这样的计算非常简单，但是，如果各个类中分类单元个数不一样，原距离矩阵中各个距离值对新距离计算的贡献就不一样，或者说是经过“加权”的，称这样的聚类为加权分组平均（Weighted pair group method with Arithmetic mean, WPGMA）。在非加权分组平均法中，在计算新分类到其它分类之间的平均距离时按照各分类中分类单元的数目进行加权处理。令类x和类y所形成的新分类为（xy），新分类到其它类u的距离按照下述公式计算：
其中nx、ny、(nx+ny)分别为x类、y类、(xy)类的元素个数。这样，如果每个距离对最终结果的贡献一样，即是“非加权”的。
图6.4是针对表6.2中的距离矩阵利用UPGMA 而构建系统发生树。
6.2.3 距离变换法
连锁聚类和UPGMA算法的一个缺陷是假定所有家系的进化速率是相同的，但是，实际情况并不总是这样。进化速率的变化容易导致连锁聚类和UPGMA算法产生错误拓扑结构的树。假设有4个分类单元A、B、C和D ，其系统发生关系及各个分类单元之间的距离如图6.5所示，距离矩阵见表6.3。如果利用UPGMA进行分析，则首先合并A和C
基因组信息分析