基因定相（Phasing）与 SHAPEIT 原理简介

基因定相

基因定相（Genotype Phasing、Phasing、Haplotype Phasing、Haplotype Estimation），也称为单倍体分型、单倍体构建等，表示将等位基因定位到父本或者母本染色体上的过程，即将基因型数据转变为单倍型数据的过程。Estimation of haplotypes from genotype data, known as phasing（Delaneau, O., Marchini, J. & Zagury, JF. A linear complexity phasing method for thousands of genomes. Nat Methods 9, 179–181 (2012). https://doi.org/10.1038/nmeth.1785）。

如今大多数 NGS 技术在测序时不区分 reads 的染色体来源，即只能检测出基因组上存在的变异，而无法得知变异位于那条染色体上。所以当序列中存在多个变异时，无法确定群体中单倍型的种类及数量（图片来自 UCLA ZarLab）。单倍型数据的缺失 会造成下游的 基因型推断（Imputation）、连锁不平衡计算、选择清扫区检测、重组率估计 等分析难以进行。

目前，对基因型数据进行定相（Phasing）的方法常为：家系分型（Related individuals Phasing）、物理分型（Physical Phasing）、参考单倍型分型。其中参考单倍型分型一般最为常用。

家系分型 是利用父母本的基因型数据来推断子代的单倍型。家系分型简单、准确，但缺点是：1. 无法应用于父母本信息缺失的群体；2. 父母本均为杂合子时，无法对子代的杂合位点进行定相；3. 需要额外测定样本的父母本基因型，费用大幅提高。（图片来自碱基矿工）

物理分型 一般要求测序时使用的 reads 足够长且测序深度要足够深。因为当 reads 含有多个 SNP时，可以视为一个单倍体的局部单倍型；reads 的拼接成染色体的过程，即为局部单倍型拼接成单倍体的过程。准确的说，物理分型并不是一种分型方法，因为 reads 的序列数据是单倍型数据，而非基因型数据。物理分型的缺点是高深度或三代测序费用较高。

参考单倍型分型 是指利用参考单倍型库（haplotype reference panel，其中 panel 是组、库的意思）中的单倍型来指导基因定相。参考单倍型库有 2002-2007 年的国际人类基因组单体型图谱计划（International HapMap Project），通过高深度的测序，得到了 310 万个 SNP 和 420 个单倍型；2008-2015 年的千人基因组计划（1000 Genomes Project），通过高深度的测序，得到了 8800 万个 SNP 和 5008 个单倍型。参考分型的常用推断方法为 HMM（隐马可夫模型），这里通过简单介绍 SHAPEIT 软件的原理来阐述如何利用 HMM 来进行参考分型。参考分型的应用有：2016 年科研人员利用 5008 个参考单倍型指导 20 个以欧洲人为主的低覆盖度（4-8 X）的全基因组测序项目的基因分型，最终得到了 3923 万个 SNP 和 64976 个单倍型，扩充了参考单倍型库，并构建了人类参考单倍型联盟（Haplotype Reference Consortium，HRC）—— A reference panel of 64,976 haplotypes for genotype imputation 。

SHAPEIT 原理简介

参考单倍型分型 HMM 原理

HMM 进行单倍型分型的 基本假设 ：

待分型个体的单倍型来自于参考库中已有单倍型；
单倍型间的重组或位点的突变是造成待分型个体的单倍型与参考库中已有单倍型间存在差异的原因；
基因型 mmm 位点的单倍型仅与 m−1m-1m−1 位点的单倍型相关；

HMM 与单倍体分型之间的 对应关系 （例子红线路径）：

隐藏状态（ZiZ_iZi）：基因型数据 GGG 在 mmm 位点处 SNP 的参考单倍型。
输出状态（XiX_iXi）：基因型数据 GGG 在 mmm 位点处 SNP 的候选单倍型。
隐藏状态数量（kkk）：参考单倍型库 HHH 中单倍型的数量。如图中位点 Z1Z_1Z1 含有 4 个隐藏状态 H1,H2,H3,H4H_1, H_2, H_3, H_4H1,H2,H3,H4。
输出状态数量（hhh）：本文中讨论的均为双等位型 SNP，单倍体上仅存在 000、111 两种输出状态。
转移（transition）概率：mmm 位点的 SNP 从参考单倍型 kmk_mkm 变化到 m+1m+1m+1 位点 SNP 的参考单倍型 km+1k_{m+1}km+1 的概率。如图中，红线路径内 Z1=H1→Z2=H2Z_1=H_1 \rightarrow Z_2=H_2Z1=H1→Z2=H2 的转移涉及染色体重组 H1→H2H_1 \rightarrow H_2H1→H2，所以概率为 ρ/3\rho/3ρ/3 ，其中 ρ\rhoρ 为重组率；红线路径内 Z2=H2→Z3=H2Z_2=H_2 \rightarrow Z_3=H_2Z2=H2→Z3=H2 的转移不涉及染色体重组 H2→H2H_2 \rightarrow H_2H2→H2，所以概率为 1−ρ1-\rho1−ρ 。
输出（emission）概率：mmm 位点的 SNP 从参考单倍型 kmk_mkm 变化到候选单倍型 hmh_mhm 的概率。如图中，红线路径内 Z1=H1=0→X1=S1=1Z_1=H_1=0 \rightarrow X_1=S_1=1Z1=H1=0→X1=S1=1 的输出涉及突变 0→10 \rightarrow 10→1，所以概率为 Θ\ThetaΘ ，其中 Θ\ThetaΘ 为突变率；红线路径内 Z4=H3=0→X4=S1=0Z_4=H_3=0 \rightarrow X_4=S_1=0Z4=H3=0→X4=S1=0 的输出不涉及突变 0→00 \rightarrow 00→0，所以概率为 1−Θ1-\Theta1−Θ 。
隐藏状态链（z⃗\vec{z}z）：由各位点隐藏状态所组成的序列 z⃗=(Z1,Z2,...,ZM)\vec{z}=(Z_1, Z_2,...,Z_M)z=(Z1,Z2,...,ZM)。如图中 红线路径 即为一种隐藏状态链，z⃗=(H1,H2,H2,H3)\vec{z}=(H_1, H_2, H_2, H_3)z=(H1,H2,H2,H3)
输出状态链（x⃗\vec{x}x）：由各位点输出状态所组成的序列 x⃗=(X1,X2,...,XM)\vec{x}=(X_1, X_2,...,X_M)x=(X1,X2,...,XM)。如图中 红线路径 即为一种输出状态链，x⃗=(S1,S1,S1,S1)\vec{x}=(S_1, S_1, S_1, S_1)x=(S1,S1,S1,S1)。基因型数据 GGG 的候选单倍型库 SSS 罗列了所有符合 GGG 的可能单倍型，即所有可能的输出状态链。

HMM 进行单倍体分型的步骤：

根据基因型 GGG 计算候选单倍型库 SSS，其中每种候选单倍型 SiS_iSi 都是一条输出状态链 x⃗\vec{x}x 。当基因型 GGG 中包含 sss 个杂合位点时，SSS 中将包含 2s2^s2s 个候选单倍型。
利用 HMM ，根据 H,ρ,ΘH, \rho, \ThetaH,ρ,Θ 计算每种候选单倍型 SiS_iSi 的概率 P(Si∣H)P(S_i|H)P(Si∣H)：
P(Si∣H)=P(X1,...,XM∣H)=P(X1∣H)∏m=2MP(Xm∣Xm−1,H)P(S_i|H)=P(X_1,...,X_M|H)=P(X_1|H)\prod_{m=2}^{M} P(X_m|X_{m-1}, H)P(Si∣H)=P(X1,...,XM∣H)=P(X1∣H)m=2∏MP(Xm∣Xm−1,H)
如 P(S1∣H)=P(X1=1∣H)×P(X2=1∣H,X1=1)×P(X3=1∣H,X2=1)×P(X4=0∣H,X3=1)P(S_1|H)=P(X_1=1|H)×P(X_2=1|H, X_1=1)×P(X_3=1|H, X_2=1)×P(X_4=0|H, X_3=1)P(S1∣H)=P(X1=1∣H)×P(X2=1∣H,X1=1)×P(X3=1∣H,X2=1)×P(X4=0∣H,X3=1)
通过归一化使概率合为 111，∑P(Si∣H,G)=1\sum P(S_i|H, G)=1∑P(Si∣H,G)=1 ，得到输出状态链的概率密度分布函数。
根据 SSS 的概率密度分布进行抽样，抽出的 SiS_iSi 被作为 GGG 的 1 个单倍型，再根据 SiS_iSi 和 GGG 推导出另 1 个单倍型，至此完成了对 GGG 的单倍型分型，得到了双倍型 DDD。

HMM 进行单倍体分型的 计算量 ：

如图所示，参考单倍型库 HHH 中包含 4 个单倍型 H1、H2、H3、H4H_1、H_2、H_3、H_4H1、H2、H3、H4，SNP 的基因型用 0、10、10、1 表示。现在有某样本的基因型序列数据 GGG，SNP 用 0、1、20、1、20、1、2 表示，候选单倍型 SSS 有 S1、S2、S3、S4S_1、S_2、S_3、S_4S1、S2、S3、S4，我们将 HHH 中每个 SNP 位点标记为 ZiZ_iZi，SSS 中标记为 XiX_iXi。下面我们利用 HMM 算法根据 HHH 计算基因型数据 GGG 为单倍型 S1S_1S1 的概率 P(S1∣G,H)P(S_1|G, H)P(S1∣G,H)。

图中红线展示了 S1S_1S1 一种可能的隐藏状态链：H1→H2→H2→H3H_1 \rightarrow H_2 \rightarrow H_2 \rightarrow H_3H1→H2→H2→H3

P(S1∣zi⃗=(H1,H2,H2,H3))=14Θ×ρ3Θ×(1−ρ)(1−Θ)×ρ3(1−Θ)P(S_1|\vec{z_i}=(H1, H2, H2, H3)) = \frac{1}{4} \Theta × \frac{\rho}{3}\Theta × (1-\rho)(1-\Theta) × \frac{\rho}{3}(1-\Theta)P(S1∣zi=(H1,H2,H2,H3))=41Θ×3ρΘ×(1−ρ)(1−Θ)×3ρ(1−Θ)

我们可以发现，如果 HHH 库中有 KKK 个参考单倍型，则每个位点 XiX_iXi 都有 KKK 个可能的隐藏状态。所以当基因型 GGG 有 MMM 个位点时，S1S_1S1 总共有 MKM^{K}MK 种可能的隐藏链。若 GGG 中总共有 LLL 个杂合位点，则 SSS 中总共有 2L2^L2L 个候选单倍型（符合 GGG 的输出状态链）且每个单倍型都包含 MKM^{K}MK 种可能的隐藏链，总计算量为：

O=2L×MKO=2^L×M^{K}O=2L×MK

计算完 2L×MK2^L×M^{K}2L×MK 种可能情况的概率后，对已有的概率密度分布进行采样，才能得到 GGG 的双倍型（dipoltype） DDD。

SHAPEIT 加速单倍型推断

SHAPEIT 软件通过压缩隐藏状态数量（参考单倍型库 HHH）和输出状态链的数量（候选单倍型库 SSS）来加速单倍型推断。

压缩参考单倍型库 HHH

上图中 HgH_gHg 为 SHAPEIT 对参考单倍型库 HHH 压缩后的结果：

空心矩阵表示 0，实心矩阵表示 1；
片段内（intra-segment）在 HHH 中相邻的 SNP 之间用实线连接，片段间（inter-segment）用虚线连接；
片段内相同的单倍型会被合并，边 c(km,km+1)c(k_m,k_{m+1})c(km,km+1) 上数字表示单倍型片段在原 HHH 库中的数量，如 c(12,13)=4c(1_2,1_3)=4c(12,13)=4，其中 kkk、mmm 分别表示 HgH_gHg 中的单倍型和 SNP 序号；

如图中 HHH 序列长度 M=8M=8M=8，含有 K=8K=8K=8 种参考单倍型，SHAPEIT 在第 4、5 SNP 间对单倍型进行切割，单倍型被分为 2 个片段（segment），并对片段内相同的单倍型进行压缩。第 1 个片段内 8 种单倍型被压缩成为 3 种，将 Z1→Z4Z_1 \rightarrow Z_4Z1→Z4 在 HHH 中 84=40968^4=409684=4096 种可能隐藏链压缩称为 HgH_gHg 中 34=813^4=8134=81 种计算量大幅降低。SHAPEIT 中根据 HgH_gHg 推导出的起始概率、转移概率、输出概率参见附录。

综上，通过 先分割再压缩，大幅减少单倍型库的 隐藏状态数量，牺牲了重组部分的精度来加速 HMM 的推断（计算公式参见附录）。

压缩候选单倍型库 SSS

∵P(Si∣H)=P(X1,...,XM∣H)=P(X1∣H)∏m=2MP(Xm∣Xm−1,H)\because P(S_i|H)=P(X_1,...,X_M|H)=P(X_1|H)\prod_{m=2}^{M} P(X_m|X_{m-1}, H)∵P(Si∣H)=P(X1,...,XM∣H)=P(X1∣H)∏m=2MP(Xm∣Xm−1,H)

∴\therefore∴ 当候选单倍型之间存在相同片段时，计算是相同的。

S1,S2S_1, S_2S1,S2 中相同片段（X1=1,X2=1,X3=1X_1=1, X_2=1, X_3=1X1=1,X2=1,X3=1）的计算是相同的。

∴\therefore∴ 合并相同的候选单倍型片段，可以大幅减少冗余计算。

SHAPEIT 利用与 HHH 相同的压缩方法压缩 SSS 中候选单倍型的数量。如下图中，基因型 GGG 含有 4 个杂合位点，理论上包含 16 种候选单倍型。SHAPEIT 通过在第 5、6 SNP 间对单倍型进行切割并合并相同的单倍型片段，使片段 1（X1,X2,X3,X4,X5X_1, X_2, X_3, X_4, X_5X1,X2,X3,X4,X5）、片段 2（X6,X7,X8X_6, X_7, X_8X6,X7,X8）内的计算量从 161616 减少至 444，片段间的计算量不变。SHAPEIT 软件默认 1 个片段包含 3 个杂合位点，即 GGG 中每个 片段内 的计算量从 2s2^s2s 缩减到 888。

附录

CHMM 起始状态概率函数 P(z1=k1)P(z_1=k_1)P(z1=k1) ：
P(z1=k1)=c(k1)KP(z_1=k_1)=\frac{c(k_1)}{K}P(z1=k1)=Kc(k1)

CHMM 状态转移概率函数 P(zm+1=km+1∣zm=km)P(z_{m+1}=k_{m+1} | z_m=k_m)P(zm+1=km+1∣zm=km)，其中 ρm\rho_mρm 为 mmm 位点发生重组的概率：
P(zm+1=km+1∣zm=km)=(1−ρm)c(km,km+1)c(km)+ρmc(km+1)KP(z_{m+1}=k_{m+1} | z_m=k_m)=(1-\rho_m)\frac{c(k_m,k_{m+1})}{c(k_m)}+\rho_m\frac{c(k_{m+1})}{K}P(zm+1=km+1∣zm=km)=(1−ρm)c(km)c(km,km+1)+ρmKc(km+1)

如图中紫虚线 P(z5=35∣z4=24)P(z_5=3_5 | z_4=2_4)P(z5=35∣z4=24)：
P(z5=35∣z4=24)=13(1−ρm)+28ρmP(z_5=3_5 | z_4=2_4)=\frac{1}{3}(1-\rho_m)+\frac{2}{8}\rho_mP(z5=35∣z4=24)=31(1−ρm)+82ρm

但紫虚线 P(z5=35∣z4=24)P(z_5=3_5 | z_4=2_4)P(z5=35∣z4=24) 的严谨概率为：
P′(z5=35∣z4=24)=13(1−ρm)+ρm(23×28+13×18){P}'(z_5=3_5 | z_4=2_4)=\frac{1}{3}(1-\rho_m)+\rho_m(\frac{2}{3}×\frac{2}{8} + \frac{1}{3}×\frac{1}{8})P′(z5=35∣z4=24)=31(1−ρm)+ρm(32×82+31×81)

CHMM（PPP）与 HMM（P′{P}'P′ ）之间状态转移概率的差值为：
P(zm+1=km+1∣zm=km)−P′(zm+1=km+1∣zm=km)=ρmc(km,km+1)c(km)c(km,km+1)KP(z_{m+1}=k_{m+1} | z_m=k_m)-{P}'(z_{m+1}=k_{m+1} | z_m=k_m)=\rho_m\frac{c(k_m,k_{m+1})}{c(k_m)}\frac{c(k_m,k_{m+1})}{K}P(zm+1=km+1∣zm=km)−P′(zm+1=km+1∣zm=km)=ρmc(km)c(km,km+1)Kc(km,km+1)

CHMM 相比于 HMM 的状态转移概率误差会随着 c(km,km+1)c(k_m,k_{m+1})c(km,km+1) 值的增加而增加。

在计算出隐藏状态 kmk_mkm 的概率后，还要计算隐藏状态 kmk_mkm 到观察状态 hmh_mhm 之间的 输出概率，其中 Θ\ThetaΘ 为突变概率：
P(hm∣zm=km)=KK+Θδ(hm,km)+Θ2(K+Θ)P(h_m|z_m=k_m)=\frac{K}{K+\Theta}\delta(h_m,k_m)+\frac{\Theta}{2(K+\Theta)}P(hm∣zm=km)=K+ΘKδ(hm,km)+2(K+Θ)Θ

当 hmh_mhm 和 kmk_mkm 的碱基相同时，δ(hm,km)=1\delta(h_m,k_m)=1δ(hm,km)=1 ，否则为 000 。

参考文章

人类基因组的Phasing原理是什么？
A linear complexity phasing method for thousands of genomes
A reference panel of 64,976 haplotypes for genotype imputation
“文献解读 | 最大参考序列集（reference panel）助力基因型推断”
SHAPEIT算法简介
一文搞懂 HMM 隐马尔可夫模型

基因定相（Phasing）与 SHAPEIT 原理简介相关推荐

遗传算法（GA）的原理简介与应用【python实现】
遗传算法(GA)的原理简介与应用[python实现] 算法原理简介与实现算法框架关于GA求解过程中一些结果的讨论 SA-GA总结与体会算法原理简介与实现遗传算法(Genetic Algorit ...
3D绘图过程及原理简介
3D绘图过程及原理简介收藏 Standard Primitives(标准几何体) 在创建命令面板的Geometry(几何体)对象类型中有如下几个次级分类项目: • Standard Primitiv ...
javascript原理_JavaScript程序包管理器工作原理简介
javascript原理 by Shubheksha 通过Shubheksha JavaScript程序包管理器工作原理简介 (An introduction to how JavaScript pa ...
Nginx 反向代理工作原理简介与配置详解
Nginx 反向代理工作原理简介与配置详解测试环境 CentOS 6.8-x86_64 nginx-1.10.0 下载地址:http://nginx.org/en/download.html 安装 ...
DeepLearning tutorial（1）Softmax回归原理简介+代码详解
FROM: http://blog.csdn.net/u012162613/article/details/43157801 DeepLearning tutorial(1)Softmax回归原理简介 ...
DeepLearning tutorial（3）MLP多层感知机原理简介+代码详解
FROM:http://blog.csdn.net/u012162613/article/details/43221829 @author:wepon @blog:http://blog.csdn.n ...
DeepLearning tutorial（4）CNN卷积神经网络原理简介+代码详解
FROM: http://blog.csdn.net/u012162613/article/details/43225445 DeepLearning tutorial(4)CNN卷积神经网络原理简介 ...
【Android 异步操作】Handler ( 主线程中的 Handler 与 Looper | Handler 原理简介 )
文章目录一.主线程中的 Handler 与 Looper 二.Handler 原理简介一.主线程中的 Handler 与 Looper Android 系统中 , 点击图标启动一个应用进程 , 就 ...
量子计算机编程原理简介和机器学习
量子计算机编程原理简介和机器学习本文翻译自D-Wave公司网站 www.dwavesys.com/en/dev-tutorial-intro.html D-wave公司在2007年就声称实现了1 ...
DL之CNN：卷积神经网络算法简介之原理简介——CNN网络的3D可视化(LeNet-5为例可视化)
DL之CNN:卷积神经网络算法简介之原理简介--CNN网络的3D可视化(LeNet-5为例可视化) CNN网络的3D可视化 3D可视化地址:http://scs.ryerson.ca/~aharley ...

基因定相（Phasing）与 SHAPEIT 原理简介

基因定相