基因与医疗和健康

 找回密码
 立即注册

扫一扫,访问微社区


这是面向专业人士的交流平台,以基因科技发展及其转化应用为核心,连接上下游相关方,倡导有序交流和知识共享,沉淀精华,增进相互了解与互信,提高科技转化应用效率,推动生命科学相关行业的发展。本论坛定位于专业人员交流,不做患者和公众咨询。

查看: 1611|回复: 1

遗传病测序数据分析解读科普系列 第3.1期 变异发现细节:...

[复制链接]

31

主题

23

帖子

357

积分

版主

Rank: 7Rank: 7Rank: 7

积分
357
发表于 2018-9-2 10:35 | 显示全部楼层 |阅读模式
本帖最后由 孙继国 于 2018-9-2 10:43 编辑

遗传病测序数据分析解读科普系列 第3.1期 变异发现细节:走进HaplotypeCaller内部

一、背景介绍
上期文章介绍了HaplotypeCaller的过程后,发现大多数读者还是愿意听四个步骤的详细解释,那么本期内容就来讲讲每个步骤背后的原理,由于涉及的数学原理较多,我也很难保证每个点都理解的很到位,因此,如果有错误的地方,还请读者指出。


二、简单回顾
HaplotypeCaller的主要步骤
       定义活跃区域(寻找要做组装的区域)
       组装活跃区域并分析单倍体型(确定哪些片段能合并)
       确定单倍体型的可能性(合并片段的可靠性)
       确定分析样本的变异信息(生成变异文件)
图片1.jpg


三、定义活跃区域 ActiveRegion determination
图片2.jpg
主要目的
上期内容讲过,HaplotypeCaller的核心主要是通过类似搭桥的局部组装过程,这个过程通常非常耗费时间,而且很多区域也没有做局部组装的必要(如组装marker非常稀少的区域),因此,如果测序的每个区域都做搭桥计算,在时间和计算成本上都会非常不划算,因此,需要需要一个预处理来决定哪些区域做组装处理。
这个步骤就是:
定义活跃区域(ActiveRegion determination)

实现步骤
确定变异坐标(Calculating the raw activity profile):这一步的目的是先确定样本可能存在的变异的位置并记录成一个profile文件,为分析活跃区域做好准备。
对这些变异位置做平滑处理(Smoothing the activity profile):分析每个变异和邻近的变异的关系,即通过特定公式计算变异之间的距离值,让变异之间不再是孤立关系。
根据上一步得到的每个变异的距离值确定活跃区域:如果距离值大于一定数值,说明这个区域可以作为组装的区域,反之则丢弃。
图片3.jpg
计算距离使用了gaussian smoothing模型,简单的说,就是将区域内数据做正态分布处理,然后用概率评估变异之间的距离。


四、组装活跃区域并分析单倍体型 Local re-assembly and haplotype determination
图片4.jpg
局部组装过程
接触过基因组从头组装的读者肯定听说过目前最常用的directed DeBruijn graph组装算法,局部组装过程虽然也使用了directed DeBruijn graph,但与从头组装不同的是,HaplotypeCaller的局部组装将局部基因组序列做了DeBruijn graph处理。
图片5.jpg
DeBruijn graph图示

图片6.jpg
局部基因组DeBruijn graph处理

局部基因组经过DeBruijn graph处理后,将比对到这个区域的reads做kmer化后,和基因组DeBruijn graph做比较,当发现不一致的kmer序列时,就在graph上生成一个新的分支,并做局部比对,增加indel的准确性,直到所有的reads被处理完,得到一个新的DeBruijn graph。
得到新的DeBruijn graph后寻找所有的欧拉路径,即可能的单倍体型,根据可信度留下较为可信的单倍体型。
图片7.jpg

如何理解这个过程
DeBruijn graph的概念可能对于大多数读者都是陌生的,可以理解为:
基因组从头到尾是一条固定的路,每个碱基就是一座检查站,测序片段可以比作需要在这条路上通过的车,当碱基与对应基因组一致时,顺利通过检查站,当序列和基因组不一样时,是通不过检查站的。
因此只能绕开这个检查站,在这个检查站旁边留下自己的碱基记录,每辆车的行驶路径,就是测序片段对应的单倍体型。


五、确定单倍体型的可能性 Evaluating the evidence for haplotypes and variant alleles
图片8.jpg
这个步骤主要是通过pairHMM模型来计算每种haplotype的可能性,由于测序深度和可能的测序错误的原因,上一步中分析的haplotype中可能会有假阳性出现,因此,需要一个判断haplotype的步骤,这个步骤就是pairHMM。

pairHMM如何工作
pairHMM基于HMM,即隐马氏模型,不同的是,由于比较的是测序片段和基因组之间的关系,是序列pair的比较,通过HMM判断变异的连续状态,即单倍体型,因此叫做pairHMM。
图片9.jpg
这个过程,大家只要知道,通过pairHMM模型,就可以计算出特定单倍体型的概率,通过概率来判断单倍体型的可信度。


六、确定分析样本的变异信息 Assigning per-sample genotypes
图片10.jpg
这个步骤通过贝叶斯统计的方法,确定每个变异位置的基因型,贝叶斯的思想是结合假设的先验概率和观测到的数据,通过贝叶斯公式来推断后验概率。
回到基因型判断问题:我们假设的先验概率是基因组指定位置的特定基因型出现的概率P(G), 我们观测到的是比对到这个位置的测序片段碱基数据D,因此,需要根据这两个信息来计算出对应的基因型P(G|D)。对应贝叶斯公式为
P(G|D)~P(G)P(D|G)
注:“~”代表简化的贝叶斯公式,因为分母P(D)为固定值

因此,通过计算这个公式,便可以得到每种基因型的概率,最终选择概率最大的基因型记录在vcf文件中。
需要注意的是,贝叶斯公式中的P(D|G)的实际计算过程比较复杂,需要计算每条测序片段对应的假设单倍体型,这里不多做讲述,理解目的即可。


七、小结
从文中可以感受到HaplotypeCaller的过程是比较复杂的,因此,HaplotypeCaller也是GATK流程中比较耗费时间的一步。
HaplotypeCaller每一步都包含不同的数学模型,由于笔者不是完全熟悉所有的模型,因此如果有认为写的不清楚或错误之处,还请讨论指出。




上一篇:遗传病测序数据分析解读科普系列 第3期 变异发现之旅:G...
下一篇:遗传病测序数据分析解读科普系列 第4期 群体的力量大:...

91

主题

339

帖子

2697

积分

管理员

双双娘+宝咪阿姨

Rank: 9Rank: 9Rank: 9

积分
2697

金牌特邀

QQ
发表于 2018-9-3 15:03 | 显示全部楼层
HaplotypeCaller博大精深,很值得去仔细理解,以免照猫画虎。很感谢您的科普!
不要去吐槽。觉得对方做的不到位,去补充去帮助提高。
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋| ( 京ICP备18008006号-2

免责声明:用户发帖仅代表个人意见,不代表本论坛立场。 GMT+8, 2019-12-10 09:57

Powered by 顾大夫®工作室  技术支持:北京众易科技

© 2018 基因与医疗和健康

快速回复 返回顶部 返回列表