决策该树收缩的截止点算法？
我有一棵Newick树，这棵树是通过比较一些DNA调控序列的相似性（欧几里得距离）来构建的。这些调控序列的长度在4到9个碱基对之间。
这棵树的互动版本可以在iTol上找到（这里），你可以随意操作——设置好参数后，只需点击“更新树”按钮：
我的具体目标：如果某些调控序列到最近的父类的平均距离小于X，就把这些序列合并在一起（这些序列在树的末端）。这在生物学上很有趣，因为一些基因调控的DNA序列可能是同源的（即它们之间有亲缘关系）。这种合并可以通过上面提到的iTol图形界面来完成，比如如果你选择X = 0.001，那么一些调控序列就会合并成三角形（调控序列家族）。
我的问题：有没有人能建议一个算法，能够输出或帮助可视化哪个X值最适合“最大化合并调控序列的生物学或统计相关性”？理想情况下，当把某个属性与X值绘制在一起时，树的某些特性会出现明显的变化，这样算法就能找到一个合理的X值。有没有已知的算法、脚本或软件包可以做到这一点？也许代码会将某个统计数据与X值进行绘图？我尝试绘制X与平均聚类大小的关系（matplotlib），但没有看到明显的“跃升”来告诉我该使用哪个X值：
我的代码和数据：我的Python脚本链接在[这里][8]，我在代码中做了很多注释，它会为你生成树数据并绘制上面的图（使用参数d_from、d_to和d_step来探索距离的截止值X）。如果你有easy-install和Python，只需执行这两个bash命令即可安装ete2：
apt-get install python-setuptools python-numpy python-qt4 python-scipy python-mysqldb python-lxml
easy_install -U ete2
2 个回答
我觉得在给出具体建议之前，我需要了解更多信息。不过，也许这些内容能帮到你。我假设每个终端节点代表一个序列，而每个内部节点代表一个PSSM（位置特异性得分矩阵）。
X的计算是根据具体应用而定的。举个例子，如果你想合并超平行基因（ultraparalogs），得到的X和你想合并所有同源基因（homologs）时得到的X是不一样的。
由于基因是通过复制和物种形成不断产生的，所以没有一个固定的X值能根据进化关系来区分序列。因此，我不认为仅仅通过查看聚类统计数据就能找到一个令人满意的方法来判断序列之间的进化关系。
一个更严谨的方法是从每个调控基序的基因构建一个基因树，并将其与物种树进行对比。市面上有相关的软件和一些额外的技巧可以帮助识别直系同源基因和旁系同源基因。
如果你这样做，你的树的内部节点将会标注上推测出的进化事件（例如，基因复制、物种形成）。然后你可以沿着树向上走，合并那些你不关心的节点。
你可以尝试使用类似于树的对比方法，正如@Jeff提到的。不过，标准的树对比方法实际上可能会失败。
树的对比首先是要在目标树上添加一些代表“损失”的分支，这些损失是指在进化过程中某些特征的消失。然后要标明在哪些节点上发生了特征的“重复”。损失和重复的加权总和提供了一个成本函数，用来进行优化。
但在你的情况下，你想解决的问题是“把这个超级树拆分成合适大小的同源子树”。这意味着你不太想给损失打分，而是更关注重复。你需要一种打分方法，能够显示有多少个同源子树合并到了你的超级树中。因此，你可以尝试以下的打分方法：
- 从超级树中，计算重复物种的数量，记为S1。
- 将所有的终端叶子（即末端节点）合并，这些叶子是平行同源物种，然后计算新的重复物种数量，记为S2。
- S1和S2之间的差值大致显示了你的超级树中有多少个子树。
- 为了修正由于超级树大小不同而造成的偏差，用超级树中独特物种的数量N来进行除法。
如果我们把这个分数称为“子树因子”，那么它的计算公式就是：
S1 - S2 / N
推论：
如果 S1 - S2 = S1，这意味着你的超级树大约只有一个真实的子树，所有重复物种的出现只是由于最近的平行同源现象。
如果 S1 - S2 = 0，这意味着你的超级树大约有 S1 个真实的子树。