搜索
查看: 2169|回复: 0

把vcf文件转为23andme芯片基因型文件

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-2-21 00:00:40 | 显示全部楼层 |阅读模式
这个难度有点高,所以我放在了这里,不知道有没有人可以做出来!

首先,找到一个全基因组测序的结果vcf文件,这个很容易
然后找到一个23andme的芯片基因型文件,我这里给大家:链接:https://pan.baidu.com/s/1hsLYOYS 密码:gsxd

把这两个文件搞清楚,就大概知道如何转换了,我这里给出思路和代码:


第一步,把芯片设计的rsID全部拿出来
第二步,根据rsID从我的VCF文件中挑取位点,并赋予纯合杂合基因型
第三步,去dbSNP数据库文件里面映射我VCF文件没有记录的点为野生型

[Perl] 纯文本查看 复制代码
perl -alne '{print if /^rs/}' dm_23andme_v3_110219.txt  |cut -f 1 >23andme.rsID.list
cat ../variation/autochr.highQuali.dbsnp.vcf  23andme.rsID.list |perl -alne '{if($F[2]=~/^rs/){if(/1\/1/){$gt=$F[4].$F[4]}else{$gt=$F[3].$F[4]};$h{$F[2]}="$F[0]\t$F[1]\t$gt" }  print "$_\t$h{$_}" if /^rs/}' >my_23andme.1.txt
zcat ~/annotation/variation/human/dbSNP/All_20160601.vcf.gz |perl -alne 'BEGIN{ open FH,"my_23andme.1.txt";while(<FH>){chomp;@F=split;if(/^rs/){ $pos{$.}=$_;if($F[3]){$h{$F[0]}=$_}else{$tmp{$F[0]}=1}  }} }{if(exists $tmp{$F[2]}){ $tmp{$F[2]}="$F[0]\t$F[1]\t$F[2]$F[2]"  }}END{foreach(sort{$a<=>$b} keys %pos){ if(exists $h{$pos{$_}} ){$value=$h{$pos{$_}}}else{$value=$tmp{$pos{$_}} } ;print "$pos{$_}\t$value" }}' >my_23andme.2.txt

说老实话,这个代码,我自己都看不下去了!

不过,需求是解决了!






上一篇:用dbsnp数据库进行基因组版本坐标转换
下一篇:rMATS 报错:AttributeError: 'csamtools.AlignedRead' object has no attribut...
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 08:29 , Processed in 0.037026 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.