||
http://www.zd200572.com/2017/10/24/snp_database_summary/
数据下载来自ncbi的ftp,想做这个的原因来自于jimmy的生信技能树的论坛帖子:
http://www.biotrainee.com/thread-265-1-1.html
我们都知道,一个人的全基因组测序会与参考基因组hg19/hg38有300~600万的snp差异,那么现在已有测序的所有样本与参考基因组hg19/hg38不同的地方是多少呢?全基因组测序已经超过10万例了,对snp位点记录最详细的当属 dbsnp啦,很容易去NCBI的ftp里面下载最新版的dbsnp文件,约15G大小,最新版的话!
如果只需要看统计信息,其实在dbsnp的官网主页即可看到。但是我们这个数据处理任务,所以必须去下载该文件,对文件进行解析。除了可以看有多少位点,还可以探索更多!欢迎爱学习的小伙伴在下面跟帖讨论自己的学习结果。
我下载的是下面这个地址的vcf.gz格式的,实在找不到txt格式的在哪,于是就下这个了。压缩文件大小是7个多G,解压之后有50多G,估计是包含了很多的注释信息。也没有选择对应于哪个参考基因组的snp,估计这个是最全的。对于这样的大文件,我习惯于IBM的ascp下载了(当然,wgetg一样可以),代码如下:
shell脚本应该也是可以搞定的,我对于shell是门外汉,先用python解决问题,尽管代码要多好几倍,估计shell也就一行的事。哪天学会了再补上。
处理思路是,发现文件里除了#开头注释行,其余均为数据,所以排除#开头的行,然后把染色体号,rs号,正常基因和变异基因提取到一个新的文件,代码如下:
可以看出,总数竟然高达325,174,796,这是三亿多呀,远远超出了我的想象,我以为只会有千万以内呢!但是,问题是现在能确定功能的snp位点只有上千个左右,也就是基因检测报告应用的数目,科技仍需加油呀!
用个字典解决问题吧!下面是是我的代码:
这个运行时间还算可以接受哈,500多秒,也就是十分钟不到。那么,排个序,哪个染色体最多,我猜应该和染色体大小一致。
可以看出2号染色体最多,那么看一下染色体的大小情况,来自论坛编程直播课程作业贴:
发现规律并不完全一致。好了,探索先到这里。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-20 04:31
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社