GWAS--Plink分析

概述：

 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的调用

1. 打开terminal，cd到解压到的文件夹（把plink拖入终端便可）
2. 以后每次使用须要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

4.  3  同时对个体和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个个体被移除。