基因定相（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 。

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

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

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

根据基因型 G G G 计算候选单倍型库 S S S，其中每种候选单倍型 S i S_i Si​ 都是一条输出状态链 x ⃗ \vec{x} x 。当基因型 G G G 中包含 s s s 个杂合位点时， S S S 中将包含 2 s 2^s 2s 个候选单倍型。

利用 HMM ，根据 H , ρ , Θ H, \rho, \Theta H,ρ,Θ 计算每种候选单倍型 S i S_i Si​ 的概率 P ( S i ∣ H ) P(S_i|H) P(Si​∣H)： P ( S i ∣ H ) = P ( X 1 , . . . , X M ∣ H ) = P ( X 1 ∣ H ) ∏ m = 2 M P ( X m ∣ X m − 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∏M​P(Xm​∣Xm−1​,H) 如 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(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)

通过归一化使概率合为 1 1 1， ∑ P ( S i ∣ H , G ) = 1 \sum P(S_i|H, G)=1 ∑P(Si​∣H,G)=1 ，得到输出状态链的概率密度分布函数。

根据 S S S 的概率密度分布进行抽样，抽出的 S i S_i Si​ 被作为 G G G 的 1 个单倍型，再根据 S i S_i Si​ 和 G G G 推导出另 1 个单倍型，至此完成了对 G G G 的单倍型分型，得到了双倍型 D D D。

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

如图所示，参考单倍型库 H H H 中包含 4 个单倍型 H 1 、 H 2 、 H 3 、 H 4 H_1、H_2、H_3、H_4 H1​、H2​、H3​、H4​，SNP 的基因型用 0 、 1 0、1 0、1 表示。现在有某样本的基因型序列数据 G G G，SNP 用 0 、 1 、 2 0、1、2 0、1、2 表示，候选单倍型 S S S 有 S 1 、 S 2 、 S 3 、 S 4 S_1、S_2、S_3、S_4 S1​、S2​、S3​、S4​，我们将 H H H 中每个 SNP 位点标记为 Z i Z_i Zi​， S S S 中标记为 X i X_i Xi​。下面我们利用 HMM 算法根据 H H H 计算基因型数据 G G G 为单倍型 S 1 S_1 S1​ 的概率 P ( S 1 ∣ G , H ) P(S_1|G, H) P(S1​∣G,H)。

图中红线展示了 S 1 S_1 S1​ 一种可能的隐藏状态链： H 1 → H 2 → H 2 → H 3 H_1 \rightarrow H_2 \rightarrow H_2 \rightarrow H_3 H1​→H2​→H2​→H3​

P ( S 1 ∣ z i ⃗ = ( H 1 , H 2 , H 2 , H 3 ) ) = 1 4 Θ × ρ 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−Θ)

我们可以发现，如果 H H H 库中有 K K K 个参考单倍型，则每个位点 X i X_i Xi​ 都有 K K K 个可能的隐藏状态。所以当基因型 G G G 有 M M M 个位点时， S 1 S_1 S1​ 总共有 M K M^{K} MK 种可能的隐藏链。若 G G G 中总共有 L L L 个杂合位点，则 S S S 中总共有 2 L 2^L 2L 个候选单倍型（符合 G G G 的输出状态链）且每个单倍型都包含 M K M^{K} MK 种可能的隐藏链，总计算量为：

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

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

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

上图中 H g H_g Hg​ 为 SHAPEIT 对参考单倍型库 H H H 压缩 后的结果：

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

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

∵ P ( S i ∣ H ) = P ( X 1 , . . . , X M ∣ H ) = P ( X 1 ∣ H ) ∏ m = 2 M P ( X m ∣ X m − 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=2M​P(Xm​∣Xm−1​,H)

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

如图中 S 1 S_1 S1​， 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(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) 如图中 S 2 S_2 S2​， P ( S 2 ∣ 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 = 1 ∣ H , X 3 = 1 ) P(S_2|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=1|H, X_3=1) P(S2​∣H)=P(X1​=1∣H)×P(X2​=1∣H,X1​=1)×P(X3​=1∣H,X2​=1)×P(X4​=1∣H,X3​=1)

S 1 , S 2 S_1, S_2 S1​,S2​ 中相同片段（ X 1 = 1 , X 2 = 1 , X 3 = 1 X_1=1, X_2=1, X_3=1 X1​=1,X2​=1,X3​=1）的计算是相同的。

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

SHAPEIT 利用与 H H H 相同的压缩方法压缩 S S S 中候选单倍型的数量。如下图中，基因型 G G G 含有 4 个杂合位点，理论上包含 16 种候选单倍型。SHAPEIT 通过在第 5、6 SNP 间对单倍型进行切割并合并相同的单倍型片段，使片段 1（ X 1 , X 2 , X 3 , X 4 , X 5 X_1, X_2, X_3, X_4, X_5 X1​,X2​,X3​,X4​,X5​）、片段 2（ X 6 , X 7 , X 8 X_6, X_7, X_8 X6​,X7​,X8​）内的计算量从 16 16 16 减少至 4 4 4，片段间的计算量不变。SHAPEIT 软件默认 1 个片段包含 3 个杂合位点，即 G G G 中每个 片段内 的计算量从 2 s 2^s 2s 缩减到 8 8 8。

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

CHMM 状态转移概率函数 P ( z m + 1 = k m + 1 ∣ z m = k m ) P(z_{m+1}=k_{m+1} | z_m=k_m) P(zm+1​=km+1​∣zm​=km​)，其中 ρ m \rho_m ρm​ 为 m m m 位点发生重组的概率： P ( z m + 1 = k m + 1 ∣ z m = k m ) = ( 1 − ρ m ) c ( k m , k m + 1 ) c ( k m ) + ρ m c ( k m + 1 ) K P(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​)​+ρm​Kc(km+1​)​

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

但 紫虚线 P ( z 5 = 3 5 ∣ z 4 = 2 4 ) P(z_5=3_5 | z_4=2_4) P(z5​=35​∣z4​=24​) 的严谨概率为： P ′ ( z 5 = 3 5 ∣ z 4 = 2 4 ) = 1 3 ( 1 − ρ m ) + ρ m ( 2 3 × 2 8 + 1 3 × 1 8 ) {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（ P P P）与 HMM（ P ′ {P}' P′ ）之间状态转移概率的差值为： P ( z m + 1 = k m + 1 ∣ z m = k m ) − P ′ ( z m + 1 = k m + 1 ∣ z m = k m ) = ρ m c ( k m , k m + 1 ) c ( k m ) c ( k m , k m + 1 ) K P(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​)=ρm​c(km​)c(km​,km+1​)​Kc(km​,km+1​)​

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

在计算出隐藏状态 k m k_m km​ 的概率后，还要计算隐藏状态 k m k_m km​ 到观察状态 h m h_m hm​ 之间的 输出概率，其中 Θ \Theta Θ 为突变概率 ： P ( h m ∣ z m = k m ) = K K + Θ δ ( h m , k m ) + Θ 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+Θ)Θ​

当 h m h_m hm​ 和 k m k_m km​ 的碱基相同时， δ ( h m , k m ) = 1 \delta(h_m,k_m)=1 δ(hm​,km​)=1 ，否则为 0 0 0 。

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