GWAS分析<五>之群体结构分析一
通过前面的几篇推文，相信大家对GWAS的质控和分析原理已经初步了解。接下来，我们需要将之前学习到的知识运用到实践中。为了更好的进行实践，选择一个好的参考数据集就就成为了一个必不可少的环节。千人基因组数据集，作为一个旨在绘制最详尽的、最有医学应用价值的人类基因组遗传多态性图谱数据集，就成为了我们这次分析的首选参考数据集。本篇推文，将通过对千人基因组数据集进行深入而细致的GWAS分析，挖掘出不同种属之间的进化关系。那么，准备好了么，让我们开始GWAS应用分析之旅！！！
一 GWAS分析之旅
1.1 前人基因组的下载与数据整理
实际上，本教程并未下载人体所有染色体上的snp数据，仅下载了第21号染色体上的snp数据。这是因为如果所有染色体上的snp数据将近有T级别的数据，下载和分析都需要耗费大量时间。但是，我们分析所用21号染色体数据仅8个G的数据左右，大大加快了分析相关的进程。这也对应之前的推文中<GWAS分析<二>之数据质控的原理>所提到的如何大家如何使用最少的数据完成本次项目的分析的策略。
下文的两行代码分别为将千人基因组文件（vcf）转为plink格式（bed、fam、bim），并为缺失rs标识符的SNP分配唯一的标识符。虽然，对于千人基因组而言，数据中缺失的rs标识符不是问题。
1.2 移除不符合标准的snp类型
分别通过‘--geno 0.2’、‘--mind 0.2’、‘--geno 0.02’、‘--mind 0.02’和‘--maf 0.05’选项完成对符合标准的snp进行过滤。过滤的标准详见之间的推文。
1.3 提取snp变异
一方面，本教程将从千人基因组数据集中提取HapMap数据集中存在的snp变异，另一方面，本教程也将从HapMap数据集中提取千人基因组数据集中存在的snp变异。
1.4 hapmap文件的构建
利用上一步提取的HapMap_MDS.map文件构建hapmap.txt，该文件每行包含一个SNP id和物理位置。随后，利用plink软件将hapmap.txt中snp信息更新1kG_MDS6中，输出到1kG_MDS7文件中。
1.5 等位基因的校准
本步骤是确保HapMap和1000基因组项目数据集中的参考等位基因被重新指定或调整。
1.6 链校正问题
本步骤将检查潜在的链校正问题，并为不应该存在的A1等位基因分配生成警告输出。
1.7 snp翻转
本步骤翻转（正义或反义链上）SNP以解决链校正问题。本步骤将输出SNP标识符并删除重复项，再生成一个包含812个SNP的文件。文件中的这些snp来自于两个文件之间不对应的SNP。
1.8 snp再次检查
检查翻转后仍有问题的SNP。
1.9 去除有问题的snp
本步骤将从HapMap和1000个基因组中删除有问题的SNP。上一步骤生成了一个42个SNP的列表，这些SNP在翻转和参考等位基因被重新指定或调整后，在HapMap和1000基因组数据集之间造成了84处差异。因此，本步骤将分别从两个数据集中删除42个有问题的SNPs，并合并HapMap与千人基因组数据。值得注意，虽然HapMap和千人基因组的数据集之间存在样本重叠，但是这对本次分析并不重要。在由千人基因组数据上锚定的HapMap-CEU数据上执行MDS。
1.10 种群异常个体去除
本步骤将人口标识转换为更高级的人口标识（即AFR、AMR、ASN和EUR），并创建一个自己的数据文件。通过对结果进行可视化，显示我们的数据属于欧洲1000个基因组的数据组（输出文件是MDS.pdf）。并且，出于教程训练的目的，本步骤给出脚本也将过滤出人口分层异常值，从而为下一步骤生成适当的文件，用以排除种群中异常的个体。本步骤将在HapMap数据中选择低于过滤阈值的个体。需要注意的是，异常水平不是固定的阈值，而是基于前两个维度的可视化来确定（读者可自行确定）。
1.11 plink协变量文件创建
本步骤将在单倍体图谱数据中提取这些个体（排除异常值的个体）。注意，由于我们的单体型图数据确实包含任何种族异常值，因此在这一步中没有删除任何个体（如果我们的数据包括超出我们设定的阈值的个人，那么这些个体将被删除）。随后，我们将基于MDS创建协变量。注意，仅对没有种族异常值的HapMap数据执行MDS分析。更改文件的格式。将mds文件转换为plink协变量文件。注意：10个MDS维度的值随后被用作关联分析中的协变量。
进行到这里，本公众号也要恭喜读者，实现了人口分层数据的质控工作！！！
1.12 关联分析
对于关联分析，本步骤使用的文件来自于上一步骤中（人口分层）中生成的文件（HapMap_3_r3_13.bed、HapMap_3_r3_13.bim、HapMap_3_r3_13.fam和covar_mds.txt）。需要注意的是，本教程中的--assoc选项不允许校正主成分（PC）/MDS成分等协变量，这使得它不太适合关联分析。因此，本步骤将使用10个主成分（covar_MDS. txt）作为协变量进行logistic分析。
本步骤使用ide-covar选项仅在输出文件中显示SNP的相加结果。本步骤也将删除NA值，因为这些值会在后续生成图片时引发问题。需要注意的是，本步骤GWAS分析中所产生的结果将在随后的步骤用于可视化分析，显示数据集中是否包含基因组范围内的重要SNP。另外，如果是定量结果测定，则本步骤的选项--logistic应替换为--linear。使用--assoc选项也可以用于定量结果测量。除了常规的全基因组显著性阈值5.0E-8之外，还有多种方法可以处理多重检测。
1.13 p值计算
上一步给出了Bonferroni校正的p值，以及FDR和其他文件。本步是一个计算密集的步骤。本教程相对应的文章中描述了这种方法的进一步优缺点，它可以用于关联和处理多个测试(https://www.ncbi.nlm.nih.gov/pubmed/29484742)为了减少计算时间，我们只对22号染色体的一部分SNP进行了测试。EMP2 列为多次测试提供了校正的p值。根据上述步骤生成的SNP子集筛选您的bfile。执行1000000次测试。按从最低到最高的p值顺序排列数据。检查有序排列结果
1.14 曼哈顿图和QQ图的绘制
本步骤的话主要是通过可视化脚本绘制曼哈顿图和QQ图。
二 惯例小结
恭喜读者，到了这里，我们共同学习了如何对GWAS数据进行质控以及如何运用GWAS进行群体分析。事实上，就如同开篇所提到的那样，GWAS的分析远远不仅如此，他还可以与其他分析进行联合，进一步挖掘信息。如近期发表在Nature medicine上的一篇文章就利用MR整合GWAS和eQTL发现复杂疾病的causal genes（具体内容可见文献Actionable druggable genome-wideMendelian randomization identifies repurposing opportunities for COVID-19）。因此，读者可以在GWAS分析的基础上结合其他组学或者实验，实现文章在质上的飞跃。也希望各位读者天天有大paper。
本公众号开发的相关软件，Multi-omics Hammer软件和Multi-omics Visual软件欢迎大家使用。文末是本公众号在其他平台的账户，也欢迎大家关注并多提意见。
简书：WJ的生信小院
公众号：生信小院
博客园：生信小院
最后，也欢迎各位大佬能够在本平台上：1传播和讲解自己发表的论文；2：发表对某一科研领域的看法；3：想要达成的合作或者相应的招聘信息；4：展示自己以寻找博后工作或者博士就读的机会;5:博导提供博后工作或者博士攻读机会，都可以后台给笔者留言。希望本平台在进行生信知识分享的同时，能够成为生信分析者的交流平台，能够实现相应的利益互补和双赢（不一定能实现，但是梦想总得是有的吧）。
另外，怎么说呢，投币也可，不强求，但奢求。
