【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏FreeBuf

详细分析使用Certutil解码的Office恶意软件

最近发现使用certutil进行解码的office恶意软件,这种恶意软件通过OLE机制落地通过BASE64编码的恶意软件,然后通过调用certutil来进行解码...

2414
来自专栏Y大宽

RNA-seq(5):序列比对:Hisat2

1 HISAT2官网下载 人类和小鼠的索引有现成的,HISAT2官网可以直接下载进行序列比对。如下图所示:选择hg19和mm10的index,文章中RNA-S...

1732
来自专栏哈雷彗星撞地球

iOS Bluetooth 打印小票(一)蓝牙打印机命令打印命令一览表打印的各个命令详解1.初始化命令2.打印命令3.行间距设置命令4.对齐方式设置5.字符设置命令6.钱箱控制命令7.按键控制命令8.

在iOS app中连接蓝牙打印机打印商品小票,在没有电脑只有手机的情况下,感觉非常实用,而且最近经常最近公司正好也在做这个功能,所以就研究了下。这一篇主要讲一下...

813
来自专栏Vamei实验室

Linux信号基础

Linux进程基础一文中已经提到,Linux以进程为单位来执行程序。我们可以将计算机看作一个大楼,内核(kernel)是大楼的管理员,进程是大楼的房客。每个进程...

1865
来自专栏吴伟祥

Linux命令缩写英文对照记忆(〇) 转

662
来自专栏散尽浮华

GlusterFS分布式存储系统中更换故障Brick的操作记录

前面已经介绍了GlusterFS分布式存储集群环境部署记录,现在模拟下更换故障Brick的操作: 1)GlusterFS集群系统一共有4个节点,集群信息如下: ...

5634
来自专栏about云

hbase Normalizer解决预分区错误,在不动数据的情况下完美解决热点问题

问题导读 1.对于预分区错误,hbase使用什么功能解决? 2.Region Normalizer的功能是什么? 3.在什么情况下运行Normalizer 比...

1061
来自专栏张善友的专栏

易学易用的Windows PowerShell

Windows PowerShell 是微软为 Windows 环境所开发的 shell 及脚本语言技术,这项全新的技术提供了丰富的控制与自动化的系统管理能力;...

1836
来自专栏流浪猫的golang

MongoDB 中文的全文索引

MongoDB 从3.2 版本以后添加了对中文索引的支持: 官网链接:https://docs.mongodb.com/manual/reference/t...

763
来自专栏Java成神之路

sublime text _注册码

打开 Sublime Text 3 的 “Help”–“Enter Licence”,然后根据版本选择输入下面的注册码。

1216

扫码关注云+社区