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

相关文章

来自专栏杨建荣的学习笔记

服务器硬件问题整理的一点总结 (r7笔记第70天)

之前写过一篇通过shell来监控磁盘坏块的文章 http://blog.itpub.net/23718752/viewspace-1872978/ 从使用情况来...

3807
来自专栏小狼的世界

IBM的双机切换HACMP方案介绍

HACMP,全称为IBM High Availablity Cluster Multiprocessing。

1103
来自专栏bboysoul

linux 服务器带宽测试脚本ZBench

很郁闷,今天我的vultr服务器的ip被ban了,无奈只能换服务器,今天给大家推荐一个vps的带宽测速脚本ZBench可以一键测试你的服务器到国内和国外的速度

3232
来自专栏SDNLAB

SDN中”软件”如何定义”网络”

SDN(Software Defined Network)软件定义网络,字面释义都说了是“软件”来定义“网络”,但有心之人会想:这个“软件”到底是如何定义了我们...

3637
来自专栏james大数据架构

NET实现微信公共平台上传下载多媒体文件

      举个例子,有人对着我们的公共微信号拍个照片发送过来,然后我们处理这个照片,比如进行ocr识别字(随后就会降到这个例子),或者人脸识别,或者拍照取证等...

2887
来自专栏张善友的专栏

在MongoDB中实现聚合函数

随着组织产生的数据爆炸性增长,从GB到TB,从TB到PB,传统的数据库已经无法通过垂直扩展来管理如此之大数据。传统方法存储和处理数据的成本将会随着数据量增长而显...

2587
来自专栏腾讯Bugly的专栏

黑科技:程序员如何打造属于自己的分体键盘

前言 作为一名程序员,键盘在手,天下我有啊,不整个高大上的键盘怎么提升B 格。之前一直想买个机械键盘,听说机械键盘敲代码时格外舒爽,实在是提升效率与 B 格的神...

3425
来自专栏24K纯开源

Premiere&After Effects的实时预览插件开发

一、介绍         Adobe Premiere和After Effects在影视编辑、渲染领域已经得到广泛应用。全景视频在相应工具拼接好后也可以导入Pr...

2349
来自专栏贾老师の博客

网络数据帧与 MTU

1744
来自专栏cmazxiaoma的架构师之路

蛋疼的ElasticSearch(四)之基本用法、高级查询

ElasticSearch系列文章: 1.蛋疼的ElasticSearch(一)之安装ElasticSearch 2.蛋疼的ElasticSearch(二)...

814

扫码关注云+社区