基因定相（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 进行单倍型分型的 基本假设 ：

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

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

1. 隐藏状态（$Z_i$）：基因型数据 $G$ 在 $m$ 位点处 SNP 的参考单倍型。
2. 输出状态（$X_i$）：基因型数据 $G$ 在 $m$ 位点处 SNP 的候选单倍型。
3. 隐藏状态数量（$k$）：参考单倍型库 $H$ 中单倍型的数量。如图中位点 $Z_1$ 含有 4 个隐藏状态 $H_1,H_2,H_3,H_4$。
4. 输出状态数量（$h$）：本文中讨论的均为双等位型 SNP，单倍体上仅存在 $0$、$1$ 两种输出状态。
5. 转移（transition）概率：$m$ 位点的 SNP 从参考单倍型 $k_m$ 变化到 $m+1$ 位点 SNP 的参考单倍型 $k_{m+1}$ 的概率。如图中，红线路径内 $Z_1=H_1\to Z_2=H_2$ 的转移涉及染色体重组 $H_1\to H_2$，所以概率为 $\rho /3$ ，其中 $\rho$ 为重组率；红线路径内 $Z_2=H_2\to Z_3=H_2$ 的转移不涉及染色体重组 $H_2\to H_2$，所以概率为 $1-\rho$ 。
6. 输出（emission）概率：$m$ 位点的 SNP 从参考单倍型 $k_m$ 变化到候选单倍型 $h_m$ 的概率。如图中，红线路径内 $Z_1=H_1=0\to X_1=S_1=1$ 的输出涉及突变 $0\to 1$，所以概率为 $\Theta$ ，其中 $\Theta$ 为突变率；红线路径内 $Z_4=H_3=0\to X_4=S_1=0$ 的输出不涉及突变 $0\to 0$，所以概率为 $1-\Theta$ 。
7. 隐藏状态链（$\vec{z}$）：由各位点隐藏状态所组成的序列 $\vec{z}=(Z_1,Z_2,...,Z_M)$。如图中 红线路径 即为一种隐藏状态链，$\vec{z}=(H_1,H_2,H_2,H_3)$
8. 输出状态链（$\vec{x}$）：由各位点输出状态所组成的序列 $\vec{x}=(X_1,X_2,...,X_M)$。如图中 红线路径 即为一种输出状态链，$\vec{x}=(S_1,S_1,S_1,S_1)$。基因型数据 $G$ 的候选单倍型库 $S$ 罗列了所有符合 $G$ 的可能单倍型，即所有可能的输出状态链。

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

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

2. 利用 HMM ，根据 $H,\rho ,\Theta$ 计算每种候选单倍型 $S_i$ 的概率 $P(S_i|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(S_1|H)=P(X_1=1|H)\times P(X_2=1|H,X_1=1)\times P(X_3=1|H,X_2=1)\times P(X_4=0|H,X_3=1)$

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

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

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

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

图中红线展示了 $S_1$ 一种可能的隐藏状态链：$H_1\to H_2\to H_2\to H_3$

$$P(S_1|\overrightarrow{z_i}=(H1,H2,H2,H3))=\frac{1}{4}\Theta \times \frac{\rho }{3}\Theta \times (1-\rho )(1-\Theta )\times \frac{\rho }{3}(1-\Theta )$$

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

$$O=2^L\times M^K$$

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

SHAPEIT 加速单倍型推断

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

压缩参考单倍型库 $H$

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

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

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

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

压缩候选单倍型库 $S$

$\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)$

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

如图中 $S_1$，$P(S_1|H)=P(X_1=1|H)\times P(X_2=1|H,X_1=1)\times P(X_3=1|H,X_2=1)\times P(X_4=0|H,X_3=1)$
如图中 $S_2$，$P(S_2|H)=P(X_1=1|H)\times P(X_2=1|H,X_1=1)\times P(X_3=1|H,X_2=1)\times P(X_4=1|H,X_3=1)$

$S_1,S_2$ 中相同片段（$X_1=1,X_2=1,X_3=1$）的计算是相同的。

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

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

附录

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

CHMM 状态转移概率函数 $P(z_{m+1}=k_{m+1}|z_m=k_m)$，其中 ${\rho }_m$ 为 $m$ 位点发生重组的概率：
$$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(z_5=3_5|z_4=2_4)$：
$$P(z_5=3_5|z_4=2_4)=\frac{1}{3}(1-{\rho }_m)+\frac{2}{8}{\rho }_m$$

但 紫虚线 $P(z_5=3_5|z_4=2_4)$ 的严谨概率为：
$$P^{\prime }(z_5=3_5|z_4=2_4)=\frac{1}{3}(1-{\rho }_m)+{\rho }_m(\frac{2}{3}\times \frac{2}{8}+\frac{1}{3}\times \frac{1}{8})$$

CHMM（$P$）与 HMM（$P^{\prime }$ ）之间状态转移概率的差值为：
$$P(z_{m+1}=k_{m+1}|z_m=k_m)-P^{\prime }(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}$$

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

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

当 $h_m$ 和 $k_m$ 的碱基相同时，$\delta (h_m,k_m)=1$ ，否则为 $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 隐马尔可夫模型