本文内容摘取自 论文: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 算法:
- diploS/HIC 使用 卷积神经网络(CNN) 作为模型框架,替换了 S/HIC 中的 随机森林(RF) 框架。
- diploS/HIC 从 多个尺度 上捕获不同参数、窗口之间的组合,挖掘出更有效的信息来提高预测精度。
- S/HIC 的输入特征包含 定相基因型数据(phased) 计算出的群体遗传统计参数,如 单倍型 等参数,但 diploS/HIC 的输入特征仅使用 非定相基因型数据(unphased) 计算出的群体遗传统计参数,扩大了模型的应用范围。
- 由参考单倍型推测出的单倍型数据,可能会由于错误率较高(error rate>0.05)而使定相数据训练出的模型精度较低,作者推荐使用非定相数据作为特征训练模型。
PS:(非)定相基因型数据参见:基因定相(Phasing) 与 SHAPEIT 原理简介
数据集
模拟数据集(Coalescent simulation)
窗口数量:
11
11
11
模拟次数(训练集):
2000
2000
2000
模拟次数(测试集):
1000
1000
1000
平衡群体
窗口大小: L = 110 K B L=110KB L=110KB
突变率:
4
N
0
μ
L
=
U
(
528
,
1584
)
4N_0μL=U(528, 1584)
4N0μL=U(528,1584)
重组率:
4
N
0
r
L
=
880
4N_0rL=880
4N0rL=880
选择强度:
2
N
0
s
=
U
(
25
,
2.5
×
1
0
2
)
2N_0s=U(25, 2.5×10^2)
2N0s=U(25,2.5×102)、
U
(
2.5
×
1
0
2
,
2.5
×
1
0
3
)
U(2.5×10^2, 2.5×10^3)
U(2.5×102,2.5×103)、
U
(
2.5
×
1
0
3
,
2.5
×
1
0
4
)
U(2.5×10^3, 2.5×10^4)
U(2.5×103,2.5×104)
软清扫起始频率:
f
=
U
(
0
,
0.2
)
f=U(0, 0.2)
f=U(0,0.2)
固定与观察间的时间差:
τ
=
U
(
0
,
0.05
)
τ=U(0, 0.05)
τ=U(0,0.05)
群体大小:保持不变
非平衡群体(模拟 Ag1000G)
窗口大小: L = 55 K B L=55KB L=55KB
突变率:
4
N
0
μ
L
=
U
(
1750
,
17500
)
4N_0μL=U(1750, 17500)
4N0μL=U(1750,17500)
重组率:
4
N
0
r
L
=
U
(
1925
,
57756
)
4N_0rL=U(1925, 57756)
4N0rL=U(1925,57756)
选择强度:
2
N
0
s
=
U
(
2.5
×
1
0
2
,
2.5
×
1
0
3
)
2N_0s=U(2.5×10^2, 2.5×10^3)
2N0s=U(2.5×102,2.5×103)、
U
(
2.5
×
1
0
3
,
2.5
×
1
0
4
)
U(2.5×10^3, 2.5×10^4)
U(2.5×103,2.5×104)、
U
(
2.5
×
1
0
4
,
2.5
×
1
0
5
)
U(2.5×10^4, 2.5×10^5)
U(2.5×104,2.5×105)
软清扫起始频率:
f
=
U
(
0
,
0.2
)
f=U(0, 0.2)
f=U(0,0.2)
固定与观察间的时间差:
τ
=
U
(
0
,
0.00004
)
τ=U(0, 0.00004)
τ=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
×
1
0
−
9
∼
5.5
×
1
0
−
9
2.8×10^{-9} \sim 5.5×10^{-9}
2.8×10−9∼5.5×10−9
有效群体大小变化区间(预测):
1
0
4
∼
1
0
7
10^4 \sim 10^7
104∼107
特征参数(Unphased)
总计 12 个参数: π π π、 θ w \theta_w θw、Tajima‘s D D D、 g k l g_{kl} gkl 方差、 g k l g_{kl} gkl 偏度、 g k l g_{kl} gkl 丰度、 Z n s Z_{ns} Zns (unphased)、 ω ω ω (unphased)、 N J N_J NJ、 J 1 J_1 J1、 J 12 J_{12} J12、 J 2 / J 1 J_2/J_1 J2/J1
- g k l g_{kl} gkl 表示 个体 k k k 和 个体 l l l 间的序列差异: g k l = ∑ i = 0 m 1 x k i ≠ x l i g_{kl}=\sum_{i=0}^{m}1_{x_{ki}\ne x_{li}} gkl=∑i=0m1xki=xli。其中 m m m 为序列中 SNP 数量, x x x 为基因型, x ∈ { 0 , 1 , 2 } x\in\{ 0, 1, 2 \} x∈{0,1,2}。如 x k i = 2 x_{ki}=2 xki=2 为 个体 k k k 的第 i i i 个 SNP 的基因型为 2 2 2。当个体 k k k 、 l l l 在位点 i i i 的基因型不同时, g k l g_{kl} gkl 的值加 1 1 1。
- 作者发现 g k l g_{kl} gkl 的方差、偏度、丰度在软、硬清扫中存在差异(下图 1),有利于提高模型对软、硬清扫的区分能力。
- Z n s Z_{ns} Zns、 ω ω ω 是利用基因型数据计算出来的,用于提取 LD 信息。
- 为了提取 SFS 信息,作者将单倍型参数 N H , H 1 , H 12 , H 2 / H 1 N_H, H_1, H_{12}, H_2/H_1 NH,H1,H12,H2/H1 替换成为 N J , J 1 , J 12 , J 2 / J 1 N_J, J_1, J_{12}, J_2/J_1 NJ,J1,J12,J2/J1。其中 J J J 表示个体的 基因型(012 形式), 100 个个体含有 100 个基因型。 N J N_J NJ 表示群体中基因型种类的数量, J 1 J_1 J1、 J 2 J_2 J2 类似 H 1 H_1 H1、 H 2 H_2 H2 表示群体中数量最多、次多的基因型的占比(下图 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)。整体结构为(详细结构参见附录):
( C o n v 1 → R e l u → C o n v 2 → R e l u → P o o l → D r o p o u t → F l a t t e n ) ∗ 3 → F u l l → R e l u → D r o p o u t → F u l l → R e l u → D r o p o u t → F u l l → S o f t m a x (Conv1 \rightarrow Relu \rightarrow Conv2 \rightarrow Relu \rightarrow Pool \rightarrow Dropout \rightarrow Flatten)*3 \rightarrow Full \rightarrow Relu \rightarrow Dropout \rightarrow Full \rightarrow Relu \rightarrow Dropout \rightarrow Full \rightarrow Softmax (Conv1→Relu→Conv2→Relu→Pool→Dropout→Flatten)∗3→Full→Relu→Dropout→Full→Relu→Dropout→Full→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 一致。
结果
- 使用 CNN 模型的预测精度高于 RF,使用 定相基因型数据 作为输入训练的预测精度高于 非定相基因型数据(下图 1)。
- 选择强度的 增加 会 降低 硬清扫区 与 硬清扫链接区 之间的分类精度,选择强度的 减少 会 降低 软清扫区 与 软清扫链接区 之间的分类精度。虽然选择强度的变化会降低模型多分类的准确性,但是可以比较准确的二分类,区分清扫(hard & soft)与非清扫(neutral & linked)(下图 2)。
- 一般数据的单倍型是通过 HMM 方法根据参考单倍型推测出来的,而不是通过直接测定或亲本基因型推测。所以推测的单倍型可能是错误的,基于错误的单倍型计算出来的 L D LD LD、 H 1 H_1 H1 等值很可能会导致 S/HIC 模型预测精度下降。作者为了验证上述观点,在基因型不变的基础上,对部分模拟的单倍型进行了修改,变成错误的单倍型并训练。结果显示(下图 3),使用 phased 参数作为特征的模型的预测精度随单倍型错误率的提升而下降,当错误率达到 0.05 时,unphased 模型预测精度超过 phased 。
- 实际案例是 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 倍,大幅提高了学习成本。
![](https://i-blog.csdnimg.cn/blog_migrate/c1b1e7f8ab22dc45e4f4511ec0d4ba3c.png)
![](https://i-blog.csdnimg.cn/blog_migrate/328acbb829cd47003210ca0ff18dc466.png)
空洞卷积优势:
- 增加卷积过程所捕获的特征范围(感受野);
- 避免 Pool 层减少数据维度的同时造成信息的丢失;
PS:个人认为作者在这里使用 D-CNN 只是一个噱头,因为 D-CNN 相比 CNN 的特征是捕获更大范围的信息,但在本文中没有必要。因为本文数据组成的图片大小为 12×11,图片极小,没有必要使用 D-CNN。从作者每次卷积输入与输出格式相同(12×11)即可看出。
参考来源:知乎:如何理解空洞卷积(dilated convolution)?(刘诗昆回答)
CNN 详细结构(以 Filter 1 为例):
- 12×11 的输入图片经过 128 个 3×3 filter size 的卷积层 1(conv1),生成 128 个 12×11(padding=”same“,strides=1,”same“ 表示通过 padding 填补使输出结构与输入相同)的图片并计算均值,得到 1 个 12×11 的图片,Relu 激活;
- 经过 64 个 3×3 filter size 的卷积层 2(conv2),生成 64 个 12×11(padding=”same“,strides=1)的图片并计算均值,得到 1 个 12×11 的图片,Relu 激活;
- 经过 Pool 层(pool_size=3,strides=(3×3),padding=”same“),生成 1 个 4×4 的图片;
- 经过 Dropout 层(rate=0.15,仅在训练时发挥作用,预测时不会 dropout),将 4×4 中随机 15% 的值替换成 0。
- 将 4×4 的矩阵展开成 1×16 的一维数据(Flatten);
- 同理,输入图片经过两次 2×2(dilation 1×3)filter size 的卷积并展开后,生成 1×16 的一维数据;
- 同理,输入图片经过两次 2×2(dilation 1×4)filter size 的卷积并展开后,生成 1×16 的一维数据;
- 将 3 个 filter 的输出合并(h=concatenate([filter1, filter2, filter3])),生成 1×48 的一维数据作为全连接层 1(1×512)的输入;
- 经过全连接层 1(1×512)后,生成 1×512 的一维数据,Relu 激活;
- 经过 Dropout 层(rate=0.2);
- 经过全连接层 2(1×128)后,生成 1×128 的一维数据,Relu 激活;
- 经过 Dropout 层(rate=0.1);
- 经过全连接层 3(1×5)后,生成 1×5 的一维数据,Softmax 激活并输出。