【直播】我的基因组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$gt\t$h{$F[2]}\n";    print OUT "$F[0]\t$F[1]\t$F[2]\t$F[3]\t$alt\t100\t.\tDP=100\tGT:DP:RO:AO\t$gt:100:$ro_po\n";}close FH;close OUT;

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

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-04-27

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏数据和云

奇异故障:SQL执行反复一慢两快

何剑敏 Oracle ACS华南区售后团队,首席技术工程师。多年从事第一线的数据库运维工作,有丰富项目经验、维护经验和调优经验,专注于数据库的整体运维。 今天...

2724
来自专栏张戈的专栏

《Vimtutor的中文版》快速学习Linux的vim命令

注:本资源收集与网络,版权归原作者所有。 下载地址 ---- = 欢 迎 阅 读 《 V I M 教 程 》 —— 版本 1.5 = vim 是一个具有很多命令...

3248
来自专栏北京马哥教育

Arch Linux的正确使用方法

谈起我的 Linux 学习之路,时间其实并不长。但是我却花了相对很少的时间,已经能达到把 Linux 当作自己的桌面系统的程度了。 Ubuntu 的体验令我有点...

3107
来自专栏知晓程序

什么?微信也可以听文章!看公众号再也不怕伤眼睛

春节假期刚刚结束,很多人(被迫)回到自己的工作岗位。双眼长时间盯着电脑屏幕,你是否也感觉到,眼睛越来越干涩,越来越酸?

754
来自专栏开源项目

码云小课堂 | 主流的开源协议有哪些?我们该如何选择?

主流的开源协议有哪些?我们该如何选择? License是软件的授权许可,里面详尽表述了你获得代码后拥有的权利,可以对别人的作品进行何种操作,何种操作又是被禁止的...

31510
来自专栏大数据挖掘DT机器学习

数据爬取、清洗到分析,如何徒手研究上海二手房房价

本文由作者:孙培培 原创投稿 声明:本文所公布代码及数据仅作学习用,若别有用途则后果自行承担。 提到上海,不得不提上海的高房价,最近一篇上海各市辖区均价的文章...

2776
来自专栏Android先生

从小白到独立开发Android和IOS两种平台app过程与总结

16年上半年在帮老师,帮外面随便做点东西以便得到些生活费养活自己。下半年去外面公司待了一段时间,然后选择回来帮自己做个东西,历时三个月,独自完成安卓及IOS版本...

541
来自专栏Python中文社区

用Python对鹿晗、关晓彤微博进行情感分析

專 欄 ❈大吉大利小米酱,Python中文社区专栏作者,Python爱好者,顽强地自学中,18线灵魂画手/段子手/脑洞女王。 简书: http://www....

2559
来自专栏知晓程序

帮你快速抢红包,微信聊天记录竟有这些隐藏操作? | 晓技巧

993
来自专栏Crossin的编程教室

工欲善其事必先利其器:用什么写Python?

通常来说,每个程序员都有自己趁手的兵器:代码编辑器。你要是让他换个开发环境,恐怕开发效率至少下降三成。然而,每个人对编辑器的喜好各不相同,甚至引发出诸如“神的编...

662

扫描关注云+社区