S/HIC 系列软件：S/HIC 利用随机森林识别软/硬清扫

内容翻译整理自文章 S/HIC: Robust Identification of Soft and Hard Sweeps Using Machine Learning. Schrider DR, Kern AD. PLoS Genet. 2016;12(3):e1005928. Published 2016 Mar 15. doi:10.1371/journal.pgen.1005928，被引用次数 118。

注意：本文中的 软、硬清扫 不是由 选择强度 定义，而是由受选择突变的 起始频率 定义。当选择位点为 从头突变 时，产生的图谱称为 硬清扫，当选择位点为 常态突变 时，产生的图谱称为 软清扫。

S/HIC 软件特色

相比已有软件，S/HIC 软/硬 清扫图谱的识别精度 得到了提高。
S/HIC 软/硬 清扫图谱中受选择位点的识别更为精准，相比于其他软件，S/HIC 设立了新的分类：软清扫 link 区、硬清扫 link 区。通过 link 区来降低受选择位点的可能区间，尽可能排除软肩效应（soft shoulder），减少不必要的分析（图）。同时，也可以准确的识别中性区间。
通过 大区间信息的纳入 和遗传学参数 归一化 两种方法，使 S/HIC 模型对经历了不同历史事件的群体的图谱识别具有较高的 鲁棒性。由于掌握信息有限，群体是否经历了瓶颈、快速扩张、迁移等事件主要是通过推测。所以难免会出现推测错误的情况，S/HIC 使用错误的模型进行图谱识别时，虽然精度会明显下降，但相比其他软件，仍然是精度最高的。
特征、训练集都可以不断拓展。如果未来有新的群体遗传统计参数，也可以方便的纳入模型中。如果未来有新的演化理论，也可以通过添加对应的模拟数据集来提高模型的全面性。

S/HIC 模型结构

各种群体遗传统计参数包含：πππ、θwθ_wθw、θHθ_HθH、不同单倍型的数量、平均单倍型纯合度、H12H_{12}H12、 H2/H1H_2/H_1H2/H1（单倍型频率 top1 的频率 / 单倍型频率 top2 的频率）、ZnS(LD)Z_{nS}(LD)ZnS(LD)、ωωω 最大值。
群体数量大幅增长（founder effect，S 型曲线）、群体数量大幅降低（bottleneck effect）等事件会造成遗传统计参数大幅改变。作者使用 归一化 （π1∑πi\dfrac{π_1}{\sumπ_i}∑πiπ1）后的统计参数作为模型的输入。使窗口内统计参数反映的不是绝对值，而是趋势、曲线形状。通过比较趋势而非绝对值，来判断窗口是否受到选择，尽量避免群体数量变化对模型预测的影响。
通过窗口内曲线的形状判断中心子窗口的类别，是利用 整个窗口更丰富的信息，来避免因随机因素导致位点上参数值的偏差而降低分类预测的精度，提高模型的鲁棒性。
以滑窗的方式读取 SNP 序列的全部信息，窗口（2200 kb）均等的分割为 11 个小的子窗口（200 kb），每次滑动步长是一个子窗口。每个子窗口（200kb）会计算出上述的 9 个群体遗传统计参数，总共得到 99 个参数值，并根据大窗口（2200kb）的和进行归一化。窗口（2200 kb）的输入特征向量形式如下：
[π1∑πi,π2∑πi,…,π11∑πi,θw1∑θwi,θw2∑θwi,…,θw11∑θwi,θH1∑θHi,…,ω1∑ωi,…,ω11∑ωi],∑πi=∑i=111πi[\dfrac{π_1}{\sumπ_i}, \dfrac{π_2}{\sumπ_i}, …, \dfrac{π_{11}}{\sumπ_i}, \dfrac{θ_{w1}}{\sumθ_{wi}}, \dfrac{θ_{w2}}{\sumθ_{wi}}, …, \dfrac{θ_{w11}}{\sumθ_{wi}}, \dfrac{θ_{H1}}{\sumθ_{Hi}}, …, \dfrac{ω_1}{\sumω_i}, …, \dfrac{ω_{11}}{\sumω_i}], {\sumπ_i} = {\sum_{i=1}^{11}π_i} [∑πiπ1,∑πiπ2,…,∑πiπ11,∑θwiθw1,∑θwiθw2,…,∑θwiθw11,∑θHiθH1,…,∑ωiω1,…,∑ωiω11],∑πi=i=1∑11πi
S/HIC 使用 RF（随机森林） 模型，执行 5 分类 问题，标签为软、硬、中、软侧（soft-linked）、硬侧（hard-linked）。样本的标签取决于窗口的 中心子窗口（6） 的类别。当中心子窗口发生了软、硬清扫时，窗口标签为软、硬清扫；当中心子窗口为中性但链接区发生软、硬清扫时，窗口标签为软、硬清扫链接区；当窗口内所有子窗口均为中性时，窗口标签为中性。

下图展示了 S/HIC 的预测示意图，通过识别 11 个子窗口中所蕴含的 πππ 的信息，来预测中心子窗口的类型。不同类型用不同颜色表示：硬蓝、软橙、中灰、硬测紫、软侧黄。

平衡群体模拟

平衡群体模型（群体大小不变）

使用 discoal 软件模拟群体数据，模拟参数包含：

受选择位点的选择系数（ααα，α=2N0sα = 2N_0sα=2N0s）：强 ~ U(2.5×103,2.5×104)U(2.5×10^3, 2.5×10^4)U(2.5×103,2.5×104)、中 ~ U(2.5×102,2.5×103)U(2.5×10^2, 2.5×10^3)U(2.5×102,2.5×103)、弱 ~ U(25,2.5×102)U(25, 2.5×10^2)U(25,2.5×102)
初始基因型频率（fff）：中频 ~ U(0.05,0.2)U(0.05, 0.2)U(0.05,0.2)、低频 ~ U(2/2N0,0.05)U(2/2N_0, 0.05)U(2/2N0,0.05)、中低频 ~ U(2/2N0,0.2)U(2/2N_0, 0.2)U(2/2N0,0.2)
染色体长度：L=110KBL=110 KBL=110KB
群体突变率：4N0μL=U(528,1584)4N_0μL=U(528, 1584)4N0μL=U(528,1584)
染色体重组率：4N0rL=8804N_0rL=8804N0rL=880
突变固定距观测的时间：t=U(0,0.1)×2N0t=U(0, 0.1)×2N_0t=U(0,0.1)×2N0
模拟完成后从模拟群体中随机抽取 100 个染色体作为代表，计算 πππ、θwθ_wθw、θHθ_HθH 等参数

PS：作者在设定模拟数据集时，没有将参数固定在某一具体数值上，而是从 均匀分布 UUU 中随机抽样得到，目的是为了提高模型的鲁棒性。

对于下面多种情况的测试，训练集（5000）均有 5 类大小相同（1000）的模拟集构成：

1000 个 硬清扫 模拟窗口，即 从头突变 的选择固定发生在中央子窗口。
1000 个 软清扫 模拟窗口，即 常态突变（ f=U(0.05,0.2)f=U(0.05, 0.2)f=U(0.05,0.2)）的选择固定发生在中央子窗口。
1000 个 硬清扫链接 模拟窗口，即 从头突变 的选择固定发生在 非中央 子窗口。discoal 对每个非中央子窗口的固定都进行了 1000 次硬清扫模拟，总共生成 10000 个硬清扫链接模拟窗口。作者从 10000 个模拟集中随机抽取 1000 个模拟作为硬清扫链接训练集。
1000 个 软清扫链接 模拟窗口，即 常态突变（ f=U(0.05,0.2)f=U(0.05, 0.2)f=U(0.05,0.2)）的选择固定发生在 非中央 子窗口。discoal 用相同的策略生成 10000 个软清扫链接模拟窗口。作者从中随机抽取 1000 个模拟作为软清扫链接训练集。
1000 个 无清扫 的中性模拟窗口。

PS：对于不包含软清扫链接和硬清扫链接分类的模型，链接区与无清扫中性区都统计为中性区。

比较软件：

SweepFinder（composite-likelihood ratio，CLR，复合似然比）；
SFselect（SVM）： SFS；
evolBoosting（boosting）： Tajima’s D，Fay and Wu’s H，integrated haplotype homozygosity（iHH）；
Approximate Bayesian Computation-like approach（Bayesian）：Bayes Factors；

判断指标：ROC 曲线（receiver operating characteristic curve）

测试 1 ：二分类训测同选择强度区分硬清扫与中性（中性）

模拟强、中、弱选择强度的 3 组训练集，每组训练集的样本数量为 5000，训练得到 3 个模型。测试集也为强、中、弱选择强度的 3 组，其中每组包含 2000 个样本，1000 个硬清扫模拟窗口与 1000 个无清扫的中性模拟窗口。用在各自选择强度训练集训练下的模型预测各自选择强度的测试集。结果显示，模型可以准确的区分硬清扫区与中性区。

测试 2 ：二分类训测异选择强度区分硬清扫与中性（中性）

因为实际染色体上区间的选择强度是未知的，所以匹配型预测没有实际价值，需要使用统一训练好的模型来区分各个选择强度下的清扫区与中性区。作者将测试 1 中的强、中、弱选择强度训练集组合，构成样本数量为 15000 的大训练集，训练得到 1 个模型。测试集与测试 1 相同，分 3 组，用同 1 个模型分别在 3 种选择强度下的测试集验证。结果显示，模型异选择强度下可以准确的区分硬清扫区与中性区。

测试 3 ：二分类训测异选择强度区分清扫（软 + 硬）与中性（中性）

训练集与测试 2 相同。测试集结构与测试 1 一致，分 3 组，每组为 1000 个清扫模拟 + 1000 个无清扫的中性模拟。但 1000 个清扫模拟变为 500 个 f=0f = 0f=0 的从头突变和 500 个 f=U(0.05,0.2)f = U(0.05, 0.2)f=U(0.05,0.2) 的常态突变。测试集中样本中心子窗口无论被预测为硬清扫区还是软清扫区，都统计为清扫区。结果显示，模型可以准确的区分清扫区（软+硬）与中性区。

测试 4 ：二分类训测异选择强度区分软、硬清扫

训练集与测试 2 相同。测试集 3 组，每组为 1000 个硬清扫模拟 + 1000 个软清扫 f=U(0.05,0.2)f = U(0.05, 0.2)f=U(0.05,0.2) 模拟。结果显示，模型可以在清扫区中准确的区分从头突变与常态突变（即区分软、硬清扫区）。

测试 5 ：二分类训测异选择强度区分软、硬清扫

训练集、测试集与测试 4 相同。但训练集与测试集中的软清扫模拟参数由 f=U(0.05,0.2)f = U(0.05, 0.2)f=U(0.05,0.2) 修改为 f=U(2/2N0,0.05)f = U(2/2N_0, 0.05)f=U(2/2N0,0.05) 。结果显示，模型无法准确的区分从头突变与低频常态突变。

测试 6 ：二分类训测异选择强度区分清扫（软 + 硬）与中性区（中性 + 链接）

训练集与测试 2 相同。测试集 3 组，每组为 1000 个清扫模拟 + 1000 个中性模拟。其中 1000 个清扫模拟包含 500 个硬清扫模拟和 500 个 f=U(0.05,0.2)f = U(0.05, 0.2)f=U(0.05,0.2) 的软清扫模拟；1000 个中性模拟包含 333 个无清扫中性模拟， 333 个硬链接模拟和 333 次软链接模拟（硬、软链接模拟为从 10000 次非中心硬、软清扫模拟中随机抽样 333 次模拟得到）。测试集中样本中心子窗口无论被预测为硬清扫区还是软清扫区，都统计为清扫区；同样，无清扫中性区和软、硬清扫链接区都统计为中性区。结果显示，模型可以准确的区分清扫区（从头突变+常态突变）与中性区（中性+链接）。

测试 7 ：多分类训测异选择强度区分软清扫、硬清扫、中性（中性 + 链接）

训练集与测试 2 训练集相同。测试集分 3 组，每组与测试 2 训练集相同，但重新生成。对测试集的预测由测试 1-5 的二分类改为了多分类。结果显示，模型可以准确的区分软清扫区、硬清扫区与中性区（中性+链接）。

非平衡群体模拟

非洲人口模型（不断收缩）

模拟参数：

受选择位点的选择系数：强 ~ U(5.0×104,5.0×105)U(5.0×10^4, 5.0×10^5)U(5.0×104,5.0×105)、中 ~ U(5.0×103,5.0×104)U(5.0×10^3, 5.0×10^4)U(5.0×103,5.0×104)
初始基因型频率：f=U(2/2N0,0.2)f=U(2/2N_0, 0.2)f=U(2/2N0,0.2)
染色体长度：L=110KBL=110 KBL=110KB
群体突变率：4N0μL=U(10560,31680)4N_0μL=U(10560, 31680)4N0μL=U(10560,31680)
染色体重组率：4N0rL=176004N_0rL=176004N0rL=17600
群体不断缩小，最终在 0.00024×4N00.00024×4N_00.00024×4N0 代缩减至 0.034N00.034 N_00.034N0
突变固定距观测的时间：t=U(0,2.44×10−4)×2N0t=U(0, 2.44×10^{-4})×2N_0t=U(0,2.44×10−4)×2N0
模拟完成后从模拟群体中随机抽取 100 个染色体作为代表，计算 πππ、θwθ_wθw、θHθ_HθH 等参数

欧洲人口模型（不断收缩 + 突然扩张）

模拟参数：

受选择位点的选择系数：2N0s=U(5.0×103,5.0×105)2N_0s=U(5.0×10^3, 5.0×10^5)2N0s=U(5.0×103,5.0×105)
初始基因型频率：f=U(2/2N0,0.2)f=U(2/2N_0, 0.2)f=U(2/2N0,0.2)
染色体长度：L=110KBL=110 KBL=110KB
群体突变率：4N0μL=U(5406.72,54067.2)4N_0μL=U(5406.72, 54067.2)4N0μL=U(5406.72,54067.2)
染色体重组率：4N0rL=176004N_0rL=176004N0rL=17600
群体先不断缩小，在 0.0009×4N00.0009×4N_00.0009×4N0 代缩减至 0.002N00.002 N_00.002N0，然后不断扩张，在 0.002×4N00.002×4N_00.002×4N0 代扩增至 0.028N00.028 N_00.028N0
突变固定距观测的时间：t=U(0,2.0×10−3)×2N0t=U(0, 2.0×10^{-3})×2N_0t=U(0,2.0×10−3)×2N0
模拟完成后从模拟群体中随机抽取 100 个染色体作为代表，计算 πππ、θwθ_wθw、θHθ_HθH 等参数

瓶颈效应（突然收缩 + 突然扩张）

模拟参数：

受选择位点的选择系数：2N0s=U(1.0×102,1.0×104)2N_0s=U(1.0×10^2, 1.0×10^4)2N0s=U(1.0×102,1.0×104)
初始基因型频率：f=U(2/2N0,0.2)f=U(2/2N_0, 0.2)f=U(2/2N0,0.2)
染色体长度：L=110KBL=110 KBL=110KB
群体突变率：4N0μL=11004N_0μL=11004N0μL=1100
染色体重组率：4N0rL=11004N_0rL=11004N0rL=1100
突变固定距观测的时间：t=(0,0.022,0.044)×2N0t=(0, 0.022, 0.044) × 2N_0t=(0,0.022,0.044)×2N0
模拟完成后从模拟群体中随机抽取 100 个染色体作为代表，计算 πππ、θwθ_wθw、θHθ_HθH 等参数

严重瓶颈效应模型：群体在 0.0084×4N00.0084 × 4N_00.0084×4N0 代时发生瓶颈效应，群体大小缩小为原群体大小的 2.9 % ；在 0.044×4N00.044 × 4N_00.044×4N0 代时群体大小恢复为原群体大小。

非严重 瓶颈效应模型：群体在 0.0084×4N00.0084 × 4N_00.0084×4N0 代时发生瓶颈效应，群体大小缩小为原群体大小的 29 % ；在 0.044×4N00.044 × 4N_00.044×4N0 代时群体大小恢复为原群体大小。

测试 1 ：二分类训测同非平衡区分清扫（软 + 硬）与中性区（中性 + 链接）

测试模型：非洲、欧洲、瓶颈

针对非洲、欧洲、瓶颈 3 个模型的群体数量变化情况，分别设计 5 组（瓶颈、非洲各 2 组）训练集和测试集，并训练出 5 个模型。如欧洲，根据模型参数，使用 discoal 模拟 1000 个清扫（500 软 + 500 硬）群体和 1000 个中性（333 中性 + 333 软链 + 333 硬链），总计 2000 个群体作为训练集；测试集同理。

结果显示（AUC）：

非洲 α=U(5.0×104,5.0×105)α = U(5.0×10^4, 5.0×10^5)α=U(5.0×104,5.0×105)，S/HIC = 0.9865
非洲 α=U(5.0×103,5.0×104)α = U(5.0×10^3, 5.0×10^4)α=U(5.0×103,5.0×104)，S/HIC = 0.8060
欧洲 α=U(5.0×103,5.0×105)α = U(5.0×10^3, 5.0×10^5)α=U(5.0×103,5.0×105)，S/HIC = 0.8188
严重瓶颈 t=0t=0t=0，S/HIC = 0.9182
严重瓶颈 t=0.022×2Nt=0.022×2Nt=0.022×2N，S/HIC = 0.6750
严重瓶颈 t=0.044×2Nt=0.044×2Nt=0.044×2N，S/HIC = 0.6672
非严重瓶颈 t=0t=0t=0，S/HIC = 0.9966
非严重瓶颈 t=0.022×2Nt=0.022×2Nt=0.022×2N，S/HIC = 0.9857
非严重瓶颈 t=0.044×2Nt=0.044×2Nt=0.044×2N，S/HIC = 0.9458

测试 2 ：二分类训平衡测非平衡区分清扫（软 + 硬）与中性区（中性 + 链接）

测试模型：非洲、欧洲、瓶颈

训练集使用平衡模拟中测试 2 的训练集，总计 15000 个群体；测试集与测试 1 相同，总共 5 类；使用相同的训练模型测试在不同测试集上的预测精度。

结果显示（AUC）：

非洲 α=U(5.0×104,5.0×105)α = U(5.0×10^4, 5.0×10^5)α=U(5.0×104,5.0×105)，S/HIC = 0.9844
非洲 α=U(5.0×103,5.0×104)α = U(5.0×10^3, 5.0×10^4)α=U(5.0×103,5.0×104)，S/HIC = 0.8122
欧洲 α=U(5.0×103,5.0×105)α = U(5.0×10^3, 5.0×10^5)α=U(5.0×103,5.0×105)，S/HIC = 0.8127
严重瓶颈 t=0t=0t=0，S/HIC = 0.7817
严重瓶颈 t=0.022×2Nt=0.022×2Nt=0.022×2N，S/HIC = 0.6790
严重瓶颈 t=0.044×2Nt=0.044×2Nt=0.044×2N，S/HIC = 0.6497
非严重瓶颈 t=0t=0t=0，S/HIC = 0.9927
非严重瓶颈 t=0.022×2Nt=0.022×2Nt=0.022×2N，S/HIC = 0.9841
非严重瓶颈 t=0.044×2Nt=0.044×2Nt=0.044×2N，S/HIC = 0.9491

测试 3 ：多分类训平衡测非平衡区分软清扫、硬清扫、中性（中性 + 链接）

测试模型：非洲、欧洲

训练集与测试 2 训练集相同。测试集与测试 2 训练集相同，但重新生成。对测试集的预测由测试 1-5 的二分类改为了多分类。下图为非洲模型的预测结果。

真实数据测试

数据使用的是 1000 Genomes Project’s CEU 群体的 18 号染色体数据。在使用数据前对数据进行了过滤，去除了染色体上的重复片段（ReapeatMasker）与测序质量较低的区间（ 50 bp read mappability scores < 1）。过滤后，18 号染色体总共被分为 344 个 200kb 的小窗口，其中 9.9% 的窗口被识别为硬清扫区，6.4% 为硬清扫链接区，14% 为软清扫区，25.9% 为软清扫链接区，43.9% 为无清扫中性区。

最后的案例分析上，S/HIC 选择了 L3MBTL4 基因所处的区间，基因已经被研究过了，但该区间先前没有被其他软件识别出受到了选择。下图中红色表示硬清扫，粉红表示硬清扫链接，淡蓝表示软清扫链接，黑色表示中性。

实验结果讨论

在非平衡实验中的瓶颈模拟中可以发现，H12 对时间敏感，随着时间的增加，单倍型会被迅速破坏，H12 等特征适合挖掘最近产生的 sweep。
从各种模拟的实验结果来看，S/HIC 的预测能力相比于已有的 sweep 挖掘软件来说，无疑是更好的。但预测精度提高的原因是否与使用 RF 构建模型有关，还需要进一步讨论。因为在对比中，其他软件使用非归一化的参数值作为特征，并且在训练中没有 link 类标签。
附件 S2 Table 中显示了各个输入特征对不同 S/HIC 模型的重要性，值得注意的是， S2 Table 中所列的特征与文中给出的特征存在出入。表中特征为 πππ、θHθ_HθH、Fay and Wu’s H、H1H_1H1、H12H_{12}H12、 H2/H1H_2/H_1H2/H1、ss（不知道代表什么指标，估计为 ZnSZ_{nS}ZnS ）、Tajima’s D。以平衡模拟训练集训练出的模型为例，可以发现，重要性 top5 的特征都是中性窗口（win5）中遗传统计参数，top10 中有 6 个是 πππ、4 个边缘窗口（win 0-1、9-10），top15 中有 9 个是 πππ、3 个 ss、6 个边缘窗口，top25 中有 11 个是 πππ、10 个 ss、边缘窗口 9 个。综上，中心窗口和边缘窗口对模型预测贡献较大， πππ 和 ss 遗传参数对模型贡献较大。作者没有检测其他参数删除后模型精度的变化情况，所以这里对特征是否存在冗余存疑。

S/HIC 方法应用于作物基因组的局限性

挖掘受选择区间的窗口过大，分辨率为 200 kb，适用于基因片段较大的物种，如人类（基因平均长度为 72 kb，L3MBTL4 基因长度为 460 kb），不适用于基因片段小的物种，如水稻（基因平均长度为 4.5 kb）、玉米等。所以若要确定水稻中某基因是否受到选择，则需要 大幅缩小窗口大小，如 10 kb。
S/HIC 使用相同的模型处理基因区和非基因区。基因中的外显子区因为表达蛋白质，相较于非表达区更为保守，即表达区与非表达区之间的群体遗传学参数存在差异。但在人类基因组中这种差异不需要着重考虑，因为人类基因中的表达区极少，被庞大的非表达区稀释了。如人类 L3MBTL4 基因长度为 460 kb，蛋白质长度 623 Aa，相当于 1.9 kb（0.4%），L3MBTL4 基因中有约 458.1 kb 的非表达区。但玉米基因中非表达片段较少，如 ZCN8 基因（pebp8）长度为 4.8 kb，蛋白质长度 175 Aa，相当于 525 bp（10.9%）。所以，当为了定位植物基因组上受选择基因而缩小窗口后，如果窗口扫描到基因区，会因为窗口内群体遗传学参数的改变，而误导模型将基因区识别为受选择区。所以，对植物基因组受选择区挖掘时，需要考虑构建 基因区与非基因区两种模型 。
S/HIC 因为窗口较大（200 kb）所以不用担心窗口内无 SNP 的问题。但当扫描植物基因组而缩小窗口大小后，需要考虑窗口内无 SNP 的情况，这会引出新问题：窗口内子窗口群体遗传参数缺失 。
S/HIC 在预测非平衡群体时，预测精度下降。自然群体中，群体数量大规模变化的次数不多，而作物等受人工选择的群体会经常性的出现数量大规模变化的情况，S/HIC 的预测精度可能会发生显著下降。所以，针对自然群体设计的 sweep 挖掘软件可能不适合人工群体。