plink软件是GWAS分析中常用的软件，它也是一个数据格式，plink里面有很多非常强大的功能，运算速度很快，是我日常分析中常用的软件之一。
之前写了一系列的GWAS教程，点击这里查看，这里继续进行。看到我的学习笔记帮助了一些同学，我也由衷的感到高兴。
这里，我将plink软件分为三部分：
- 格式转换
- 常用质控
- 文件提取
1. 格式转换
第一种常用的格式：plink格式
- 正常格式
map和ped：比如a.ped，a.map - 二进制文件
bim，bed，fam：比如a.bed, a.bim, a.fam
第二种常用的格式：vcf格式
第三种常用的格式：hapmap格式
1.1 plink正常格式转二进制格式
比如这里有plink格式的文件，前缀为a的plink文件：
$ ls
a.map a.ped
将其转化为二进制文件：b.bed, b.bim, b.fam
plink --file a --out b
结果：
$ ls b*
b.bed b.bim b.fam b.log
注意：
- 如果染色体超过23，比如30对染色体，需要设定
--chr-set 30 - 如果有非数字染色体，比如性染色体，需要设定
--allow-extra-chr - 常用的动物都有对应的参数，直接设定相关动物就行，比如牛的
--cow，下面是其它动植物的。如果没有对应的物种，直接设置染色体的条数以及允许非数字染色体即可。
--cow
--dog
--horse
--mouse·
--rice
--sheep
1.2 plink二进制格式转为正常格式（map和ped）
这里有plink格式的文件，前缀为b的plink二进制文件：
$ ls b*
b.bed b.bim b.fam b.log
将其转化文件：c.map, c.ped
plink --bfile b --recode --out c
注意：
- –bfile，因为输入文件b*为二进制，所以用–bfile，如果是一般格式，用–file即可
- –recode，要输出正常格式，所以用–recode指定，如果不加这个参数，默认是输出二进制文件
- –out，输出文件的前缀
结果：
$ ls *c*
c.hh c.log c.map c.ped
1.3 正常plink文件转为vcf文件
这里有plink格式的文件，前缀为c的plink二进制文件：
$ ls *c*
c.hh c.log c.map c.ped
将其转化文件：d.vcf
plink --file c --recode vcf --out d
注意：
- –file，用–file指定正常plink格式的文件
- –recode vcf，要输出vcf文件格式
- –out，输出文件的前缀
文件预览：
1.4 二进制plink文件转为vcf文件
和正常plink文件类似，除了--file 变为--bfile即可。
现有文件：
$ ls b*
b.bed b.bim b.fam b.log
将二进制文件转化为vcf文件：
plink --bfile b --recode vcf --out e
结果预览：
1.5 vcf文件转化为plink文件
转化为正常plink文件：
现有文件：
$ ls e.vcf
e.vcf
plink --vcf e.vcf --recode --out f
注意：
- –vcf 需要文件名完整，不能只写前缀，所以这里要写
--vcf e.vcf - –recode 保存plink文件
保存为二进制文件：
plink --vcf e.vcf --out g
结果：
$ ls g*
g.bed g.bim g.fam g.log
2. 常用质控
2.1 SNP缺失质控
无论是测序还是芯片，得到的基因型数据要进行质控，而对缺失数据进行筛选，可以去掉低质量的数据。如果一个个体，共有50万SNP数据，发现20%的SNP数据（10万）都缺失，那这个个体我们认为质量不合格，如果加入分析中可能会对结果产生负面的影响，所以我们可以把它删除。同样的道理，如果某个SNP，在500个样本中，缺失率为20%（即该SNP在100个个体中都没有分型结果），我们也可以认为该SNP质量较差，将去删除。当然，这里的20%是过滤标准，可以改变质控标准。
现有文件：
$ ls a*
a.map a.ped
某个SNP在样本中缺失大于10%，删除该SNP：--geno
plink --file a --geno 0.1 --recode --out re
某个在某个样本中，SNP缺失大于10%，删除该样本：--mind
plink --file a --mind 0.1 --recode --out re
2.2 最小等位基因频率过滤
最小等位基因频率怎么计算？比如一个位点有AA或者AT或者TT，那么就可以计算A的基因频率和T的基因频率，qA + qT = 1，这里谁比较小，谁就是最小等位基因频率，比如qA = 0.3, qT = 0.7， 那么这个位点的MAF为0.3. 之所以用这个过滤标准，是因为MAF如果非常小，比如低于0.02，那么意味着大部分位点都是相同的基因型，这些位点贡献的信息非常少，增加假阳性。更有甚者MAF为0，那就是所有位点只有一种基因型，这些位点没有贡献信息，放在计算中增加计算量，没有意义，所以要根据MAF进行过滤。
现有文件：
$ ls a*
a.map a.ped
某个SNP在的MAF小于0.01，那么该SNP删掉：--maf 0.01
plink --file a --maf 0.01 --recode --out re
2.3 哈温平衡过滤
「卡方适合性检验！」 ，一个群体是否符合这种状况，即达到了遗传平衡，也就是一对等位基因的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），符合遗传平衡定律.
现有文件：
$ ls a*
a.map a.ped
某个SNP在哈温平衡检验中p值小于1e-5，那么该SNP删掉：--maf 0.01
plink --file a --hwe 1e-5 --recode --out re
3. 文件提取
文件提取，可以提取plink个数中的样本信息，也可以提取特定的SNP位点信息。
3.1 样本提取--keep和-- remove
- –keep， 提取样本ID
- –remove，删除样本ID
提取样本文件的格式：
- 第一列：FID，家系ID
- 第二列：IID，个体ID
1328 NA06989
1377 NA11891
1349 NA11843
1330 NA12341
1344 NA10850
1328 NA06984
1463 NA12877
1418 NA12275
13291 NA06986
1418 NA12272
样本提取
plink --file a --keep id_sample.txt --recode --out re
完成。
$ wc -l re*
2 re.hh
32 re.log
1431211 re.map
10 re.ped
样本删除
plink --file a --remove id_sample.txt --recode --out re
完成。
3.2 SNP提取--extract和-- exclude
- –extract， 提取SNP ID
- –exclude，删除SNP ID
提取样本文件的格式：
- 一列：SNP名称ID
rs2185539
rs11240767
rs3131972
rs3131969
rs1048488
rs12562034
rs12124819
rs4040617
rs2905036
rs4245756
SNP提取
plink --file a --extract id_snp.txt --recode --out re
完成。
$ wc -l re*
179 re.hh
30 re.log
10 re.map
164 re.ped
可以看到，map共10行，共提取10个SNP
SNP删除
plink --file a --exclude id_snp.txt --recode --out re
完成。