生信教程:使用拓扑加权探索基因组进化(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-40],
                    [0,    1e-40,    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=(21)), # mig stop after merge
                      msprime.MigrationRateChange(time=t_01, rate=0, matrix_index=(12)),
                      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 多平台发布

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/105331.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

基于nodejs+vue教学辅助管理系统

学生;首页、个人中心、本课程设计了线上教学辅助系统 ,学生可以此系统实现在线学习,作业提交管理、作业成绩管理。随着社会的快速发展,计算机的影响是全面且深入的。教师:首页、个人中心、课程信息管理、教学资料管理、作业信息管…

IDEA的使用(四)创建不同类型的工程(IntelliJ IDEA 2022.1.3版本)

1. 创建Java工程 创建之后,src下是空的。可以在src下创建软件包Package,命名采用域名倒序。在软件包下再创建Java类。Java类运行后出现中文乱码,就到控制台和文件编码这两个地方设置编码。 2. 创建JavaWeb工程 2.1 在win11和IDEA中配置Tomca…

《UnityShader入门精要》学习4

一个最简单的顶点/片元着色器 一个最简单的顶点/片元着色器 Unity Shader的基本结构。它包含了Shader、Properties、SubShader、Fallback等语义块。顶点/片元着色器的结构与之大体类似 Shader "MyShaderName" {Properties {// 属性}SubShader {// 针对显卡A的S…

React之事件机制与事件绑定

一,时间机制 #是什么 React基于浏览器的事件机制自身实现了一套事件机制,包括事件注册、事件的合成、事件冒泡、事件派发等 在React中这套事件机制被称之为合成事件 #合成事件(SyntheticEvent) 合成事件是 React模拟原生 DOM…

opencv图形绘制2

目录 制作宣传语(中文) 制作宣传语(英文) 绘制标记 鼠标交互绘制十字线 鼠标交互绘制图形 鼠标交互制作几何画板 滚动条控制 鼠标事件练习 制作宣传语(中文) import cv2 import numpy as np from …

React之组件通信

#一、是什么 我们将组件间通信可以拆分为两个词: 组件通信 回顾Vue系列 (opens new window)的文章,组件是vue中最强大的功能之一,同样组件化是React的核心思想 相比vue,React的组件更加灵活和多样,按照不同的方式可…

实现Element Select选择器滚动加载

<template><el-selectpopper-class"more-tag-data"v-model"tagId"filterableplaceholder"请选择"focus"focusTag"><el-optionv-for"(item, index) in taskTagLists":key"index":label"item.n…

互联网Java工程师面试题·Java 并发编程篇·第七弹

目录 16、CAS 的问题 17、什么是 Future&#xff1f; 18、什么是 AQS 19、AQS 支持两种同步方式&#xff1a; 20、ReadWriteLock 是什么 21、FutureTask 是什么 22、synchronized 和 ReentrantLock 的区别 23、什么是乐观锁和悲观锁 24、线程 B 怎么知道线程 A 修改了…

蓝桥杯 枚举算法 (c++)

枚举就是根据提出的问题&#xff0c;——列出该问题的所有可能的解&#xff0c;并在逐一列出的过程中&#xff0c;检验每个可能解是否是问题的真正解&#xff0c; 如果是就采纳这个解&#xff0c;如果不是就继续判断下一个。 枚举法一般比较直观&#xff0c;容易理解&#xff0…

完美解决lftp遇到put: Access failed: 553 Could not create file.

目录 一、问题 二、原因 三、解决方法 一、问题 put: Access failed: 553 Could not create file. 二、原因 &#xff08;1&#xff09;没有关闭SeLinux &#xff08;2&#xff09;linux默认安装vsftp服务之后只允许匿名用户的访问和下载&#xff0c;不支持上传。 三、解决方…

源码解析FlinkKafkaConsumer支持punctuated水位线发送

背景 FlinkKafkaConsumer支持当收到某个kafka分区中的某条记录时发送水位线&#xff0c;比如这条特殊的记录代表一个完整记录的结束等&#xff0c;本文就来解析下发送punctuated水位线的源码 punctuated 水位线发送源码解析 1.首先KafkaFetcher中的runFetchLoop方法 public…

matlab 图像均值滤波

目录 一、算法原理二、代码实现三、结果展示本文由CSDN点云侠翻译,放入付费专栏只为防不要脸的爬虫。专栏值钱的不是本文,切勿因本文而订阅。 一、算法原理 均值滤波是一种常用的线性滤波方法,用于平滑图像并减少噪声。它的实现过程如下: 确定滤波器的大小:选择一个固定的…

flutter dio 请求封装(空安全)

一、添加依赖 dio: ^5.3.2二、请求封装 class HttpHelper {static Dio? mDio;static BaseOptions? options;static HttpHelper? httpHelper;CancelToken cancelToken CancelToken();static const String GET get;static const String POST post;static const String PU…

js冒泡排序的几种写法?

冒泡排序是一种简单的排序算法&#xff0c;其基本思想是重复地遍历要排序的数列&#xff0c;一次比较两个元素&#xff0c;如果他们的顺序错误就把他们交换过来。以下是几种常见的 JavaScript 实现方式&#xff1a; 使用基本的 for 循环实现冒泡排序&#xff1a; function bub…

P1443 马的遍历

#include <iostream> #include <queue> using namespace std; #define M 400 int arr[M 5][M 5]; typedef struct Node {int x, y; } Node; //将马能走的8个方向封装成一个二维数组 int dir[8][2] {{2, 1}, {2, -1}, {-2, 1}, {-2, -1},{1, 2}, {-1, 2}, {1, -2…

nginx的location的优先级和匹配方式

nginx的location的优先级和匹配方式 在http模块中有server&#xff0c;server模块中有location&#xff0c;location匹配的是uri 在一个server中&#xff0c;会有多个location&#xff0c;如何来确定匹配哪个location niginx的正则表达式 ^ 字符串的起始位置 $ 字符串的…

IDEA中查看整个项目代码行数

近期正在手撸Spring源码&#xff0c;想要看下自己写了多少行代码。记录下如何查看项目的代码行数&#xff0c;方便日后查阅

Hadoop3教程(六):HDFS中的DataNode

文章目录 &#xff08;63&#xff09;DataNode工作机制&#xff08;64&#xff09;数据完整性&#xff08;65&#xff09;掉线时限参数设置参考文献 &#xff08;63&#xff09;DataNode工作机制 DataNode内部存储了一个又一个Block&#xff0c;每个block由数据和数据元数据组…

【云计算】相关解决方案介绍

文章目录 1.1 云服务环境 Eucalyptus1.1.1 介绍1.1.2 开源协议及语言1.1.3 官方网站 1.2 开源云计算平台 abiCloud1.2.1 开源协议及语言1.2.2 官方网站 1.3 分布式文件系统 Hadoop1.3.1 开源协议及语言1.3.2 官方网站 1.4 JBoss云计算项目集 StormGrind1.4.1 开源协议及语言1.4…

【数据库——MySQL(实战项目1)】(3)图书借阅系统——存储函数

目录 1. 简述2. 功能代码2.1 创建存储函数&#xff0c;根据图书编号查借阅人姓名&#xff0c;并调用该函数查询‘ **小邓在森林** ’已借未还的图书情况&#xff1b;2.2 创建存储函数&#xff0c;计算某借阅人还能借阅的图书数目&#xff0c;学生限额 5 本&#xff0c;教师限额…