生信教程：使用拓扑加权探索基因组进化（1）

数据科学工厂

动动发财的小手，点个赞吧！

使用 Twisst 探索整个基因组的进化关系的拓扑加权教程。

简介

拓扑加权是量化不一定是单系群之间关系的一种方法。它通过考虑更简单的“分类单元拓扑”并量化与每个分类单元拓扑匹配的子树的比例，提供了复杂谱系的摘要。我们用来计算权重的方法称为 Twisst：通过子树迭代采样进行拓扑权重。

在本次实践中，我们将使用模拟数据来探索拓扑权重如何提供谱系历史。然后，我们将尝试使用针对窄窗口推断的邻居连接树来推断整个模拟染色体的拓扑权重。

工作流程

我们将分析一组谱系，这些谱系代表了在相当复杂的历史（包括种群细分、基因流动和选择）下进化的染色体部分的历史。我们将使用 twist 计算该基因组区域的拓扑权重，然后在 R 中探索结果。

1. 模拟家谱分析

下载代码和数据

这部分实践的脚本和示例数据位于 github 上的 twisst 包中。

#download the twisst package zip file from github
wget https://github.com/simonhmartin/twisst/archive/v0.2.tar.gz

#extract the files from the zipped file
tar -xzf v0.2.tar.gz

#remove the zipped file
rm v0.2.tar.gz

我们将使用的示例数据由编码为 Newick 树的家谱文本文件组成。在本例中，树木是使用模拟器 msms 进行模拟的。如果我们有真实数据，我们将不知道这些树，并且必须使用 Relate、tsinfer 等工具来推断它们，或者仅在狭窄的窗口上运行系统发育推断。

我们可以查看文件中的第一棵树：

zcat twisst-0.2/examples/msms_4of10_l50k_r500_sweep.trees.gz | head -n 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

twist 将考虑所有可能的样本组合，其中每组只有一个样本。例如，检查的第一个组合将是样本 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 -n 30
#total number of lines
zcat msms_4of10_l50k_r500_sweep.weights.tsv.gz | wc -l

权重文件中的三列代表三种拓扑，这三种拓扑也在文件中定义。这些数字不是比例，而是代表每个拓扑的组合总数。每行的总和应为 104 = 10,000，因为每行有四组样本。

您将看到一些相邻线具有相同的权重。这与一些重组事件改变样本之间的关系，但不影响权重的方式有关。

分析结果

打开 R 或 RStudio，如有必要，将工作目录设置为保存文件的位置。您可以使用 setwd() 命令，或者在 RStudio 中使用菜单。启动一个新的R脚本来记录命令首先，我们将导入一组随 twistt 一起分发的函数，这将有助于绘图。

source("twisst-0.2/plot_twisst.R")

请注意上面脚本的酷名称。我们定义包含每个谱系的权重以及染色体上每个块的开始和结束位置的文件。

#weights file with a column for each topology
weights_file <- 'msms_4of10_l50k_r500_sweep.weights.tsv.gz'

#coordinates file for each window
window_data_file <- 'twisst-0.2/examples/msms_4of10_l50k_r500_sweep.data.tsv.gz'

我们已经知道这两个文件的结构，但我们将使用方便的 import.twisst 函数，而不是直接读入它们并使用它们。

twisst_data <- import.twisst(weights_file, window_data_file)

使用提供的plot.twisst 函数绘制原始权重。

pdf("simulted_trees_weights.pdf", width=8, height=6)
plot.twisst(twisst_data)
dev.off()

这会将 pdf 文件写入您正在工作的目录 - 打开它！

图顶部的树显示了我们加权的 3 种不同拓扑。下图显示了权重。您将看到不同宽度的颜色列。每列对应一个具有独特谱系的块。有些块都是一种颜色，并且值达到 1。这表明该块中的所有子树具有相同的拓扑，表明谱系一致且完全排序。其他柱子有两种或多种颜色叠加，表明家谱具有更复杂的进化历史，个体在群体之间跳跃。完全随机的谱系，其中不存在按组聚类的情况，对于所有三种拓扑具有相同的权重。

通常需要平滑权重，以便我们可以更清楚地看到它们在染色体上的变化。使用 smooth.twisst 函数创建平滑权重并重新绘图。

twisst_data_smooth <- smooth.twisst(twisst_data, span_bp=5000)

pdf("simulted_trees_weights_smooth.pdf", width=8, height=6)
plot.twisst(twisst_data_smooth)
dev.off()

这对 5 kb 窗口的权重进行了平均。

现在我们更清楚地看到，占主导地位的拓扑是 topo1 和 topo3。在本例中，模拟涉及根据 topo1 进行种群分裂，但模拟了从 C 到 B 的适应性渗入，这就是为什么 topo3 比 topo2 更普遍，也是为什么 topo3 在该区域中部有一个大尖峰的原因。这是所选轨迹的位置。我们也可以查看权重的总体分布，或者只检查平均值。

plot.twisst.summary.boxplot(twisst_data)

twisst_data$weights_mean

正如预期的那样，模拟历史遵循 topo1，这是最丰富的拓扑。基因渗入产生了过量的 topo3。

如果类群最近分裂，谱系排序通常是不完整的，即使谱系完整，我们也可以发现由于过去谱系排序的随机性而与“物种树”不一致的谱系。如果仔细观察我们制作的第一个图，您可能会发现一个狭窄的窗口，其中 topo2 的权重为 1。这表明谱系完全排序，但不一致。

未完待续!

微信扫一扫分享文章