前言

王向峰老师2020年发表在《Genomics Proteomics Bioinformatics》（IF=6.597）上的文章。对于做数据分析的人来说，如何挖掘公共数据，如何从海量SNP中挖掘目标SNP等问题都是每天要面对的，这篇文章给了一个参考，很值得学习。

image.png

文章从水稻变异数据库RVD的子库IC4R（http://variation.ic4r.org/）中5152份水稻材料的18m SNPs进行层层过滤，以不同条件筛选到不同大小的SNP集，每个数据集可应用于不同场景和目的。最后还构建了数据库，提供web工具、代码及数据下载。

四个SNP集

image.png

用途：

• hapmapSNPs可用于GWAS
• tagSNPs可用于群体遗传研究和GS
• fixedSNPs可用于种子纯度和遗传背景分析
• barcodeSNPs可用于指纹图谱进行品种分类

各SNP集处理标准如下：

hapmapSNPs

• 过滤基因型缺失率大于20%的样本，剩余2556个样本。
• 过滤缺失率大于0.1，次等位基因频率MAF小于0.05的SNP。
• Beagle填充2556个样本的基因型。
• 最终包含2,097,405个SNPs，无任何缺失值。

tagSNPs

• 采用基于LD的SNP修剪步骤从hapmapSNPs来推断单倍型标记SNP（tagSNPs）。
• 水稻的LD长度40-500kb。
• Plink --indep命令，参数基于方差膨胀因子（VIF），使用滑动窗口50个SNPs，步长5个SNPs来连续过滤SNP。

fixedSNPs

• 通过比较栽培稻亚群和野生稻的Fst和θπ，鉴定选择清除区域。
• 使用100kb和10kb窗口分别鉴定大和小的基因组选择清除信号区域，6个亚群被选择区域的Tajima'D显著小于其他区域，具体为：227 (cultivated vs. wild), 381 (Ind vs. wild), 333 (Aus vs. wild), 296 (Aro vs. wild), 256 (TrJ vs. wild) and 269 (TeJ vs. wild)。
• 鉴定受选择清除区域的基因，这些基因内共1180 SNPs。

barcodeSNPs

• 使用MinimalMarker算法来详尽遍历所有可能的基因型组合，来区分这2556份材料。
• MinimalMarker算法生成最少标记组合的三个集合，每个集合包含28个SNPs。
• 合并三个集合后，共得到38个barcodeSNPs。

hapmapSNPs的指标统计

• 每步处理的指标统计
• ARNOVAR注释

image.png

tagSNPs的群体结构验证

高密度SNP对于GWAS的功能位点鉴定是有用，但对于群体遗传分析是不合适的，因为SNP的高冗余会带来不必要的计算成本，也会对结果带来偏差。
在同一LD block区，一个有代表性的SNP（tagSNP）可解决冗余问题。

156,502个tagSNPs来验证2556份材料的亚群分类和起源验证。
K=3时，能明显区分籼粳和Aus三类亚群；当K=8时，能清晰划分6个亚群；当K=4-7时，籼稻亚群能划分6个子群S1-S6。

image.png

tagSNPs的遗传多样性

5个亚群多项遗传多样性指标统计比较分析：

• Identity by state (IBS) 分析等位基因相似性。
• Runs of homozygosity (ROH)：连续性纯合片段分析（基因组中出现的连续不间断的纯合现象）。
• LD衰减速率。
• Fst分化指数。
• θπ核酸多样性。
• Tajma'D中性进化检验

image.png

tagSNPs用于GS

156,502个tagSNPs已经去除高度冗余的SNP，因此可作为GS应用的标记池。
使用rrBLUP模型比较了水稻9个性状的不同5个SNP集的准确性：

• set1：水稻44K芯片的原始29,434个SNPs。
• set2：156,502个tagSNPs与set1的交集，共1090个SNPs。
• set3：set1中随机选取的1090个SNPs。
• set4：根据set1中基因组距离（每350kb一个SNP）选取的1090个SNPs。
• set5：根据set1中随机基因组区域的1090个连续SNPs。

image.png

结论：从tagSNPs池中选择大约1000个SNPs可能是降低GS应用成本的方法。

fixedSNPs验证

• 图A：100kb窗口内计算的θπ和Fst，红点为潜在的强选择清除信号。
• 图B：潜在的强选择清除信号和其他基因组区域的Tajma'D值分布。
• 图C：栽培亚群中共有和特有的选择信号（括号内外数值分别为基因和GSEA term数目，GSEA使用PlantGSEA分析）。
• 图D：fixedSNPs的2556个材料的进化树。
• 图E：Affymetrix 700K芯片的880个材料的进化树。
• 图F：Illumina 44K芯片的351个材料的进化树。

image.png

barcodeSNPs指纹图谱

使用MinimalMarker算法筛选到的38个barcodeSNPs可作为特征来区分水稻品种（即指纹图谱）。
对2556份材料应用了7种经典机器学习算法来建模，10折CV，五类编码（10000, 01000, 00100, 00010, 00001）：

• 决策树DT
• K近邻KNN
• 朴素贝叶斯NB
• 人工神经网络ANN
• 随机森林RF
• 一对多法逻辑回归：one-vs-rest logistic regression(LR-O)
• 多元逻辑回归：multivariate logistic regression (LR-M)

5个水稻栽培亚群分类精度最佳的模型是LR-M（AUC为0.99）。并使用Affymetrix 700K芯片的880份材料进行了验证，证明其鲁棒性。

image.png

barcodeIndel

除了SNP，他们还做了Indel分析（<50bp）。

• 5152份材料种共4,217,174 raw Indels。
• 2556份材料按missing rate <0.01，MAF>=0.05 过滤后剩余109,898 Indels。
• 根据水稻6个亚群和籼稻内的6个子类，进一步鉴定62个亚群特异性Indels，即barcodeIndels。
• SR4R数据库中可下载供个性化分析。

SR4R数据库

• 地址：http://sr4r.ic4r.org/
• 在线分析包括亚群分类和指纹图谱：http://sr4r.ic4r.org/onlineTools/ml
• 提供了数据预处理、群体多样性分析和品种分类与鉴定的脚本：http://sr4r.ic4r.org/tools/bgp
• 提供了4类SNPs集的基因型和注释文件，还有Indel信息、脚本打包工具等：http://sr4r.ic4r.org/download

更加详细的方法可参考Method部分