搜索
查看: 2787|回复: 0

用dbsnp数据库进行基因组版本坐标转换

[复制链接]

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-2-20 23:43:59 | 显示全部楼层 |阅读模式
这个不是正规的基因组版本坐标转换需求,因为真正的基因组版本坐标转换一定是用软件,去我博客随便搜索一下就明白了。
我这里提出这个题目,是基于一个实际的需求!
我下载了23 and me 的芯片数据基因型文件,如下:
[mw_shl_code=applescript,true]# This data file generated by 23andMe at: Sat Feb 19 04:12:15 2011
#
# Below is a text version of your data. Fields are TAB-separated
# Each line corresponds to a single SNP.  For each SNP, we provide its identifier
# (an rsid or an internal id), its location on the reference human genome, and the
# genotype call oriented with respect to the plus strand on the human reference
# sequence.     We are using reference human assembly build 36.  Note that it is possible
# that data downloaded at different times may be different due to ongoing improvements
# in our ability to call genotypes. More information about these changes can be found at:
# https://www.23andme.com/you/download/revisions/
#
# More information on reference human assembly build 36:
# http://www.ncbi.nlm.nih.gov/proj ... taxid=9606&build=36
#
# rsid  chromosome      position        genotype
rs4477212       1       72017   AA
rs3094315       1       742429  AG
rs3131972       1       742584  AG
rs12124819      1       766409  --
rs11240777      1       788822  AG
[/mw_shl_code]

很明显,上面说了是基于 human assembly build 36坐标系的,也就是hg18,但是目前主流是hg19,所以需要转换坐标!
我这里希望大家去下载dbSNP数据库里面的ren的数据文件,然后基于hash这样的关联数组来做转换!
[mw_shl_code=shell,true]mkdir -p ~/annotation/variation/human/dbSNP
cd ~/annotation/variation/human/dbSNPwget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/All_20160601.vcf.gz
[/mw_shl_code]

我的脚本如下:
[mw_shl_code=applescript,true]zcat ~/annotation/variation/human/dbSNP/All_20160601.vcf.gz |perl -alne 'BEGIN{ open FH,"dm_23andme_v3_110219.txt";while(<FH>){chomp;@F=split;if(/^rs/){$pos{$.}=$_;$h{$F[0]}=$F[3]} } }{if(exists $h{$F[2]}){ $h{$F[2]}="$F[0]\t$F[1]\t$h{$F[2]}"  }}END{print "$pos{$_}\t$h{$pos{$_}}" foreach sort{$a<=>$b} keys %pos}' >dm_23andme_v3_hg19.txt
[/mw_shl_code]

你不需要看我的脚本,只需要实现需求即可!

这个dm_23andme_v3_110219.txt我以百度云方式给出!链接:https://pan.baidu.com/s/1hsLYOYS 密码:gsxd 生信技能树下载的23andme数据





上一篇:非模式生物转录组GO注释该怎么办?
下一篇:把vcf文件转为23andme芯片基因型文件
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-4-3 07:02 , Processed in 0.024767 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.