[直播]我的基因组 72:将基因检测微阵列数据转换为 vcf 格式
最编程
2024-03-18 20:27:36
...
这个需求比较少见,主要是因为有很多朋友都做了基因检测芯片数据,而芯片检测的结果只有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文件,如下所示:
下一篇: dbglog