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删掉:--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
完成。
标签:文件,初体验,plink,vcf,--,re,SNP 来源: https://blog.csdn.net/yijiaobani/article/details/110006924