S/HIC 系列软件：diploS/HIC 利用 CNN 和非定向基因型数据识别 软/硬 清扫

本文内容摘取自 论文：diploS/HIC: An Updated Approach to Classifying Selective Sweeps. Kern AD, Schrider DR. G3 (Bethesda). 2018;8(6):1959-1970. Published 2018 May 31. doi:10.1534/g3.118.200262，引用次数 63。

diploS/HIC 与 S/HIC 算法区别

文章提出了新的清扫区识别算法 diploS/HIC，相比于原先的 S/HIC 算法：

1. diploS/HIC 使用 卷积神经网络（CNN） 作为模型框架，替换了 S/HIC 中的 随机森林（RF） 框架。
2. diploS/HIC 从 多个尺度 上捕获不同参数、窗口之间的组合，挖掘出更有效的信息来提高预测精度。
3. S/HIC 的输入特征包含 定相基因型数据（phased） 计算出的群体遗传统计参数，如 单倍型 等参数，但 diploS/HIC 的输入特征仅使用 非定相基因型数据（unphased） 计算出的群体遗传统计参数，扩大了模型的应用范围。
4. 由参考单倍型推测出的单倍型数据，可能会由于错误率较高（error rate>0.05）而使定相数据训练出的模型精度较低，作者推荐使用非定相数据作为特征训练模型。

PS：（非）定相基因型数据参见：基因定相（Phasing） 与 SHAPEIT 原理简介

数据集

模拟数据集（Coalescent simulation）

窗口数量：$11$
模拟次数（训练集）：$2000$
模拟次数（测试集）：$1000$

平衡群体

窗口大小：$L=110KB$

突变率：$4N_0\mu L=U(528,1584)$
重组率：$4N_{0}rL=880$
选择强度：$2N_{0}s=U(25,2.5\times 10^2)$、$U(2.5\times 10^2,2.5\times 10^3)$、$U(2.5\times 10^3,2.5\times 10^4)$
软清扫起始频率：$f=U(0,0.2)$
固定与观察间的时间差：$\tau =U(0,0.05)$

群体大小：保持不变

非平衡群体（模拟 Ag1000G）

窗口大小：$L=55KB$

突变率：$4N_0\mu L=U(1750,17500)$
重组率：$4N_{0}rL=U(1925,57756)$
选择强度：$2N_{0}s=U(2.5\times 10^2,2.5\times 10^3)$、$U(2.5\times 10^3,2.5\times 10^4)$、$U(2.5\times 10^4,2.5\times 10^5)$
软清扫起始频率：$f=U(0,0.2)$
固定与观察间的时间差：$\tau =U(0,0.00004)$

群体大小：先增加至 102%，再减少至 12%，再增加至 21%，再减少至 1%，再增加至 26%，再减少至 0.2%，再增加至 4%

真实数据集（Ag1000G）

真实数据集为 Anopheles gambiae 1000 Genomes Project (Ag1000G) ，其中包含 Anopheles gambiae 和 Anopheles coluzzii 两种疟蚊，这两种疟蚊都是疟疾的主要传播昆虫。在过去的 20 年中，作为遏制疟疾计划的一部分，蚊子种群一直受到大规模杀虫剂的影响，已经开始出现了对杀虫剂的抗性。作者希望使用 diploS/HIC 可以识别出 Ag1000G 中已验证的对抗杀虫剂有关的基因，如 Gste2（glutathione S-transferase genes，谷胱甘肽 S-转移酶基因）。

Ag1000G 数据集信息 ：

样本来源地：15
样本数量：765
基因组大小：141 MB
SNP 数量：5252w（密度很高）
基因组位点的突变概率范围（预测）：$2.8\times 10^{-9}\sim 5.5\times 10^{-9}$
有效群体大小变化区间（预测）：$10^4\sim 10^7$

特征参数（Unphased）

总计 12 个参数：$\pi$、${\theta }_w$、Tajima‘s $D$、$g_{kl}$ 方差、$g_{kl}$ 偏度、$g_{kl}$ 丰度、$Z_{ns}$ (unphased)、$\omega$ (unphased)、$N_J$、$J_1$、$J_{12}$、$J_2/J_1$

1. $g_{kl}$ 表示 个体 $k$ 和 个体 $l$ 间的序列差异：$g_{kl}={\sum }_{i=0}^m{1}_{x_{ki}\ne x_{li}}$。其中 $m$ 为序列中 SNP 数量，$x$ 为基因型， x ∈ { 0 , 1 , 2 } x\in\{ 0, 1, 2 \} x∈{0,1,2}。如 $x_{ki}=2$ 为 个体 $k$ 的第 $i$ 个 SNP 的基因型为 $2$。当个体 $k$ 、$l$ 在位点 $i$ 的基因型不同时，$g_{kl}$ 的值加 $1$。
2. 作者发现 $g_{kl}$ 的方差、偏度、丰度在软、硬清扫中存在差异（下图 1），有利于提高模型对软、硬清扫的区分能力。
3. $Z_{ns}$、$\omega$ 是利用基因型数据计算出来的，用于提取 LD 信息。
4. 为了提取 SFS 信息，作者将单倍型参数 $N_H,H_1,H_{12},H_2/H_1$ 替换成为 $N_J,J_1,J_{12},J_2/J_1$。其中 $J$ 表示个体的 基因型（012 形式）， 100 个个体含有 100 个基因型。$N_J$ 表示群体中基因型种类的数量，$J_1$、$J_2$ 类似 $H_1$、$H_2$ 表示群体中数量最多、次多的基因型的占比（下图 2）。

数据输入格式

输入矩阵格式 12×11 （12行，11列，如下图），每列为一个窗口，每行为一种群体遗传学参数。图中展示了 3 种情况下各窗口内参数的平均值，可以发现不同选择强度下，图片的图案是存在差异的。作者希望通过 CNN 识别图案来确定区间的选择强度。

CNN 结构

diploS/HIC 使用 D-CNN 架构（上图，Dilation CNN 称为 空洞卷积，详情参见 附录），设置了 3 个不同尺寸的卷积层（filter size），分别是 3×3、2×2（dilation 1×3）、2×2（dilation 1×4）。整体结构为（详细结构参见附录）：

$$(Conv1\to Relu\to Conv2\to Relu\to Pool\to Dropout\to Flatten)\ast 3\to Full\to Relu\to Dropout\to Full\to Relu\to Dropout\to Full\to Softmax$$

使用 Adam 算法对 D-CNN 参数进行优化，当预测精度变化差异小于 0.001 时优化停止。作者实践中发现，通常迭代次数 epoch < 20。

PS：个人认为作者在这里使用 D-CNN 只是一个噱头，因为 D-CNN 相比 CNN 的特征是捕获更大范围的信息，但在本文中没有必要。因为本文数据组成的图片大小为 12×11，图片极小，没有必要使用 D-CNN。

作者设置 3 种尺寸的卷积层目的是从 多个尺度 上进行参数联合，捕获不同参数、窗口之间的组合，挖掘出更有效的信息来提高预测精度。作者证明，这种方式的预测精度优于单一 3×3 的卷积层；还测试了不同的特征排序下模型的预测精度，结果是一致的，遗传学参数的排序差异不影响模型预测精度 。作者将 filter size 的分别更改为 12 × 3 和 12 × 2，以消除行顺序差异对模型预测的影响，结果中模型的预测精度与 fiter size=3×3、2×2 一致。

结果

1. 使用 CNN 模型的预测精度高于 RF，使用 定相基因型数据 作为输入训练的预测精度高于 非定相基因型数据（下图 1）。
2. 选择强度的 增加 会 降低 硬清扫区 与 硬清扫链接区 之间的分类精度，选择强度的 减少 会 降低 软清扫区 与 软清扫链接区 之间的分类精度。虽然选择强度的变化会降低模型多分类的准确性，但是可以比较准确的二分类，区分清扫（hard & soft）与非清扫（neutral & linked）（下图 2）。
3. 一般数据的单倍型是通过 HMM 方法根据参考单倍型推测出来的，而不是通过直接测定或亲本基因型推测。所以推测的单倍型可能是错误的，基于错误的单倍型计算出来的 $LD$、$H_1$ 等值很可能会导致 S/HIC 模型预测精度下降。作者为了验证上述观点，在基因型不变的基础上，对部分模拟的单倍型进行了修改，变成错误的单倍型并训练。结果显示（下图 3），使用 phased 参数作为特征的模型的预测精度随单倍型错误率的提升而下降，当错误率达到 0.05 时，unphased 模型预测精度超过 phased 。
4. 实际案例是 A. gambiae 基因组的 Gste 基因（glutathione S-transferase gene，谷胱甘肽 S-转移酶基因），已被证明与拟除虫菊酯和 DDT 的解毒有关，帮助 A. gambiae 抵抗杀虫剂（下图 4）。作者设置了不同物理窗口大小（10KB、20KB、50KB）的模型，而不再是先前训练中的 110 KB。结果发现，diploS/HIC 可以有效缩小选择区间的范围。

附录

空洞卷积（Dilated Convolution，D-CNN）

D-CNN 是一种专门用于语义分割的 CNN 算法。由于语义任务往往需要涉及 全局信息等大范围信息，所以要想准确判断语义，需要在卷积中 捕获较大范围的信息。CNN 中捕获较大信息的方法可以是 扩大 kernel 尺寸，但这样会 大幅增加学习负担；也可以是 增加层数，层与层之间 使用 pool 层对关键信息进行提取，但这样会造成 信息的丢失。为此，Fisher Yu 等人提出了 D-CNN 算法，以实现对大范围信息的提取。

D-CNN 相比 CNN 在卷积层矩阵中注入了权重为 0 的点。如下图中，左图 为 3×3 kernel 的 CNN，右图 为 5×5 kernel 且 dilation rate=2 的 D-CNN（dilation rate=1 即为 CNN）。当 dilation rate=2 时，5×5 的 kernel 中虽然有 25 个权重，但其中只有 9 个是非 0 的，剩下 16 个权重都是 0。相比于 3×3 kernel 的 CNN，D-CNN 在没有增加权重数量的情况下扩大了 kernel 的捕获范围 。或者说，D-CNN 在不减少信息的情况下降低了输出维度。如 3×3 kernel，stride=1，padding=1 的 CNN 卷积 100×100 的图片时输出尺寸为 100×100，如果使用 3×3 kernel，dilation rate=5，stride=1，padding=1 的 D-CNN 卷积后输出尺寸为 92×92。避免使用 Pool 层降低矩阵大小而造成信息丢失。当然，增加 kernel 大小也可以降低输出维度，但会增加学习负担，如相比 3×3，11×11 的参数扩增了13.44 倍，大幅提高了学习成本。

空洞卷积优势：

1. 增加卷积过程所捕获的特征范围（感受野）；
2. 避免 Pool 层减少数据维度的同时造成信息的丢失；

PS：个人认为作者在这里使用 D-CNN 只是一个噱头，因为 D-CNN 相比 CNN 的特征是捕获更大范围的信息，但在本文中没有必要。因为本文数据组成的图片大小为 12×11，图片极小，没有必要使用 D-CNN。从作者每次卷积输入与输出格式相同（12×11）即可看出。

参考来源：知乎：如何理解空洞卷积（dilated convolution）？（刘诗昆回答）

CNN 详细结构（以 Filter 1 为例）：

1. 12×11 的输入图片经过 128 个 3×3 filter size 的卷积层 1（conv1），生成 128 个 12×11（padding=”same“，strides=1，”same“ 表示通过 padding 填补使输出结构与输入相同）的图片并计算均值，得到 1 个 12×11 的图片，Relu 激活；
2. 经过 64 个 3×3 filter size 的卷积层 2（conv2），生成 64 个 12×11（padding=”same“，strides=1）的图片并计算均值，得到 1 个 12×11 的图片，Relu 激活；
3. 经过 Pool 层（pool_size=3，strides=(3×3)，padding=”same“），生成 1 个 4×4 的图片；
4. 经过 Dropout 层（rate=0.15，仅在训练时发挥作用，预测时不会 dropout），将 4×4 中随机 15% 的值替换成 0。
5. 将 4×4 的矩阵展开成 1×16 的一维数据（Flatten）；
6. 同理，输入图片经过两次 2×2（dilation 1×3）filter size 的卷积并展开后，生成 1×16 的一维数据；
7. 同理，输入图片经过两次 2×2（dilation 1×4）filter size 的卷积并展开后，生成 1×16 的一维数据；
8. 将 3 个 filter 的输出合并（h=concatenate([filter1, filter2, filter3])），生成 1×48 的一维数据作为全连接层 1（1×512）的输入；
9. 经过全连接层 1（1×512）后，生成 1×512 的一维数据，Relu 激活；
10. 经过 Dropout 层（rate=0.2）；
11. 经过全连接层 2（1×128）后，生成 1×128 的一维数据，Relu 激活；
12. 经过 Dropout 层（rate=0.1）；
13. 经过全连接层 3（1×5）后，生成 1×5 的一维数据，Softmax 激活并输出。