搜索
查看: 2391|回复: 0

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

[复制链接]

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-2-20 23:43:59 | 显示全部楼层 |阅读模式
这个不是正规的基因组版本坐标转换需求,因为真正的基因组版本坐标转换一定是用软件,去我博客随便搜索一下就明白了。
我这里提出这个题目,是基于一个实际的需求!
我下载了23 and me 的芯片数据基因型文件,如下:
[AppleScript] 纯文本查看 复制代码
# 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:
# [url=https://www.23andme.com/you/download/revisions/]https://www.23andme.com/you/download/revisions/[/url]
# 
# More information on reference human assembly build 36:
# [url=http://www.ncbi.nlm.nih.gov/projects/mapview/map_search.cgi?taxid=9606&build=36]http://www.ncbi.nlm.nih.gov/proj ... taxid=9606&build=36[/url]
#
# rsid  chromosome      position        genotype
rs4477212       1       72017   AA
rs3094315       1       742429  AG
rs3131972       1       742584  AG
rs12124819      1       766409  --
rs11240777      1       788822  AG


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


我的脚本如下:
[AppleScript] 纯文本查看 复制代码
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 


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

这个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, 2019-11-21 04:59 , Processed in 0.026304 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.