找外包公司做网站,湘潭网站建设 搜索磐石网络,网站制作的电话,免费邯郸网站建设使用 Twisst 探索整个基因组的进化关系的拓扑加权教程[1]。 简介 拓扑加权是量化不一定是单系群之间关系的一种方法。它通过考虑更简单的“分类单元拓扑”并量化与每个分类单元拓扑匹配的子树的比例#xff0c;提供了复杂谱系的摘要。我们用来计算权重的方法称为 Twisst#… 使用 Twisst 探索整个基因组的进化关系的拓扑加权教程[1]。 简介 拓扑加权是量化不一定是单系群之间关系的一种方法。它通过考虑更简单的“分类单元拓扑”并量化与每个分类单元拓扑匹配的子树的比例提供了复杂谱系的摘要。我们用来计算权重的方法称为 Twisst通过子树迭代采样进行拓扑权重。 在本次实践中我们将使用模拟数据来探索拓扑权重如何提供谱系历史。然后我们将尝试使用针对窄窗口推断的邻居连接树来推断整个模拟染色体的拓扑权重。 工作流程 我们将通过使用 msprime 运行合并模拟然后直接从输出计算拓扑权重来探索人口统计参数如何影响权重所有这些都在 Python 中进行。 使用树序列格式的拓扑加权 到目前为止我们已经分析了树文件这些文件对于具有独特谱系的每个染色体“块”具有独特的树或者在推断树的情况下对于每个窗口具有独特的树。这种格式有点浪费因为相邻的树通常非常相似通常仅因单个重组事件而不同该重组事件将一个分支从树中的一个点移动到另一个点。 树序列格式非常高效因为它不仅记录树上节点之间的连接还记录每个连接所在的染色体的长度。 msprime 和 tskit 包使用此格式它还提供了更多信息。 模拟树序列 我们将使用 msprime 来模拟树序列。 msprime 是一个合并模拟器这意味着它的工作原理是计算任意两个个体在过去某个特定时间拥有共同祖先的概率。在单一种群中这是由种群规模决定的。对于多个种群这还受到种群之间的迁移率以及它们从单一祖先种群传承多久的影响。 我们将使用 Python 交互式环境来处理 msprime 和树序列并使用 twisst 中的函数来分析它。 为了确保我们可以从 python 中导入 twisst 模块请将其添加到我们的 python 路径中。 export PYTHONPATH$PYTHONPATH:twisst-0.2 现在打开一个Python交互式会话输入“python”并在文本编辑器中打开一个脚本因为我们将多次修改并重新运行其中一些行以测试不同模拟参数的效果 导入所需的模块。 import msprimeimport matplotlib.pyplot as pltimport twisst 我们将从对单个总体的 10 个样本进行简单模拟开始。我们必须指定长度和重组率因此 msprime 将为我们提供一个由多个谱系组成的序列通过重组分开。在这里我们还指定了随机种子只是为了确保在这种情况下我们都得到相同的模拟。 ts msprime.simulate(sample_size10, Ne1000, length10000, recombination_rate5e-8, random_seed 1) 我们可以检查树序列中有多少个不同的家谱。 ts.num_trees 并使用树序列对象提供的良好可视化方法来查看它们。 for tree in ts.trees(): print(interval , tree.interval) print(tree.draw(formatunicode)) 包含四个种群和基因流的模拟 我们现在将建立一个包含多个群体的更大的模拟。这需要一些不同的组件我们将分别定义它们。首先我们定义样本数量和每个总体的总体规模。 pop_n 10pop_Ne 10000population_configurations [msprime.PopulationConfiguration(sample_sizepop_n, initial_sizepop_Ne), msprime.PopulationConfiguration(sample_sizepop_n, initial_sizepop_Ne), msprime.PopulationConfiguration(sample_sizepop_n, initial_sizepop_Ne), msprime.PopulationConfiguration(sample_sizepop_n, initial_sizepop_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 1000t_02 5000t_03 10000demographic_events [msprime.MassMigration(timet_01, source1, destination0, proportion1.0), # first merge msprime.MigrationRateChange(timet_01, rate0, matrix_index(2, 1)), # mig stop after merge msprime.MigrationRateChange(timet_01, rate0, matrix_index(1, 2)), msprime.MassMigration(timet_02, source2, destination0, proportion1.0), #next merge msprime.MassMigration(timet_03, source3, destination0, proportion1.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(formatunicode)) 从树序列计算权重 然后运行 twist 中的函数来计算树序列的权重。我们不需要指定组因为该信息已由 msprime 包含在树序列对象中。但我们仍然需要告诉它使用最终的总体数字 3因为 python 从 0 开始计数。 weightsData twisst.weightTrees(ts, treeFormatts, outgroup 3, verboseFalse) 我们可以快速总结平均权重 twisst.summary(weightsData) 我们还可以直接快速保存权重图这节省了导出到文件并在 R 中绘制精美图的时间 #extract mid positions on chromosome from tree sequence fileposition [(tree.interval[0] tree.interval[1])/2 for tree in ts.trees()]#normalise weights by dividing by number of combinationsweights weightsData[weights]/10000#create a plot with all three topology weightsfor i in range(3): plt.plot(position, weights[:,i], labeltopostr(i1))plt.legend()#save plotplt.savefig(sim_ts_weights.pdf) 完 Reference [1] Source: https://github.com/simonhmartin/tutorials/blob/master/topology_weighting/README.md 本文由 mdnice 多平台发布
