S/HIC 系列软件：partialS/HIC 利用 CNN 识别不完全软/硬清扫

内容翻译整理自文章 Discovery of Ongoing Selective Sweeps within Anopheles Mosquito Populations Using Deep Learning, Molecular Biology and Evolution, March 2021。

特点

相比上一个软件 diploS/HIC ，本次作者提出的 partialS/HIC 有如下特点：

软件使用 89 个群体遗传学统计参数，使用 CNN 作为核心框架，不再使用 D-CNN。同时，遗传统计参数的排列方式遵循相似的参数放在一起，作者认为 斑驳的数据会降低卷积结果间的差异。
分类细化。从原来的 5 个扩展到 9 个：中性、硬、软、部分硬、部分软、硬侧、软侧、部分硬侧、部分软侧。
模拟细化。作者根据 Ag1000G 采样的 8 个群体进行分别的模拟，生成对应的训练集、测试集。
利用热图的数据展现形式，“稀释” 了模型的错误率。模型正确预测的精确度不高，如 Fig2 的 Completed Soft Sweep-left1 模拟下，模型正确率为 54.2%，错误率是 45.8%，但具体到每个类别上时就变成了 7.3%，0.8%，1.1%，6.6%，54.2%，0.2%，0.5%，13.9%，15.4%，可以看出预测主要还是集中在了正确的类型上。
考虑了 真实数据与模拟数据之间的部分差异，如存在 SNP 缺失、无法极化等问题。
形象的展示了各群体、不同情况下的 图谱结构，提高了说服力。

有待提升的地方：

软件使用的 89 个群体遗传统计参数过多，作者没有证明是否 参数冗余 。过多的遗传统计参数使程序准备数据的 运算时间大幅增加。
文章没有讨论部分清扫与清扫完成后被中性突变和重组而磨损的区间是否相似，即部分清扫指代的是否既包含未完成的清扫，又包含已完成许久的清扫？亦或这种图谱有更丰富的内含，如平衡选择、多基因竞争、有益基因的迁入迁出等等。
对整个染色体识别时，出现了大量受选择区间（24%），无法通过非实验的手段判断受选择区的被选择理由，即无法挖掘出与 IR 相关的基因。作者尝试通过研究选择区内基因的功能（GO 分析）、选择区与基因的位置关系来破解被选择理由，但结果并没有统一的规律。所以选择区间的挖掘，往往需要在一定先验知识下进行才有意义。
随着窗口的减小（2200KB -> 55KB），基因区保守性造成的 SNP 数量下降并没有考虑进去。这造成模型往往会将保守的基因区视为受选择区 。
使用 Masking profile 去除模拟数据中部分 SNP 的方法并不合理，会造成模拟数据中 SNP 密度偏高。
对于训练集大小为 18000 的模型来说，batch_size = 32 实在是太小了，会大幅增加模型的训练时间。是否因为 batch_size 的增加会降低模型的训练效果？

数据集

数据集模拟

因为 S/HIC 系列软件先前的工作积累，partialS/HIC 没有再模拟平衡群体下的受选择或中性区间，而是直接根据实际群体的发展历程进行模拟。同时，不再将 Ag1000G 作为 单个群体 考量，而是根据来源地和种群，将 Ag1000G 分成 8 个子群体 ：AOM (Anopheles coluzzii from Angola)、BFM (A. coluzzii from Burkina Faso)、BFS (A. gambiae from Burkina Faso)、CMS (A. gambiae from Cameroon)、GAS (A. gambiae from Gabon)、GNS (A. gambiae from Guinea)、GWA (Anopheles of uncertain species from Guinea-Bissau)、UGS (A. gambiae from Uganda)。8 子群体根据各自的群体变化史模拟中性、硬、软、部分硬、部分软、硬侧、软侧、部分硬侧、部分软侧，总计 9 种状态。

训练集 2000，每个群体的训练集 18000（2000×9），总训练集 144000（2000×9×8）。
测试集 1000，每个历史群体的测试集 45000（1000×(5+4×10)，注意测试集中作者对 10 子窗口的 linked 情况都进行了 1000 次模拟，所以 4 种 linked 总共包含 40000 次模拟），总测试集 360000（45000×8）。除了 8 个群体的模拟外，作者为群体历史背景的错误指定设计了测试集，样本数量为 495000（11=5+6，45000*11，详情参见实验设计）。

突变率（μμμ）： 3.5×10−93.5 × 10^{−9}3.5 × 10−9
模拟区间长度：555555 KB
模拟区间内子窗口数量：111111
子窗口长度：555 KB
选择强度（sss）：U(10−4,10−2)U(10^{−4}, 10^{−2})U(10−4,10−2)
固定与抽样间时间差（ttt）：U(0,2,000)U(0, 2,000)U(0,2,000)
软清扫的突变起始频率（fsoftf_{soft}fsoft）：U(1/N0,0.2)U(1/N_0, 0.2)U(1/N0,0.2)
部分清扫的突变频率（fpartf_{part}fpart）：U(0.20,0.99)U(0.20, 0.99)U(0.20,0.99)

数据预处理

Ag1000G 数据（empirical data）预处理

将 Ag1000G 基因组划分为连续的 5KB 子窗口，去除窗口内 有缺失 的 SNP、多等位 基因型的 SNP（基因型 > 2）、通过外群 无法极化（polarized）的 SNP。并将被去除的 SNP 位置信息记录在 Masking profile 中，每个文件包含了 11 个子窗口（55KB）内的被去除的 SNP 位置信息。因为 Ag1000G 基因组总共被划分为 1552 个不重合的 55KB 区间（如区间 2L：1-55000，55001-110000 等），所以总共生成 1552 个 Masking profile。
保留 剩余 SNP 数量 >= 25% 原始 SNP 数量 的子窗口。过滤会对窗口内反馈的信息产生影响，原始 SNP 越少，剩余 SNP 所表达的信息可能与原信息差异越大。本步骤 筛除了信息丢失严重，无法展现原状态的窗口。
剩下的子窗口中，满足 11 个连续 5 kb 子窗口 的形成一个完整的窗口（55KB），并用于后续的预测分析。本步骤筛除了不满足 11 个连续 5 kb 子窗口的零碎子窗口，避免了信息缺失对模型预测的影响 。

PS：极化、衍生概念参见附录

模拟数据的预处理

为了保证模拟数据与 Ag1000G 预处理后数据的结构相似，也对模拟数据的 SNP 位点进行了筛选。通过 随机抽样 Masking profile，根据文件中的 SNP 位点去除 55KB 模拟数据中的 SNP，以实现在 Ag1000G 真实数据中因 有缺失、多等位、无法极化 而被去除的 SNP 的模拟。如 discoal 生成模拟群体 PiP_iPi 后，从 1552 个文件中随机抽取 1 个剔除文件 MiM_iMi，并按照其中的位点来剔除符合的 SNP，使 PiP_iPi 中不含有 MiM_iMi 内的任何 SNP。

值得注意的是，模拟 SNP 没有使用百分比筛选，根据无限等位基因突变模型，Masking profile 中的位点大概率与模拟群体中的 SNP 不重合，所以模拟数据中 SNP 的密度会显著多于 Ag1000G。

输入特征

群体遗传统计参数：π\piπ、Tajima’s DDD、θHθ_HθH、Fay-Wu’s HHH、单倍型种类、H1H_1H1、H12H_{12}H12、H2/H1H_2/H_1H2/H1、ZnSZ_{nS}ZnS、最大 ωωω、E(iHS)E(iHS)E(iHS)、最大 iHSiHSiHS、异常 iHSiHSiHS 比例、“双倍型” 方差、“双倍型” 偏度、“双倍型” 峰度、72（6×12）个 SAFE 值。其中 “双倍型” 概念参见 S/HIC 系列软件：diploS/HIC 利用 CNN 和非定向基因型数据识别软/硬清扫

上图中，每个热图（89×11）中 row 为遗传统计参数，column 为子窗口序号，value 为群体均值。遗传统计参数的排列方式遵循相似的参数放在一起，作者认为斑驳的数据会降低卷积结果间的差异。热图反映了 AOM、BFM 两个群体不同选择情况下窗口的图谱结构。可以发现：

中性区间（Neutral Region）无规律，清扫区间有 各自特征 且主要集中在中心子窗口，为模型预测的 准确性 提供了理论基础。
除了 Completed Soft Sweep 外，其余情况下不同群体间 图谱差异不大 ，为模型预测的 鲁棒性 提供了理论基础（训练集与测试集间 群体变化历史不同）。

CNN 模型参数

epoch = 20，batch_size = 32，optimizer = Adam，CV = 10-fold

CNN 结构如下图所示（图片与代码存在出入），代码输入数据格式为 89×11，经过 Conv 1（3×6 kernel，padding=“same”，filters=256）、ReLU 后格式为 89×11，经过 Pool 1（3×3 kernel，padding=“same”）后格式为 30×4；经过 Conv 2（3×3 kernel，padding=“same”，filters=256）、ReLU 后格式为 30×4，经过 Pool 2（3×3 kernel，padding=“same”）后格式为 10×2；经过 Dropout 层（0.25）后数据展开，格式为 1×20；经过全连接层（512）、ReLU、Dropout 层（0.5）；经过全连接层（128）、ReLU、Dropout 层（0.5）；经过全连接层（9，进行 X 分类任务时更改为 X）、Softmax 输出。

实验设计

首先，作者设计了 3 种训练方案：

在每个群体中，以 9 种状态的训练样本（18000）作为训练集，构建 8 个模型。观察其在各自测试集（45000）的预测精度，评估 partialS/HIC 模型在 正确群体历史背景 下的预测精度。

仅利用方案 1 中中性、硬、硬链接、软、软链接 5 个状态的训练样本（10000）作为训练集，测试集与方案 1 相同（45000）。观察模型在不进行部分硬、部分软训练的情况下，将部分硬、部分软分类到硬、软类中的能力。评估模型 设定新分类（部分硬、部分软）的必要性。结果显示，新分类是有必要的，可以提高模型的预测能力。
训练集、测试集同方案 1，但训练目的改为判断窗口是否受到选择，partialS/HIC 为 2 分类模型（最后全连接层 nodes 由 9 更改为 2）。评估模型在区分软清扫、部分软清扫能力较差时，确保准确 区分选择与非选择区间 的能力。同时，实验还检测了部分遗传参数单独作为输入特征时预测精度，评估利用 CNN 将多个参数整合应用的必要性 。

其次，作者设计了 3 种训练集和测试集间群体历史发展背景不匹配的训练方案，评估模型在背景指定错误时的预测精度：

将在 GAS 群体（群体大小变化较为稳健）模拟训练集上训练的 CNN 应用于 CMS 群体（群体大小经历了快速扩增）的模拟测试集。
模拟了经历 2 轮瓶颈（瞬时收缩）的 5 组群体（45000×5），瓶颈的严重程度分别为 20、40、60、80、100 倍。其中 N 倍表示群体大小缩小了 N 倍。训练集为按照 BFS 模拟出的训练集（18000），测试集为 5 组瓶颈中的 1 组（45000）。
模拟了 6 组 θ\thetaθ 不同的群体（45000×6），θ=4Neμ=1000,5,000,10000,15000,20000,25000\theta = 4N_e\mu = 1000, 5,000, 10000, 15000, 20000, 25000θ=4Neμ=1000,5,000,10000,15000,20000,25000，其中 μ\muμ 为位点突变率。训练集为按照 BFS 模拟出的训练集（18000），测试集为 6 组中的 1 组（45000）。

Ag1000G 预测结果

8 个群体的 partialS/HIC 预测结果统计：

类型	硬清扫	硬清扫链接	部分硬清扫	部分硬清扫链接	软清扫	软清扫链接	部分软清扫	部分软清扫链接	中性区间
比例（中位数）%	0.03	0.14	2.84	7.55	5.01	15.84	7.24	22.82	29.67

从表中可以看出，除了中性（29.67）、链接（0.14+7.55+15.84+22.82=46.35）外剩余的约 24% 区间都受到了选择。

附录

衍生等位基因型（derived allale）与祖先等位基因型（ancestral allele）

突变会产生新的等位基因型，当新等位基因型在群体中的频率达到 5% 以上时，新的等位基因型被称为 衍生等位基因型，原始的非突变的基因型被称为 祖先基因型。如位点 MMM 在群体中只有 AAA 一种基因型，新突变产生了 GGG 基因型，当 GGG 在群体中频率达到 5% 以上时，GGG 被称为衍生等位基因型，AAA 被称为祖先等位基因型。

极化（polarized）数据

极化（polarization）是指通过外群（outgroup）确定位点中祖先和衍生等位基因型的过程。极化数据（polarized data）是指每个位点等位基因的祖先和衍生状态是已知的数据。

极化通常是使用外群（outgroup）来实现的：如果 同源位点 在外群中是单态的并且与所研究群体中的等位基因之一重合，则该变体称为祖先。如果外群中位点是 **多态 **的，或与所研究群体中等位基因都 不重合，则无法确定位点的祖先和衍生等位基因型，即此位点为无法极化的位点。

极化结果的准确与否与外群物种的选择密切相关。如果外群与研究群体在系统发育上太远，则遗传背景差异较大，导致祖先等位基因型识别错误而产生大量（伪）正选择区间。如进化关系为物种 X（外群）-> 物种 Y（祖先）-> 物种 Z（研究群体），位点 M 在外群中为基因型 AAA，在祖先和研究群体中都是 C、AC、AC、A，其中基因型 CCC 是主频且中性。将基因型 AAA 视为祖先基因型，会导致大量遗传背景的中性突变被识别为受正选择的突变。如果外群与研究群体太近，分化时间较短，大量位点在外群与研究群体间多态性是一致的，研究群体中只有少量多态位点满足外群单态的情况，大部分多态位点在外群中也是多态，造成大量的突变无法确定极性。

另一方面，全基因组范围的极化需要有研究群体的 参考基因组，将外群的基因组比对到研究群体的参考基因组上。如果仅极化部分片段，则无需参考基因组，寻找外群对应的同源片段即可。外群和研究群体的基因组间的非同源片段无法被极化。例如，千人基因组项目（1000 Genomes Project）中有约 4% 的 SNP 无法极化。