【直播】我的基因组72：把基因检测芯片数据转为vcf格式

这个需求比较少见，主要是因为有很多朋友都做了基因检测芯片数据，而芯片检测的结果只有4列，分别是dbSNP数据库ID号，染色体，坐标，还有基因型。如下：

head -100 wegene.txt |tail rs9442387    1   1110586 CCrs191284590    1   1115954 CCrs1320565    1   1119858 TCrs116321663    1   1120377 TTrs11260549    1   1121794 AGrs9729550    1   1135242 ACrs3819001    1   1138913 TTrs201786281    1   1140851 CCw01001152631    1   1152631 CCrs2887286    1   1156131 CC

但是呢，大部分的基因检测结果注释都是基于vcf文件的，vcf文件的详细介绍，我们以前讲过，就是

【直播】我的基因组28-必须要理解vcf格式记录的变异位点信息

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  unknown

这10列。

要想把基因检测芯片数据转为vcf格式就需要在充分理解vcf的基础上面再增加几个信息。

因为基因芯片的结果里面没有参考碱基是什么的信息，只有基因型，所以我们没办法判断纯合杂合或者突变。

至于 QUAL，这里统一给100，filter统一给一个占位符即可， INFO就给一个测序深度， FORMAT就给 GT:DP:RO:QR:AO信息吧，具体介绍如下：

##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth at the locus">##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=RO,Number=1,Type=Integer,Description="Reference allele observation count"> ##FORMAT=<ID=AO,Number=A,Type=Integer,Description="Alternate allele observation count">

这里我们还是借用dbSNP数据库文件，如下：

head  ~/annotation/variation/human/dbSNP/dbsnp.pos 1    10019   rs775809821 TA  T1    10055   rs768019142 T   TA1    10108   rs62651026  C   T1    10109   rs376007522 A   T1    10128   rs796688738 A   AC1    10139   rs368469931 A   T1    10144   rs144773400 TA  T1    10146   rs779258992 AC  A1    10150   rs371194064 C   T1    10165   rs796884232 A   AC

这个文件我以前讲过：

【直播】我的基因组（六）:变异位点注释数据库的准备

那么很简单的一个perl程序就可以达成这个转换啦：

open FH,"wegene.txt";while(<FH>){    chomp;    @F=split;    next if /^#/;    $h{$F[0]}=$F[3] if /[ATCG]{2}/; }close FH;open FH,"/home/jianmingzeng/annotation/variation/human/dbSNP/dbsnp.pos";open OUT,">wegene.vcf";print OUT '##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth at the locus">'."n";print OUT '##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">'."n"; print OUT '##FORMAT=<ID=RO,Number=1,Type=Integer,Description="Reference allele observation count">'."n"; print OUT '##FORMAT=<ID=AO,Number=A,Type=Integer,Description="Alternate allele observation count">'."n";print OUT '#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  unknown'."n";while(<FH>){     chomp;    @F=split;    next unless exists $h{$F[2]};    next if  $h{$F[2]} eq $F[3].$F[3];     $gt=$h{$F[2]} eq substr($h{$F[2]},0,1)x2 ?"1/1":"0/1";    $alt=substr($h{$F[2]},0,1) eq $F[3] ? substr($h{$F[2]},1,1): substr($h{$F[2]},0,1);    $ro_po=$gt eq '1/1' ? "0:100" : "50:50" ;    #print "$_t$gtt$h{$F[2]}n";    print OUT "$F[0]t$F[1]t$F[2]t$F[3]t$altt100t.tDP=100tGT:DP:RO:AOt$gt:100:$ro_pon";}close FH;close OUT;

运行完毕就可以打开我们转换好的vcf文件，如下所示：