Twisst 是一个基于树文件计算一组分类群拓扑权重的软件。它可用于使用群体基因组数据探索分类群关系在整个基因组中的变化。twisst 依赖 ete3包
D:\Python39> pip install ete3
Downloading ete3-3.1.3.tar.gz (4.8 MB)
Successfully built ete3
Installing collected packages: ete3
Successfully installed ete3-3.1.3
--- pip show ete3
Summary: A Python Environment for (phylogenetic) Tree Exploration
Home-page: http://etetoolkit.org
简介
拓扑加权是量化不一定是单系群之间关系的一种方法。它通过考虑更简单的“分类单元拓扑”并量化与每个分类单元拓扑匹配的子树的比例，提供了复杂谱系的摘要。我们用来计算权重的方法称为 Twisst：通过子树的迭代采样进行拓扑权重。
在本次实践中，我们将使用模拟数据来探索拓扑权重如何提供谱系历史。然后，我们将尝试使用针对窄窗口推断的邻居连接树来推断整个模拟染色体的拓扑权重。
工作流程
我们将分析一组谱系，这些谱系代表了在相当复杂的历史（包括种群细分、基因流动和选择）下进化的染色体部分的历史。我们将使用 twist 计算该基因组区域的拓扑权重，然后在 R 中探索结果。
1. 模拟家谱分析
下载代码和数据
这部分实践的脚本和示例数据位于 github 上的 twisst 包中。
wget https://github.com/simonhmartin/twisst/archive/v0.2.tar.gz
cd D:\生物信息学
tar xzf twisst-0.2.tar.gz
我们将使用的示例数据由编码为 Newick 树的家谱文本文件组成。在本例中，树木是使用模拟器 msms 进行模拟的。如果我们有真实数据，我们将不知道这些树，并且必须使用 Relate、tsinfer 等工具来推断它们，或者仅在狭窄的窗口上运行系统发育推断。
我们可以查看文件中的第一棵树：
zcat twisst-0.2/examples/msms_4of10_l50k_r500_sweep.trees.gz | head -1
每个 : 之前的数字是样本名称。 : 之后的数字是分支长度。在本教程中，我们将仅考虑树的形状，而不考虑分支长度。我们还可以检查染色体该区域的不同谱系的总数：
zcat twisst-0.2/examples/msms_4of10_l50k_r500_sweep.trees.gz | wc -l
为了绘图，我们还需要知道这些谱系出现在染色体上的位置。该数据在第二个文件中提供，其中包含三列：每个谱系的染色体、开始和结束。该文件与树文件具有相同的行数。
zcat twisst-0.2/examples/msms_4of10_l50k_r500_sweep.data.tsv.gz | head
zcat twisst-0.2/examples/msms_4of10_l50k_r500_sweep.data.tsv.gz | wc -l
正如您所看到的，该模拟数据集中的一些谱系占据了非常狭窄的染色体区域，小至 1 bp。经过许多代的重组，情况可能是，对于给定的一组样本，谱系在整个染色体上有细微的变化。在这种情况下，我们知道每个地区的真实谱系 - 尚未推断出来。
计算拓扑权重
我们运行 twist 来计算每个拓扑的权重。Twist 需要的唯一信息是：
树文件的名称（用 -t 指定）
输出权重文件的名称 (-w)
每个组的名称以及属于该组的样本 (-g)。
分组可以根据物种、表型或地理（或任何你喜欢的）来确定。在我们的例子中，有四组，每组 10 个单倍体样本。 A 组由 1:10 的样本组成，B 组由 11:20 的样本组成，依此类推。
python twisst-0.2/twisst.py \
-t twisst-0.2/examples/msms_4of10_l50k_r500_sweep.trees.gz \
-w msms_4of10_l50k_r500_sweep.weights.tsv.gz \
-g A 1,2,3,4,5,6,7,8,9,10 \
-g B 11,12,13,14,15,16,17,18,19,20 \
-g C 21,22,23,24,25,26,27,28,29,30 \
-g D 31,32,33,34,35,36,37,38,39,40 \
--outgroup D
D:\生物信息学> python twisst-0.2/twisst.py -t twisst-0.2/examples/msms_4of10_l50k_r500_sweep.trees.gz -w msms_4of10_l50k_r500_sweep.weights.tsv.gz -g A 1,2,3,4,5,6,7,8,9,10 -g B 11,12,13,14,15,16,17,18,19,20 -g C 21,22,23,24,25,26,27,28,29,30 -g D 31,32,33,34,35,36,37,38,39,40 --outgroup D /-A /-A /-A/-| /-| /-|| \-B | \-C | \-D
--| --| --||--C |--B |--B| | |\-D \-D \-Ctwist 将考虑所有可能的样本组合，其中每组只有一个样本。例如，检查的第一个组合将是样本 1、11、21 和 31，分别代表 A、B、C 和 D 组。忽略树中的所有其他分支，twist 记录仅包含四个感兴趣样本的“子树”的拓扑，该子树可能具有三种可能形状之一：(((A,B),C),D), (( (A,C),B),D) 或 (((B,C),A),D) （请注意，这里的树表示为有根树，其中 D 为外群。我们可以告诉 twistt 哪个是外群（--outgroup），因此它将树显示为正确的根，但这不会影响结果）。
检查结果是什么样的：
# first 10 lines
zcat msms_4of10_l50k_r500_sweep.weights.tsv.gz | head -10
# total number of lines
zcat msms_4of10_l50k_r500_sweep.weights.tsv.gz | wc -l
权重文件中的三列代表三种拓扑，这三种拓扑也在文件中定义。这些数字不是比例，而是代表每个拓扑的组合总数。每行的总和应为 10^4 = 10,000，因为每行有四组样本。
您将看到一些相邻线具有相同的权重。这与一些重组事件改变样本之间的关系，但不影响权重的方式有关。
参考：生信教程：使用拓扑加权探索基因组进化（1）
)
