内容

库(GWAS.BAYES)

1介绍

的GWAS。贝叶斯包提供了用于分析GWAS数据的统计工具。目前可用的功能允许分析自交物种，如水稻和拟南芥。

我们用于GWAS分析的贝叶斯管道分两个阶段选择重要的snp。第一阶段与典型的GWAS数据统计分析相同:该阶段拟合与SNP数量一样多的线性混合模型，每个模型包含一个SNP以及模型成分，以解释SNP之间的相关性。这些模型拟合产生每个SNP的p值，然后使用多次比较校正(如Bonferroni校正或Benjamini和Hochberg FDR校正)来获得一组显著SNP。然后，典型的GWAS分析将报告这组显著的snp。

GWAS。贝叶斯更进一步。第二阶段GWAS。贝叶斯执行贝叶斯模型搜索，其中模型由筛选阶段的一组snp组成。第二阶段允许snp相互竞争，并且通常会产生更小的snp集。第二阶段提供了更多的信息，哪些是最有希望的snp后续实验室实验。

在这个程度上GWAS。贝叶斯提供以下功能:

这篇小短文探讨了三个不同的例子，以展示中提供的函数GWAS。贝叶斯实现以上几点。我们感兴趣的三个案例研究是:

每个数据集都有一组snp以及表型。

2功能

中实现的函数GWAS。贝叶斯描述如下:

3.模型/模型的假设

完全理解包装GWAS。贝叶斯需要一些简单的模型细节。本包中使用的GWAS研究模型是

\[\begin{equation*} \textbf{y} = X \boldsymbol{\beta} + Z \textbf{u} + A \textbf{q} + \boldsymbol{\epsilon} \ \text{where} \ \boldsymbol{\epsilon} \sim N(\textbf{0}，\sigma^2 I) \ \text{and} \ \textbf{q} \sim N(\textbf{0}，\sigma^2 \tau K) \end{equation*}\]

在哪里

并非所有GWAS模型都必须由亲属矩阵(亲缘关系)或主成分(人口结构)组成，但在某些情况下，这些结构需要提供准确的结果。包含亲属关系组件的一个常见问题是此添加会增加计算成本。GWAS。贝叶斯利用高效混合模型关联(EMMA)技术(Kang et al. 2008)加快计算速度。为了得到准确的结果，还必须满足模型的假设。上面列出的模型假设\ \ (pmb{\ε}\)服从正态分布。小插图强调了如何评估这一假设。

4两个模拟数据集

本节探讨模拟自交物种的模拟数据答:芥。我们研究了两个模拟数据集，其中一个数据集没有亲属关系结构，而另一个数据集有。在两个模拟数据集中，前五个snp是影响响应的唯一snp，即前五个snp是唯一的非空snp。数据集是在假设所有snp相互独立的情况下创建的。

4.1用主成分测量人口结构的数据

以下4.4.1数据预处理

首先，我们将讨论原始数据的格式，以及如何使用中可用的函数转换数据GWAS。贝叶斯包中。

数据(“vignette_lm_dat”)头(vignette_lm_dat[1:10]) # >表型SNP1 SNP2 SNP3 SNP4 SNP5 SNP6 SNP7 SNP8 SNP9 # > 1 2.849069 G T C T C G G # > 2 2.931406 G T C T C G G # > 3 2.375999 G T T C G C 3.573713 G # > 4 G T C T C G G # > 5 2.631515 G T C T C G G 6 # > 3.322648 G T C T C G G
$表型snp <- vignette_lm_dat[，-1]

snp既可以作为字符串导入，也可以作为因子导入。

要为分析准备snp，必须将数据从原始状态(如上所示)转换为精细状态，以加快分析速度。首先，应用标准化函数将snp转换为{0,1}。这可以通过两种方式实现。第一种方法是按字母顺序排列;也就是说，如果这个字母在字母表中出现在第一个，那么这个字母就是1，而在字母表中出现在第二个的字母就是0。第二种方法是通过主要和次要等位基因，主要等位基因变成1次要等位基因变成0。注意，两种方法都会导致相同的显著snp;两者对结果的解释不同。默认方法为方法= "major-minor"可能的方法是方法= c(“长短”,“字母”)。

snp < -标准化(snp =单核苷酸多态性方法=“长短”,number_cores = 1)单核苷酸多态性(1:6,1:10]# > [1][2][3][4][5][6][7][8][9][10]# > [1]1 0 1 1 1 0 1 1 1 1 # > [2]1 0 1 1 1 0 1 1 1 1 # > [3]1 0 1 1 1 0 1 1 1 1 # > [4]1 0 1 1 1 0 1 1 1 1 # > [5]1 0 1 1 1 0 1 1 1 1 # > [6]1 0 1 1 1 0 1 1 1 1

与大多数GWAS研究一样，我们在每个物种/生态型/分类群的复制中汇总snp和表型，以加快计算速度。的aggregate_SNPs函数同时接受SNP矩阵和表型响应，并返回一个包含聚合SNP设计和聚合表型的列表。

aggregate_list <- aggregate_SNPs(SNPs = SNPs, Y = Y) SNPs <- aggregate_list$SNPs Y <- aggregate_list$Y

下一步是对标准化SNP数据进行整理，为下一步的研究做准备预选函数。如果不执行此步骤，则预选函数会失效，因为它不能给出数据集中所有物种/生态型/分类群中具有相同等位基因的SNP的p值。的level_function函数将清除所有次要等位基因频率小于或等于输入值的snp，从而清理标准化snp加。使所有小等位基因频率的snp大于零，设置Maf = 0。该函数将返回一个列表，其中第一个元素是标准化的SNP矩阵，第二个元素是TRUE和FALSE的向量，如果保留相应的SNP，则为TRUE，如果删除SNP，则为FALSE。MAF的默认值为Maf = 0.01。

level_list <- level_function(SNPs = SNPs, MAF = 0.01) SNPs <- level_list$SNPs level_list$SNPs_Dropped #> [1] "None" dim(SNPs) #> [1] 270 1000

查看上面的输出，我们可以看到没有小的等位基因频率小于0.01的snp，因此水平函数返回的新矩阵具有原始的1000个snp。

方法计算上述所有预处理步骤preprocess_SNPs函数。下面突出显示了一个示例。

fullPreprocess <- preprocess_SNPs(SNPs = vignette_lm_dat[，-1]， Y = vignette_lm_dat$Phenotype,MAF = 0.01,number_cores = 1, na。rm = FALSE) all.equal(fullPreprocess$SNPs,SNPs) #> [1] TRUE all.equal(fullPreprocess$Y,Y) #> [1] TRUE fullPreprocess$SNPs_Dropped;level_list$SNPs_Dropped #> [1] "None" #> [1] "None"

上面的结果突出了两种数据准备方法如何产生相同的结果。

其次，通过主成分的计算估计种群结构。这是GWAS研究中用于消除假阳性数量的常用工具。的pca_function函数需要三个输入:上面创建的标准化snp矩阵，您希望在模型中出现的主成分数量，以及如果您想绘制主成分解释的百分比变化。一定要利用plot_it变量，以充分了解需要多少主成分来控制误报。该函数将返回一个矩阵，该矩阵的行数与您输入的矩阵相同，列数等于主成分的数量。

principal_comp <- pca_function(SNPs = SNPs,number_components = 3, plot_it = TRUE)
主成分解释的响应变化百分比

图1:主成分解释的响应变化百分比

图1表明，一个主成分足以捕获SNP数据中的大部分变化。因此，我们使用下面的命令保存第一个主成分，以便以后进行分析。

Principal_comp <- as.matrix(Principal_comp [，1])

4.1.2筛选步骤

的预选函数是基于更典型的GWAS分析。这将着眼于每个SNP单独与相关的主要成分或潜在的亲属关系结构。可以使用p值或bic (frequentist = TRUE对于p值和frequentist = FALSEBICs);由此，可以对p值(controlrate)。bic只使用贝叶斯错误发现修正。

4.1.2.1Bonferroni调整

使用类型1的错误率。05。

显著_snps_bonf <-预选(Y = Y, SNPs = SNPs,number_cores = 1, principal_components = principal_comp,frequentist = TRUE, controlrate = "bonferroni"，threshold = 0.05,kinship = FALSE, info = FALSE) # bonferroni校正和(显著_snps_bonf $Significant) #五个显著snp #> [1] 5 which(显著_snps_bonf $Significant == 1) #> [1] 1 2 3 4 5 5

预筛选步骤发现前5个snp是显著的，它们是模拟数据集的真实snp集。

我们可以评估残差，以确保不需要表型反应的转化。调用resids_diag函数;该函数将计算具有预选步骤中显著snp的模型，并返回残差图和夏皮罗-威尔克正态性检验结果。

resids_diag(Y = Y,SNPs = SNPs, significant = Significant_SNPs_Bonf$ significant, kinship = FALSE,principal_components = principal_comp, plot_it = TRUE)

夏皮罗-威尔克正态性检验数据:残差值W = 0.9918, p值= 0.1391

Shapiro-Wilk检验的p值为0.1391。因此，没有证据表明残差不是正态分布。这表明高斯假设是满足的，没有表型的转化是有保证的。

的cor_plot函数创建一个图来检查显著snp之间的相关性。

cor_plot(SNPs = SNPs, significant = Significant_SNPs_Bonf$ significant, info = FALSE)

该图突出了SNP对之间的相关性，对角线显示了相同SNP之间的相关性。

4.1.2.2错误发现纠正

的GWAS。贝叶斯包包括其他控制多重比较的方法，包括Benjamini和Hochberg(1995)提出的虚假发现方法。要查看包用于频繁多次比较更正的方法的完整列表，请键入p.adjust ?并探索p.adjust.methods部分。在下面的例子中，使用了Benjamini-Hochberg校正(controlrate = "BH")，第一类错误率为0.05。

显著_snps_fdr <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, principal_components = principal_comp,frequentist = TRUE, controlrate = "BH"，threshold = 0.05,kinship = FALSE,info = FALSE) sum(显著_snps_fdr $Significant) # 5个显著snp #> [1] 5 which(显著_snps_fdr $Significant == 1) #> [1] 1 2 3 4 5

benjamin - hochberg FDR校正的筛选前步骤也正确地发现前5个snp是显著的。

4.1.2.3贝叶斯错误发现纠正

GWAS。贝叶斯也可以使用贝叶斯多重比较更正，如贝叶斯错误发现更正(Newton et al. 2004)。下面的例子对零假设和备选假设使用了相同的权重(nullprob和alterprob分别)。

显著_snps_bfdr <- preselection(Y = Y, snp = snp,number_cores = 1, principal_components = principal_comp,frequentist = FALSE,nullprob = .5, alterprob = .5,threshold = .05,kinship = FALSE,info = FALSE) sum(显著_snps_bfdr $Significant) # 5个显著snp #> [1] 5 which(显著_snps_bfdr $Significant == 1) #> [1] 1 2 3 4 5

贝叶斯错误发现纠正也正确地发现前五个snp是显著的。

4.1.3贝叶斯模型选择

postGWAS函数实现了GWAS的第二阶段。贝叶斯,which performs Bayesian model search where the set of candidate SNPs obtained from the screening phase compete with each other. The resulting best models typically contain a much smaller set of SNPs, which are the most promising SNPs for follow-up lab experiments.

下面的函数使用了带有Bonferroni校正的预筛选步骤的结果。注意，这是用重要的论点。

GA_results <- postGWAS(Y = Y,SNPs = SNPs,number_cores = 1, significant = Significant_SNPs_Bonf$ significant, principal_components = principal_comp,maxiterations = 100, runs_til_stop = 10,kinship = FALSE,info = FALSE) GA_results #> $Models #> [1] "Model 1: Y ~ SNP1 + SNP2 + SNP3 + SNP4 + SNP5后验概率#> SNP1 SNP2 SNP3 SNP4 SNP5 #> [1，] 1 1 1 1 1 1 1

详尽的搜索表明，最好的模型是前5个snp再次匹配这个模拟数据集的正确设置的模型。此外，该模型的后验概率为1，进一步证实了这5个snp是最有希望进行后续实验室实验的。

4.2用亲属矩阵测量人口结构的数据

另一个特点是GWAS。贝叶斯包装是塑造亲属关系的能力。亲属矩阵定义了分析中个体或物种之间的亲缘关系。包含随机效应与亲属协方差矩阵是一个非常流行的工具，以帮助减少假阳性在GWAS分析。我们将展示如何使用该函数计算亲属关系A.mat从rrBLUP包(Endelman 2011)利用我们的包裹和亲属关系结构。这个函数A.mat以标准化snp矩阵为参数，类似于pca_function以上。请注意，由于大多数对自交物种的GWAS分析具有相同生态型的复制，因此具有相同的SNP结构A.mat函数计算具有不同SNP结构的所有不同生态型之间的亲缘关系矩阵。

数据(“vignette_kinship_dat”)头(vignette_kinship_dat[1:10]) # >表型SNP1 SNP2 SNP3 SNP4 SNP5 SNP6 SNP7 SNP8 SNP9 # > 1 0.9370809 G C T C T T G # > 2 1.0153545 G C T C T T G # > 3 C 0.6490744 G T C T T 0.6758369 G # > 4 G C T C T T G # > 5 C 0.8911326 G T C T T G 6 # > 0.9906855 G C T C T T G

4.2.1数据预处理

按照上面列出的步骤:

$Phenotype snp <- vignette_kinship_dat[，-1] fullPreprocess <- preprocess_snp (SNPs = SNPs, Y = Y, MAF = 0.01, number_cores = 1,na。rm = FALSE) SNPs <- fullPreprocess$SNPs Y <- fullPreprocess$Y fullPreprocess$ snps_drop# > [1] "None"

计算亲属矩阵(k)。rrBLUP包中。我们可以用n.core论点A.mat。如上所述，rrBLUP包使用居中的IBS方法。

library(rrBLUP,quiet = TRUE) k <- A.mat(SNPs,n.;核心= 1)dim(k) #> [1] 270 270

4.2.2筛选步骤

用亲属关系结构来实现模型的变量亲属关系必须设置为亲属关系= k在哪里k是亲属矩阵。

4.2.2.1Bonferroni调整

使用0.05类型1错误率。

<- preselection(Y = Y, SNPs = SNPs,number_cores = 1, principal_components = FALSE,frequentist = TRUE,controlrate = "bonferroni"， threshold = 0.05,kinship = k, info = FALSE) sum(Significant_SNPs_Bonf$Significant)#四个显著snp #> [1] 5 which(Significant_SNPs_Bonf$Significant == 1) #> [1] 1 2 3 4 5 5

使用匹配正确设置的筛选前步骤，前5个snp是显著的(前5个snp被设置为非空)。

当存在亲属关系结构时，我们可以很容易地评估残差。只需设置亲属关系值等于亲属矩阵(亲属关系= k)。

resids_diag(Y = Y,SNPs = SNPs,significant = Significant_SNPs_Bonf$ significant, kinship = k,principal_components = FALSE,plot_it = TRUE)

夏皮罗-威尔克正态性检验数据:残差值W = 0.99422, p值= 0.3935

下面描述了控制多重比较的其他方法，以便读者了解结果的差异。

4.2.2.2错误发现纠正

使用类型1的错误率为0.05。

显著_snps_fdr <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, principal_components = FALSE,frequentist = TRUE,controlrate = "BH"， threshold = 0.05,kinship = k, info = FALSE) sum(显著_snps_fdr $Significant) # 6个显著snp #> [1] 5 which(显著_snps_fdr $Significant == 1) #> [1] 1 2 3 4 5

benjamin - hochberg FDR校正的筛选前步骤表明SNPs 1、2、3、4和5是显著的。这正确地识别了前5个snp。

4.2.2.3贝叶斯错误发现纠正

对零假设和备选假设使用相同的权重。

显著_snps_bfdr <- preselection(Y = Y, snp = snp,number_cores = 1, principal_components = FALSE,frequentist = FALSE,nullprob = .5, alterprob = .5,threshold = .05,kinship = k, info = FALSE) sum(显著snps_bfdr $Significant) #4显著snp #> [1] 5 which(显著snps_bfdr $Significant == 1) #> [1] 1 2 3 4 5

与Bonferroni校正一样，贝叶斯错误发现校正的预筛选步骤表明前5个snp是重要的。

4.2.3贝叶斯模型选择

我们现在使用benjamin - hochberg FDR决定的显著snp作为候选snp进行贝叶斯模型选择。为了在模型中包含亲属关系结构，我们再次使用论证亲属关系= k。

GA_results <- postGWAS(Y = Y,SNPs = SNPs,number_cores = 1, significant = Significant_SNPs_FDR$ significant, principal_components = FALSE,kinship = k, info = FALSE) GA_results #> $Models #> [1] "Model 1: Y ~ SNP1 + SNP2 + SNP3 + SNP4 + SNP5后验概率#> SNP1 SNP2 SNP3 SNP4 SNP5 #> [1，] 1 1 1 1 1 1 1

遗传搜索正确地提出了以前5个snp为最佳模型的模型。这个模型的后验概率是所有模型中最高的。这可以解释为在所有模型中，模型1是最好的模型。

5A. Thaliana分析

我们强调了一个真实的数据分析答:芥来自《植物细胞》杂志上发表的论文《盐胁迫下根构型重塑的遗传成分》(Julkowska et al. 2017)。本文研究了答:芥对盐的压力。本图研究的表型是在低盐胁迫条件下(75 mM NaCl)平均侧根长与主根长之比。GWAS分析在他们的论文的图3中突出显示。本文使用使用亲属矩阵控制亲缘关系的模型确定了4号染色体和5号染色体之间的显着snp区域。为了说明，本插图中包含的数据集包含4号染色体的最后500个snp, 5号染色体的前500个snp，以及随机选择的500个额外的snp。

5.1数据预处理

首先我们使用preprocess_SNPs数据预处理函数:

data("RealDataSNPs_Y") Y <- RealDataSNPs_Y$Phenotype snp <-子集(RealDataSNPs_Y,select = -c(Phenotype)) fullPreprocess <- preprocess_snp (SNPs = snp,Y = Y, MAF = 0.01,number_cores = 1,na。rm = FALSE) SNPs <- fullPreprocess$SNPs Y <- fullPreprocess$Y fullPreprocess$ snps_drop# > [1] 231 924 981 1372 1481

有两件事值得注意。首先，本例中使用的数据集没有复制，因此如果使用单个函数标准化,aggregate_SNPs,level_function,aggregate_SNPs功能将不需要。第二，注意SNPs_Droppedvalue返回的数字对应于由于次要等位基因频率低而从SNP矩阵中删除的列。

矩阵RealDataInfo提供了每个SNP的位置信息。每一列对应一个SNP，第一行给出染色体，第二行给出染色体内的位置。下面的代码显示了前6个snp的信息:

数据(“RealDataInfo”)头(RealDataInfo[1:6]) # >[1][2][3][4][5][6] # >染色体“1”“1”“1”“1”“1”“1”# >位置“88300”“164815”“225129”“292581”“400593”“901328”

因为从分析中删除了一些snp，我们更新了矩阵RealDataInfo:

RealDataInfo <- RealDataInfo[，- fullpreprocess $ snps_drop]

因为我们只提取了一些snp的数据，从这组snp中计算出的亲属矩阵将不能代表真正的亲属矩阵。因此，我们计算了整个单核苷酸多态性的亲属矩阵。现在我们将亲属关系矩阵引入到我们的分析中:

data("RealDataKinship")亲属关系<- as.matrix(RealDataKinship)

5.2筛选步骤

在屏幕前的步骤中信息参数可设置为info = RealDataInfo。这样，结果表将提供更多信息。

5.2.0.1Bonferroni调整

使用0.05类型1错误率。

Significant_SNPs <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, principal_components = FALSE,frequentist = TRUE,controlrate = "bonferroni"， threshold = 0.05,kinship = kinship，info = RealDataInfo) sum(Significant_SNPs$Significant)#11 Significant SNPs #> [1] 10 Significant_SNPs[Significant_SNPs$Significant == 1,c(1,2)] #>染色体位置#> 771 4 18521230 #> 801 4 18532166 #> 819 4 18537129 #> 820 4 18537303 #> 823 4 18537997 #> 827 4 18539430 #> 830 4 18540004 #> 832 4 18540257 #> 833 4 18540401 #> 852 4 18545882

4号染色体末端有10个显著的snp，这与论文的结果相符。

绘制残差来检验正态性假设，

resids_diag(Y = Y, SNPs = SNPs, significant = Significant_SNPs$ significant, kinship =亲属关系)

夏皮罗-威尔克正态性检验#> #>数据:残差#> W = 0.75853, p值< 2.212 -16

夏皮罗-威尔克检验的p值小于0.05，有力地证明残差不是正态分布，因此违反高斯假设。应对措施的转变是必要的;log变换返回一个可接受的Shapiro-Wilk p值。

显著snp <-预选(Y = log(Y)， snp = snp,number_cores = 1, principal_components = FALSE,frequentist = TRUE,controlrate = "bonferroni"， threshold = .05,kinship = kinship,info = RealDataInfo) sum(显著snp $Significant)#4显著snp #>[1] 4显著snp[显著snp $Significant == 1,c(1,2)] #>染色体位置#> 801 4 18532166 #> 820 4 18537303 #> 833 4 18540401 #> 852 4 18545882

有了新的对数变换反应，现在4号染色体上有4个重要的snp。

resids_diag(Y = log(Y)， SNPs = SNPs, significant = Significant_SNPs$ significant,kinship =亲缘关系)

夏皮罗-威尔克正态性检验数据:残差值W = 0.99359, p值= 0.1785

现在没有证据表明违反高斯假设。因此，我们可以相信结果。

5.3曼哈顿和QQ地块使用qqman包

如果你想想象曼哈顿的普通地块，这个包裹qqman（特纳(2014))是这些情节的绝佳组合。

library(qqman,quietly = TRUE) significant_snp <- chbind (significant_snp,paste0("rs"，1:nrow(significant_snp))) colnames(significant_snp) <- c(colnames(significant_snp)[1:4]，"SNPs") manhattan(significant_snp,chr = " chromosome "，bp = "Positions"， p = "P_values"，snp = "SNPs"，suggestiveline = FALSE, genomewideline = -log10(.05/nrow(significant_snp)))

曼哈顿地块的不同颜色突出了不同的染色体。

的qqman包也可以用这个函数为p值制作QQ图qq。该函数在x轴上绘制了-log10(p值)的期望值(在无SNP效应的零假设下)，在y轴上绘制了观测到的p值。偏离直线表明存在预测表型的snp。下面是如何使用qq功能:

qq (Significant_SNPs P_values美元)

由于这条线偏离直线，有证据表明存在预测低盐胁迫下平均侧根与主根长度之比的snp。

5.4贝叶斯模型选择

接下来，我们执行贝叶斯模型选择来确定，在这4个snp中，哪些是预测表型的重要因素。

GA_results < - postGWAS (Y =日志(Y), snp = snp, number_cores = 1,重大= Significant_SNPs重大美元,principal_components = FALSE,亲属关系=亲属信息= RealDataInfo) # >[1]”“删除列:3 # >[1]”“删除列:3 # >[1]”“删除列:3 # >[1]”“删除列:3 # >[1]“删除列:3”“删除列:4 # >[1]“删除列:4 # >[1]“删除列:4 # >[1]“删除列:3”“删除列:5”GA_results # > $模型# >[1]”模型1:y ~ SNP833后探模型2:y ~ SNP801型号3:y ~ SNP820型号4:y ~ SNP801 + SNP820后置问题模型5:y ~ SNP852后验概率型号6:y ~ SNP801 + SNP852后验问题型号7:y ~ SNP820 + SNP852型号8:y ~ SNP801 + SNP820 + SNP852$ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $

对于4个候选snp，有2^4 - 1 = 15个可能的模型。上述输出中没有显示15个模型的原因是，输出仅限于显示后验概率最高累积95%的模型。由于候选snp的数目为4(<= 12)，函数postGWAS执行穷举搜索并计算所有15个可能模型的后验概率。在这个例子中，模型2到8具有完全相同的后验概率。解释这一结果的一种方法是，这三个snp(801、820和852)提供了相同的信息。另一种看待这个结果的方式是，在模型中拥有1个SNP与拥有所有3个SNP没有区别。这可能表明这3个snp之间存在高度相关性。

贝叶斯模型搜索显示，这三个snp为表型提供了完全相同的信息。这就是为什么贝叶斯模型选择提出了7种不同的模型。它告诉我们拥有1个SNP和拥有全部3个SNP没有区别，它们提供相同的信息。在实验室环境下，这三个snp应该被观察，因为这些snp在同一个区域，这个区域应该是进一步研究的优先事项。

让我们用相关图进一步研究这个问题。

cor_plot(SNPs = SNPs[，Significant_SNPs$Significant == 1]， Significant = c(1,1,1,1)， info = RealDataInfo[，Significant_SNPs$Significant == 1])

事实上，相关图表明，3个最重要的snp对之间的所有相关性都等于1。为了进一步研究，可以用以下代码计算4个显著snp的次要等位基因频率:

1 - colMeans(SNPs[，Significant_SNPs$Significant == 1]) #> [1] 0.04268293 0.04268293 0.04878049 0.04268293

其中3个显著snp的次要等位基因频率为0.04268，相关系数为1，表明这些snp是相同的。我们可以通过查看4个snp之间的snp数量来了解基因的区域。

which(significant_snp $Significant == 1) #> [1] 801 820 833 852

前两个snp之间有18个snp，第二和第三个snp之间有12个snp，第三和第四个snp之间有18个snp，该区域的snp总数为52个。这些值仅用于本分析中使用的snp，而不包括由于次要等位基因频率低而减少的snp。为了更深入地了解这52个snp的次要等位基因频率汇总统计如下:

summary(1 - colMeans(SNPs[，801:852])) #>#> 0.04268 0.06402 0.16159 0.17765 0.28430 0.45732

这些汇总统计表明，在GWAS分析中发现的4个显著snp之间，有许多snp具有高等位基因变异性。因此，3个最显著的SNPs的两两相关等于1，并且被高度可变的SNPs分开，这一事实可能表明这3个最显著的SNPs倾向于共突变，从而可能与重要的生物学功能有关。因此，这三个最重要的snp似乎是后续实验室实验中非常有希望的候选者。

参考文献

Jeffrey B. Endelman, 2011。R包rrBLUP基因组选择的岭回归和其他内核。植物基因组4(3): 250-55。https://doi.org/10.3835/plantgenome2011.08.0024。

Julkowska, Magdalena M.， Iko T. Koevoets, Selena Mol, Huub Hoefsloot, Richard Feron, Mark A. Tester, Joost J.B. Keurentjes等。2017。盐胁迫下根构型重塑的遗传成分植物细胞29(12): 3198-3213。https://doi.org/10.1105/tpc.16.00680。

Kang, Hyun Min, Noah A. Zaitlen, Claire M. Wade, Andrew Kirby, David Heckerman, Mark J. Daly和Eleazar Eskin, 2008。模式生物关联图谱中种群结构的有效控制遗传学[j](3): 1709 - 1723。https://doi.org/10.1534/genetics.107.080101。

Newton, Michael A.， Amine Noueiry, Deepayan Sarkar和Paul Ahlquist, 2004。用半参数分层混合方法检测差异基因表达。生物统计学5(2): 155-76。https://doi.org/10.1093/biostatistics/5.2.155。

Stephen D. Turner, 2014。Qqman:一个使用Q-Q和曼哈顿图可视化Gwas结果的R包。bioRxiv。https://doi.org/10.1101/005165。