GWAS -4 VCF格式文件转为Plink文件

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了GWAS -4 VCF格式文件转为Plink文件相关的知识,希望对你有一定的参考价值。

参考技术A 参考: https://cloud.tencent.com/developer/article/1556166

plink1.9版本支持转化为VCFv4.2格式

plink2.0版本支持转化为VCFv4.3格式

两个版本用到的命令不一样

对于plink1.9版本,转化为vcf文件的命令行为:

生成的vcf为4.2版本

对于plink2.0版本,转化为vcf文件的命令行为:

生成的vcf为4.3版本

参考: https://www.cnblogs.com/chenwenyan/p/8574237.html

链接: https://www.jianshu.com/p/8b4e7b3b7f5e

vcf 转为 ped/map

vcftools --vcf snp.vcf --plink --out snp

plink --vcf snp.vcf --recode --out snp
ped和map文件是Plink的基本格式。

ped文件包含以下几列:

第一列:Family ID。

第二列:Individual ID。自然群体这列和Family ID是一样的。

第三列:Paternal ID。未提供信息的话这列为0。

第四列:Maternal ID。未提供信息的话这列为0。

第五列:Sex。未提供信息的话这列为0。

第六列:Phenotype。一般来说,直接拿vcf转换的话这列为-9,也就是缺失。

第七列开始就是个体在每个标记位点的基因型。

map文件包含以下几列:

第一列:染色体编号。

第二列:SNP编号。

第三列:遗传距离。未提供信息的话这列为0。

第四列:物理位置。

ped/map 与 tped/tfam 格式互换

plink --file snp --recode --transpose --out snp_test

plink --tfile snp_test --recode --out snp
ped/map 与 bed/bim/fam互换

plink --file snp --make-bed --out snp_test

plink --bfile snp_test --recode --out snp
tped/tfam 与 bed/bim/fam互换

plink --tfile snp --make-bed --out snp_test

plink --bfile snp_test --recode --transpose --out snp
bed/bim/fam 转为 vcf

plink --bfile snp --export vcf --out snp_test
常用的Plink格式转换就是这些,大家可以根据自己实际需要相互转换。

因为PLINK默认的设置是人的染色体, 所以动物中,我们应该设置
--chr-set 19 # 猪
已有的选择:
--cow
--dog
--horse
--mouse
--rice
--sheep
参考: https://zhuanlan.zhihu.com/p/109071456

Vcf文件格式

 

Vcf文件格式是GATK钟爱的表示遗传变异的一种文件格式。

 

就拿GATK给出的vcf例子说明吧,下面这个文件只表示了一个完整vcf文件的前几个SNP。

 

技术分享

看上去确实有点复杂,那就把它分为两部分看吧,第一部分把他归为说明文件,就是每一列最前面有2个#符号的那些列所提到的就是为了解释下面“正文”INFO列中可能要出现的一些tags和和FORMAT列中对基因型的表示。第二部分可以归为下面的内容:

 

#CHROM  POS ID      REF ALT QUAL    FILTER  INFO    FORMAT  NA12878

chr1    873762  .       T   G   5231.78 PASS    AC=1;AF=0.50;AN=2;DP=315;Dels=0.00;HRun=2;HaplotypeScore=15.11;MQ=91.05;MQ0=15;QD=16.61;SB=-1533.02;VQSLOD=-1.5473 GT:AD:DP:GQ:PL   0/1:173,141:282:99:255,0,255

chr1    877664  rs3828047   A   G   3931.66 PASS    AC=2;AF=1.00;AN=2;DB;DP=105;Dels=0.00;HRun=1;HaplotypeScore=1.59;MQ=92.52;MQ0=4;QD=37.44;SB=-1152.13;VQSLOD= 0.1185 GT:AD:DP:GQ:PL  1/1:0,105:94:99:255,255,0

chr1    899282  rs28548431  C   T   71.77   PASS    AC=1;AF=0.50;AN=2;DB;DP=4;Dels=0.00;HRun=0;HaplotypeScore=0.00;MQ=99.00;MQ0=0;QD=17.94;SB=-46.55;VQSLOD=-1.9148 GT:AD:DP:GQ:PL  0/1:1,3:4:25.92:103,0,26

chr1    974165  rs9442391   T   C   29.84   LowQual AC=1;AF=0.50;AN=2;DB;DP=18;Dels=0.00;HRun=1;HaplotypeScore=0.16;MQ=95.26;MQ0=0;QD=1.66;SB=-0.98 GT:AD:DP:GQ:PL 0/1:14,4:14:60.91:61,0,255

 

CHROM:       表示变异位点是在哪个contig 里call出来的,如果是人类全基因组的话那就是chr1…chr22,chrX,Y,M了。

 

POS:          变异位点相对于参考基因组所在的位置,如果是indel,就是第一个碱基所在的位置。

 

ID:            如果call出来的SNP存在于dbsnp数据库里,就会显示相应的dbsnp里的rs编号。

 

REF和REF:     在这个变异位点处,参考基因组中所对应的碱基和研究对象基因组中所对应的碱基。

 

QUAL:         可以理解为所call出来的变异位点的质量值。Q=-10lgP,Q表示质量值;P表示这个位点发生错误的概率。因此,如果想把错误率从控制在90%以上,P的阈值就是1/10,那lg(1/10)=-1,Q=(-10)*(-1)=10。同理,当Q=20时,错误率就控制在了0.01。

 

FILTER:       理想情况下,QUAL这个值应该是用所有的错误模型算出来的,这个值就可以代表正确的变异位点了,但是事实是做不到的。因此,还需要对原始变异位点做进一步的过滤。无论你用什么方法对变异位点进行过滤,过滤完了之后,在FILTER一栏都会留下过滤记录,如果是通过了过滤标准,那么这些通过标准的好的变异位点的FILTER一栏就会注释一个PASS,如果没有通过过滤,就会在FILTER这一栏提示除了PASS的其他信息。如果这一栏是一个“.”的话,就说明没有进行过任何过滤。

 

到现在,我们就可以解释上面的例子了:

 

chr1:873762是一个新发现的T/G变异,并且有很高的可信度(qual=5231.78)。

chr1:877664是一个已知的变异为A/G 的SNP位点,名字rs3828047,并且具有很高的可信度(qual=3931.66)。

chr1:899282是一个已知的变异为C/T的SNP位点,名字rs28548431,但可信度较低(qual=71.77)。

chr1:974165是一个已知的变异为T/C的SNP位点,名字rs9442391,但是这个位点的质量值很低,被标

 成了“LowQual”,在后续分析中可以被过滤掉。

 

Vcf文件看起来很复杂,挺吓人的样子,但是里面大部分都是一些tags,而这些tags基本上都是在VASR中过滤用的,能够理解每个tags的意思最好,如果实在不理解也就不用管了。其实最关键的信息也就是那么几列:

 

chr1    873762      .       T   G   [CLIPPED]  GT:AD:DP:GQ:PL    0/1:173,141:282:99:255,0,255

chr1    877664  rs3828047   A   G   [CLIPPED]  GT:AD:DP:GQ:PL    1/1:0,105:94:99:255,255,0

chr1    899282  rs28548431  C   T   [CLIPPED]  GT:AD:DP:GQ:PL    0/1:1,3:4:25.92:103,0,26

 

其中最后面两列是相对应的,每一个tag对应一个或者一组值,如:

 

chr1:873762,GT对应0/1;AD对应173,141;DP对应282;GQ对应99;PL对应255,0,255。

 

GT:    表示这个样本的基因型,对于一个二倍体生物,GT值表示的是这个样本在这个位点所携带的两个等位基因。0表示跟REF一样;1表示表示跟ALT一样;2表示第二个ALT。当只有一个ALT 等位基因的时候,0/0表示纯和且跟REF一致;0/1表示杂合,两个allele一个是ALT一个是REF;1/1表示纯和且都为ALT;

 

AD:    对应两个以逗号隔开的值,这两个值分别表示覆盖到REF和ALT碱基的reads数,相当于支持REF和支持ALT的测序深度。

 

DP:    覆盖到这个位点的总的reads数量,相当于这个位点的深度(并不是多有的reads数量,而是大概一定质量值要求的reads数)。

 

PL:      对应3个以逗号隔开的值,这三个值分别表示该位点基因型是0/0,0/1,1/1的没经过先验的标准化Phred-scaled似然值(L)。如果转换成支持该基因型概率(P)的话,由于L=-10lgP,那么P=10^(-L/10),因此,当L值为0时,P=10^0=1。因此,这个值越小,支持概率就越大,也就是说是这个基因型的可能性越大。

 

GQ:   表示最可能的基因型的质量值。表示的意义同QUAL。

 

举个例子说明一下:

 

chr1    899282  rs28548431  C   T   [CLIPPED]  GT:AD:DP:GQ:PL    0/1:1,3:4:25.92:103,0,26

 

在这个位点,GT=0/1,也就是说这个位点的基因型是C/T;GQ=25.92,质量值并不算太高,可能是因为cover到这个位点的reads数太少,DP=4,也就是说只有4条reads支持这个地方的变异;AD=1,3,也就是说支持REF的read有一条,支持ALT的有3条;在PL里,这个位点基因型的不确定性就表现的更突出了,0/1的PL值为0,虽然支持0/1的概率很高;但是1/1的PL值只有26,也就是说还有10^(-2.6)=0.25%的可能性是1/1;但几乎不可能是0/0,因为支持0/0的概率只有10^(-10.3)=5*10-11

 

到现在为止,基本上就把vcf文件介绍完了,如果想要了解更多的关于tags的信息的话,可以参考《GATK使用方法详解》第四部分。

 

作者 葡萄糖

邮箱 [email protected]



以上是关于GWAS -4 VCF格式文件转为Plink文件的主要内容,如果未能解决你的问题,请参考以下文章

2021-03-17 在linux上将vcf文件转plink的格式bed,bim,fam

将phasing的vcf文件转化为plink格式时保持phasing状态

VCF格式转换(plink, hapmap)

生物数据格式 - vcf/bcf

plink 进行PCA分析

全基因组关联分析(Genome-Wide Association Study,GWAS)流程