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

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

简介

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

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

工作流程

我们将通过使用 msprime 运行合并模拟，然后直接从输出计算拓扑权重，来探索人口统计参数如何影响权重，所有这些都在 Python 中进行。

使用树序列格式的拓扑加权

到目前为止，我们已经分析了树文件，这些文件对于具有独特谱系的每个染色体“块”具有独特的树，或者在推断树的情况下对于每个窗口具有独特的树。这种格式有点浪费，因为相邻的树通常非常相似，通常仅因单个重组事件而不同，该重组事件将一个分支从树中的一个点移动到另一个点。

树序列格式非常高效，因为它不仅记录树上节点之间的连接，还记录每个连接所在的染色体的长度。 msprime 和 tskit 包使用此格式，它还提供了更多信息。

模拟树序列

我们将使用 msprime 来模拟树序列。 msprime 是一个合并模拟器，这意味着它的工作原理是计算任意两个个体在过去某个特定时间拥有共同祖先的概率。在单一种群中，这是由种群规模决定的。对于多个种群，这还受到种群之间的迁移率以及它们从单一祖先种群传承多久的影响。

我们将使用 Python 交互式环境来处理 msprime 和树序列，并使用 twisst 中的函数来分析它。

为了确保我们可以从 python 中导入 twisst 模块，请将其添加到我们的 python 路径中。

export PYTHONPATH=$PYTHONPATH:twisst-0.2

现在打开一个Python交互式会话（输入“python”），并在文本编辑器中打开一个脚本，因为我们将多次修改并重新运行其中一些行以测试不同模拟参数的效果 导入所需的模块。

import msprime
import matplotlib.pyplot as plt
import twisst

我们将从对单个总体的 10 个样本进行简单模拟开始。我们必须指定长度和重组率，因此 msprime 将为我们提供一个由多个谱系组成的序列，通过重组分开。在这里，我们还指定了随机种子，只是为了确保在这种情况下我们都得到相同的模拟。

ts = msprime.simulate(sample_size=10,
                      Ne=1000,
                      length=10000,
                      recombination_rate=5e-8,
                      random_seed = 1)

我们可以检查树序列中有多少个不同的家谱。

ts.num_trees

并使用树序列对象提供的良好可视化方法来查看它们。

for tree in ts.trees():
    print("interval = ", tree.interval)
    print(tree.draw(format="unicode"))

包含四个种群和基因流的模拟

我们现在将建立一个包含多个群体的更大的模拟。这需要一些不同的组件，我们将分别定义它们。首先，我们定义样本数量和每个总体的总体规模。

pop_n = 10
pop_Ne = 10000

population_configurations = [msprime.PopulationConfiguration(sample_size=pop_n, initial_size=pop_Ne),
                             msprime.PopulationConfiguration(sample_size=pop_n, initial_size=pop_Ne),
                             msprime.PopulationConfiguration(sample_size=pop_n, initial_size=pop_Ne),
                             msprime.PopulationConfiguration(sample_size=pop_n, initial_size=pop_Ne)]

接下来我们设置种群之间的迁移率。这些可以以 4x4 矩阵的形式定义，每个条目给出从另一个群体（行）迁移到另一个群体（列）的速率。这里我们设定第二和第三种群之间双向适度迁移，其他种群之间不迁移。该值代表m，即每代移民占人口的比例。

migration_matrix = [[0,    0,    0,    0],
                    [0,    0,    1e-4, 0],
                    [0,    1e-4, 0,    0],
                    [0,    0,    0,    0]]

最后，我们设置了分割时间，在合并世界观中，这些时间是时间倒流的连接。在 msprime 中，这些称为大规模迁移。因此，前两个种群之间的分裂被建模为所有个体从第二个种群到第一个种群的大规模迁移（称为 0 和 1，因为 Python 从 0 开始编号）。再往前追溯，种群 2 和 3 也大规模迁移到种群 0。在第一个事件（时间向后）之后，我们还关闭了 1 和 2 之间的迁移。

t_01 = 1000
t_02 = 5000
t_03 = 10000

demographic_events = [msprime.MassMigration(time=t_01, source=1, destination=0, proportion=1.0), # first merge
                      msprime.MigrationRateChange(time=t_01, rate=0, matrix_index=(2, 1)), # mig stop after merge
                      msprime.MigrationRateChange(time=t_01, rate=0, matrix_index=(1, 2)),
                      msprime.MassMigration(time=t_02, source=2, destination=0, proportion=1.0), #next merge
                      msprime.MassMigration(time=t_03, source=3, destination=0, proportion=1.0)] #final merge

现在我们准备模拟树序列。我们将长度设置为 50 kb，重组率为 5e-8。 msprime 速度非常快！

ts = msprime.simulate(population_configurations = population_configurations,
                      migration_matrix = migration_matrix,
                      demographic_events = demographic_events,
                      length = 50000,
                      recombination_rate = 5e-8
                      )

我们可以再次检查树的数量

ts.num_trees

而且，如果我们敢的话，我们可以看看树序列中的第一棵树：

print(ts.first().draw(format="unicode"))

从树序列计算权重

然后运行 twist 中的函数来计算树序列的权重。我们不需要指定组，因为该信息已由 msprime 包含在树序列对象中。但我们仍然需要告诉它使用最终的总体（数字 3，因为 python 从 0 开始计数）。

weightsData = twisst.weightTrees(ts, treeFormat="ts", outgroup = "3", verbose=False)

我们可以快速总结平均权重

twisst.summary(weightsData)

我们还可以直接快速保存权重图（这节省了导出到文件并在 R 中绘制精美图的时间）

#extract mid positions on chromosome from tree sequence file
position = [(tree.interval[0] + tree.interval[1])/2 for tree in ts.trees()]

#normalise weights by dividing by number of combinations
weights = weightsData["weights"]/10000

#create a plot with all three topology weights
for i in range(3): 
    plt.plot(position, weights[:,i], label='topo'+str(i+1))

plt.legend()

#save plot
plt.savefig('sim_ts_weights.pdf')

完！

Reference

[1]

Source: https://github.com/simonhmartin/tutorials/blob/master/topology_weighting/README.md

本文由 mdnice 多平台发布