plink软件初体验2--常用参数

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删掉：--hwe 1e-5」

 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

完成。