GWAS全称“全基因组关联分析”，使用统计模型找到与性状关联的位点，用于分子标记选择（MAS）或者基因定位。全基因组关联分析是对多个个体在全基因组范围的遗传变异多态性进行检测，获得基因型，进而将基因型与可观测的性状，即表型，进行群体水平的统计学分析，根据统计量或P值筛选出最有可能影响该性状的遗传变异。这次学习的教程是Plink做GWAS。
PLINK是一个免费的开源全基因组关联分析工具集，在SNP数据统计，过滤，GWAS分析中都可以用得上，而且计算速度很快。直接去百度上搜索plink就可以很容易就找到plink官网（http://www.cog-genomics.org/plink2）功能大概有以下几种：
- 数据管理： SNP数据格式的转换，合并两个或多个文件，提取SNP子集，以二进制文件格式压缩数据等。
- 质量控制的SNP数据统计： 计算丢失基因型率，等位基因，基因型频率，HWE测试，个体和个体对的近亲繁殖，IBS和IBD统计，LD区域计算等。
- GWAS关联分析
- Meta分析
GWAS 操作流程1-1：数据下载和plink配置
GWAS分析的两类性状：
- 分类性状（阈值性状，质量性状）：比如抗病性，颜色等等
- 连续性状（数量性状）：比如株高，体重，产量等等
GWAS的分析方法:
- 分类性状：logistic等等
- 连续性状：GLM，MLM模型等等
「一般线性模型（GLM）：
这里，SNP作为固定因子，可以考虑其它协变量（比如性别，PCA，群体结构等等）
「混合线性模型（MLM）：
- 固定因子：SNP + 可以考虑其它协变量（比如性别，PCA，群体结构等等），这里固定因子和前面的GLM一样
- 随机因子：亲缘关系矩阵（K矩阵或者A矩阵）
参考：
❝
教程代码和数据下载：https://github.com/MareesAT/GWA_tutorial/
❞
这个教程非常的经典，我看网上很多人推荐。
❝
相关的文章：https://onlinelibrary.wiley.com/doi/full/10.1002/mpr.1608
❞
教程中包括数据的过滤，SNP的过滤，样本的过滤，质控的标准等等，介绍的非常清楚，看完这篇文章，感觉plink的语法知识又增加了很多。
下载R语言和plink软件
- R：https://www.r-project.org/
- plink：http://zzz.bwh.harvard.edu/plink/ https://www.cog-genomics.org/plink2
plink的调用
- 打开terminal，cd到解压到的文件夹（把plink拖入终端即可）
- 之后每次使用需要cd到此路径，输入./plink再输入参数，如：
./plink --help
./plink --vcf Arab.vcf --allow-extra-chr --maf 0.05 --recode --out Arab
若需在终端中直接调用plink，则需要把plink写入全局路径
- vim .bash_profile
- #然后按i键进入编辑模式，输入以下语句，其中路径为解压plink的文件夹路径
- export PATH=/Users/Downloads/plink_mac_20200428:$PATH
- #输完后回车，按ESC
- ：w #表示保存
- ：q #表示退出
- source .bash_profile #更新即可
GWAS 操作流程2-1：缺失质控
「主要介绍」
❝
--geno筛选个体；--mind筛选SNP
❞
GWAS分析时，拿到基因型数据，拿到表型数据，要首先做以下几点：
- 1，查看自己的表型数据，是否有问题
- 2，查看自己的基因型数据，是否有问题
然后再进行建模，得到显著性SNP以及可视化结果。
「为何对缺失数据进行筛选？」
❝
无论是测序还是芯片，得到的基因型数据要进行质控，而对缺失数据进行筛选，可以去掉低质量的数据。如果一个个体，共有50万SNP数据，发现20%的SNP数据（10万）都缺失，那这个个体我们认为质量不合格，如果加入分析中可能会对结果产生负面的影响，所以我们可以把它删除。同样的道理，如果某个SNP，在500个样本中，缺失率为20%（即该SNP在100个个体中都没有分型结果），我们也可以认为该SNP质量较差，将去删除。当然，这里的20%是过滤标准，可以改变质控标准。下文中的质控标准是2%。
1. plink数据格式转化
数据使用上一篇的数据，因为数据是plink的bfile格式，二进制不方便查看，我们将其转化为文本map和ped的格式。
plink --bfile HapMap_3_r3_1 --recode --out test
结果生成：test.map test.ped
2. 查看基因型个体和SNP数量
wc -l test.map test.ped
可以看出，共有165个基因型个体，共有1447897个SNP数据。
「预览一下ped文件：」
「预览一下map文件：」
3. 查看一下个体缺失的位点数，每个SNP缺失的个体数
看一下描述：
--missing: Sample missing data report written to plink.imiss, and variant-based
missing data report written to plink.lmiss.
结果生成两个文件，分别是一个个体ID上SNP缺失的信息，另一个是每个SNP在个体ID中缺失的信息。
- 个体缺失位点的统计在plink.imiss中
- 单个SNP缺失的个体数在plink.lmiss.中
「个体缺失位点统计预览：」第一列为家系ID，第二列为个体ID，第三列是否表型缺失，第四列缺失的SNP个数，第五列总SNP个数，第六列缺失率。
「SNP缺失的个体数文件预览：」第一列为染色体，第二列为SNP名称，第三列为缺失个数，第四列为总个数，第五列为缺失率
4.对个体及SNP缺失率进行筛选
- 1， 如果一个SNP在个体中2%都是缺失的，那么就删掉该SNP，参数为：--mind 0.02
- 2，如果一个个体，有2%的SNP都是缺失的，那么就删掉该个体，参数为：--geno 0.02
4. 1 对个体缺失率进行筛选
「先过滤个体缺失率高于2%的SNP」
plink --bfile HapMap_3_r3_1 --geno 0.02 --make-bed --out HapMap_3_r3_2
「转化为map和ped的形式：」
plink --bfile HapMap_3_r3_2 --recode --out test
查看一下过滤后的行数，
「之前的为：」
1457897 test.map
165 test.ped
「现在的为：」
1430443 test.map
165 test.ped
可以看出，过滤了2万多个位点。
从当时的log日志里也可以看出这一点：
PLINK v1.90b6.5 64-bit (13 Sep 2018) www.cog-genomics.org/plink/1.9/
(
可以看到--geno，过滤了27454个位点。
4. 2 对SNP缺失率进行筛选
「过滤SNP缺失率高于2%的个体」
plink --bfile HapMap_3_r3_2 --mind 0.02 --make-bed --out HapMap_3_r3_3
「查看日志：」
「没有过滤掉个体，剩余：」个体：165 SNP：1430443
5 同时对个体和SNP的缺失率进行筛选
「两步合在一起，即过滤位点，又过滤个体：」
plink --bfile HapMap_3_r3_1 --geno 0.02 --mind 0.02 --make-bed --out HapMap_3_r3_5
plink日志：
可以看出，两者最终结果是一样的。
GWAS 操作流程2-2：性别质控
人类性别的信息的质控，主要是根据性染色上SNP的比值，判断性别，然后把性别错误的个体去掉或者更改性别信息。对其它物种参考意义不大，因为在动物中一般把性别信息的SNP去掉，植物中一般都是雌雄同体的，不涉及到这个问题，之所以会有这一篇，是因为原文中有这个信息，而且plink 也有--check-sex的参数，所以操作一下，留下笔记。
❞「原理：」检查性别差异。先验信息，女性的受试者的F值必须小于0.2，男性的受试者的F值必须大于0.8。这个F值是基于X染色体近交（纯合子）估计。不符合这些要求的受试者被PLINK标记为“PROBLEM”。
「上一步，去掉缺失信息后，现在有文件是过滤缺失后的文件：」
HapMap_3_r3_5.bed HapMap_3_r3_5.fam HapMap_3_r3_5.irem
HapMap_3_r3_5.bim HapMap_3_r3_5.hh HapMap_3_r3_5.log
1. 检查性别冲突
plink --bfile HapMap_3_r3_5 --check-sex
结果文件：plink.sexcheck第一列为家系ID，第二列为个体ID，第三列为系谱中的性别，第四列为SNP推断的性别，第五列是否正常，第六列为F值。
2. 提取错误的ID
我们使用grep过滤一下：根据STATUS列，如果有问题的话，为“PROBLEM”，我们可以根据这个关键词将有问题的行打印出来。
grep "PROBLEM" plink.sexcheck
1349 NA10854 2 1 PROBLEM 0.99
可以看出，个体NA10854是有问题的。
将相关错误的ID提取出来（家系ID，个体ID），之所以提取家系ID和个体ID，因为plink有参数remove可以根据ID进行筛选。
grep 'PROBLEM' plink.sexcheck | awk '{print $1,$2}' >sex_discrepancy.txt
我们将结果保存在sex_discrepancy.txt。
3. 使用remove去掉个体
plink --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMap_3_r3_6
当然，你也可以对个体进行判定填充，这是用--impute-sex就可以实现，这样的话那个错误的个体会根据统计量更改性别信息。这里我们选择的是删掉这个个体。
4. 过滤的关键词
去掉个体或者SNP，关键词不一样，容易混淆，这里总结一下。
「保留或去掉个体：」
--keep <filename>
--remove <filename>
--keep-fam <filename>
--remove-fam <filename>
「保留或去掉SNP：」
--extract ['range'] <filename>
--exclude ['range'] <filename>
GWAS 操作流程2-3：MAF过滤
上一次我们经过去掉缺失，去掉错误的性别信息，得到的文件为：
HapMap_3_r3_6.bed HapMap_3_r3_6.fam HapMap_3_r3_6.log
HapMap_3_r3_6.bim HapMap_3_r3_6.hh
这里，我们根据最小等位基因频率（MAF）去筛选。
「为什么要根据MAF去筛选？」
❝
最小等位基因频率怎么计算？比如一个位点有AA或者AT或者TT，那么就可以计算A的基因频率和T的基因频率，qA + qT = 1，这里谁比较小，谁就是最小等位基因频率，比如qA = 0.3, qT = 0.7， 那么这个位点的MAF为0.3. 之所以用这个过滤标准，是因为MAF如果非常小，比如低于0.02，那么意味着大部分位点都是相同的基因型，这些位点贡献的信息非常少，增加假阳性。更有甚者MAF为0，那就是所有位点只有一种基因型，这些位点没有贡献信息，放在计算中增加计算量，没有意义，所以要根据MAF进行过滤。
❞
1. 去掉性染色体上的位点
「思路：」
- 在map文件中选择常染色体，提取snp信息
- 根据snp信息进行提取
「提取常染色体上的位点名称：」
因为这里是人的数据，所以染色体只需要去1~22的常染色体，提取它的家系ID和个体ID，后面用于提取。
awk '{ if($1 >=1 && $1 <= 22) print $2}' HapMap_3_r3_6.bim > snp_1_22.txt
wc -l snp_1_22.txt
1399306 snp_1_22.txt
常染色体上共有1399306位点。
2. 提取常染色体上的位点
这里，用到了位点提取参数--extract
plink --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7
3. 计算每个SNP位点的基因频率
首先，通过参数--freq，计算每个SNP的MAF频率，通过直方图查看整体分布。可视化会更加直接。
plink --bfile HapMap_3_r3_7 --freq --out MAF_check
结果文件：MAF_check.frq预览：
4. 去掉MAF小于0.05的位点
plink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8
日志：
PLINK v1.90b6.5 64-bit (13 Sep 2018) www.cog-genomics.org/plink/1.9/
GWAS 操作流程2-4：哈温平衡检验
「什么是哈温平衡？」
❝
哈迪-温伯格(Hardy-Weinberg)法则 哈迪-温伯格（Hardy-Weinberg）法则是群体遗传中最重要的原理，它解释了繁殖如何影响群体的基因和基因型频率。这个法则是用Hardy,G.H (英国数学家) 和Weinberg，W.（德国医生）两位学者的姓来命名的，他们于同一年（1908年）各自发现了这一法则。他们提出在一个不发生突变、迁移和选择的无限大的随机交配的群体中，基因频率和基因型频率将逐代保持不变。---百度百科
❞
「怎么做哈温平衡检验？」
❝
「卡方适合性检验！」 ，一个群体是否符合这种状况，即达到了遗传平衡，也就是一对等位基因的3种基因型的比例分布符合公式：p2+2pq+q2=1,p+q=1,(p+q)2=1.基因型MM的频率为p2,NN的频率为q2,MN的频率为2pq。MN:MN：NN＝P2：2pq：q2。MN这对基因在群体中达此状态，就是达到了遗传平衡。如果没有达到这个状态，就是一个遗传不平衡的群体。但随着群体中的随机交配，将会保持这个基因频率和基因型分布比例，而较易达到遗传平衡状态。应用Hardy-Weinberg遗传平衡吻合度检验方法，把计算得到的基因频率代入，计算基因型平衡频率，再乘以总人数，求得预期值（e）。把观察数（O）与预期值（e）作比较，进行χ2检验。病例组和对照组的基因型分布的观察值和预期值差异无显著性（P>0.05），符合遗传平衡定律.
❞
「哈温平衡过滤和MAF过滤的区别？」
❝
之前，我对这两个概念有点混淆，后来明白过来了。这两个概念一个是对基因频率进行的筛选，一个是对基因型频率进行的筛选。对于一个位点“AA AT TT”，其中A的频率为基因频率，AA为基因型频率。MAF直接是对基因频率进行筛选，而哈温平衡检验，则是根据基因型推断出理想的（AA，AT，TT）的分布，然后和实际观察的进行适合性检验，然后得到P值，根据P值进行筛选。即P值越小，说明该位点越不符合哈温平衡。
❞
「两个目的：」
- 计算所有位点的哈温检测结果
- 删除SNP中不符合哈温平衡的位点
1. 计算所有位点的HWE的P值
plink --bfile HapMap_3_r3_8 --hardy
plink.hwe的数据格式：
- CHR 染色体
- SNP SNP的ID
- TEST 类型
- A1 minor 位点
- A2 major 位点
- GENO 基因型分布：A1A1, A1A2, A2A2
- O(HET) 观测杂合度频率
- E(HET) 期望杂合度频率
- P 哈温平衡的卡方检验P-value值
结果预览：
2. 提取哈温p值小于0.0001的位点
这里我们使用awk：
awk '{if($9 < 0.0001) print $0}' plink.hwe >plinkzoomhwe.hwe
共有123个位点，其中UNAFF为45个位点。
3. 设定过滤标准1e-4
plink --bfile HapMap_3_r3_8 --hwe 1e-4 --make-bed --out HapMap_3_r3_9
日志：
可以看到，共有45个SNP根据哈温的P值过滤掉了，和上面手动计算的一样。
GWAS 操作流程2-5：杂合率检验
一般自然群体，基因型个体的杂合度过高或者过低，都不正常，我们需要根据杂合度进行过滤。偏差可能表明样品受到污染，近亲繁殖。我们建议删除样品杂合率平均值中偏离±3 SD的个体。
❝
我的理解：非自然群体中，比如自交系，杂交种F1，这些群体不需要过滤杂合度。
❞
「参数过滤和手动过滤」plink有个特点，所有的过滤标准，都可以生成过滤前的文件，然后可以手动过滤，也可以用参数进行过滤。
- 比如：--missing生成结果，也可以用--geno和--mind过滤。
- 比如：--hardy生成结果，可以使用--hwe过滤
- 比如：--freq生成结果，可以用--maf过滤 但是杂合度--het，没有过滤的函数，只能通过编程去提取ID，然后用--remove去实现。
1. 计算杂合度
plink --bfile HapMap_3_r3_9 --het --out R_check
结果文件：
.het (method-of-moments F coefficient estimates)
Produced by --het.
A text file with a header line, and one line per sample with the following six fields:
FID Family ID # 家系ID
IID Within-family ID # 个体ID
O(HOM) Observed number of homozygotes # 实际纯合个数
E(HOM) Expected number of homozygotes # 期望纯合个数
N(NM) Number of non-missing autosomal genotypes # 总个数
F Method-of-moments F coefficient estimate #
所以，杂合度 = (N-O)/N
GWAS 操作流程2-6：去掉亲缘关系近的个体
「注意：」
❝
这里讲亲子关系的个体移除，不是必须要的，比如我们分析的群体里面有亲子关系的个体，想要进行分析，不需要做这一步的筛选。
❞
1. 计算pihat > 0.2的组合
plink --bfile HapMap_3_r3_10 --genome --min 0.2 --out pihat_min0.2
说明文档：
--genome invokes an IBS/IBD computation, and then writes a report with the following fields to plink.genome:
FID1 Family ID for first sample
IID1 Individual ID for first sample
FID2 Family ID for second sample
IID2 Individual ID for second sample
RT Relationship type inferred from .fam/.ped file
EZ IBD sharing expected value, based on just .fam/.ped relationship
Z0 P(IBD=0)
Z1 P(IBD=1)
Z2 P(IBD=2)
PI_HAT Proportion IBD, i.e. P(IBD=2) + 0.5*P(IBD=1)
PHE Pairwise phenotypic code (1, 0, -1 = AA, AU, and UU pairs, respectively)
DST IBS distance, i.e. (IBS2 + 0.5*IBS1) / (IBS0 + IBS1 + IBS2)
PPC IBS binomial test
RATIO HETHET : IBS0 SNP ratio (expected value 2)
2. 提取Z1大于0.9的个体
awk '{if($8>0.9) print $0}' pihat_min0.2.genome > zoom_pihat.genome
过滤出91个组合：
3.删除亲子关系的个体
plink --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11
日志：
可以看出，51个个体被移除。
