用gnomDB数据库对个人vcf变异文件进行过滤

直播我的基因组前面的上游分析到此为止了,这里是一个分界线,经过孜孜不倦的探索挖掘我已经拿到了我个人基因组跟hg19参考基因组的全部差异位点,而且可以肯定方法学上面没有毛病。现在到了解释这些差异位点的时候,或者说是注释它们。

 754755 indel.vcf3784343 snp.vcf

三百多万的snp和近100万的indel仍然是天文数字,前面我多次强调人类的hg19参考基因组并不意味着都是好的,我的DNA跟参考基因组不一样反而是好事,而且更多的位点,仅仅是多态性而已,那么我们就应该在数据分析的过程中把位点区分开来。

首先,来一个最简单的,过滤掉人群突变位点,做这个分析是基于一个显而易见的假设,如果人群中有不少人都是在某个位点跟参考基因组不一样,那么这个位点,至少不是致命的,一般来说也不会是有害的。而公共人群数据库比较出名的有,1000基因组数据库,NHLBI外显子测序数据库,EXAC数据库,gnomAD数据库等。目前 gnomAD数据库是最大最全,而且最新的一个,我们就直接用它吧。

gnomAD数据库背景介绍

GenomeAggregation Database(简称gnomAD)是由各国研究者联合发展起来的基因组突变频率数据库。其目的是汇集和协调来自众多大规模测序计划的全外显子组和全基因组测序数据,为广泛的科学研究团体汇总数据。

该数据库提供的数据集包括123,136个个体的全外显子组测序数据和15,496个个体的全基因组测序数据,这些数据来源于各种疾病研究项目及大型人群测序项目。

该数据库所有的数据都可免费下载。

下载最方便的就是 google的gsutil啦,但是墙内的朋友有点麻烦,而且数据量也的确是太大了。

gsutil -m cp -r gs://gnomad-public/release/2.0.2/vds/exomes/gnomad.exomes.r2.0.2.sites.vds gnomad_data # 16 GB          gsutil -m cp -r gs://gnomad-public/release/2.0.2/vds/exomes/gnomad.exomes.r2.0.2.sites.split.vds gnomad_datagsutil -m cp -r gs://gnomad-public/release/2.0.2/vds/genomes/gnomad.genomes.r2.0.2.sites.vds gnomad_data # 108 GB

如果我们本身只需要该数据库的人群频率信息,其实没必要下载全部的vcf文件, 这里调用 annovar 软件整理好的数据库吧:

nohup /public/biosoft/ANNOVAR/annovar/annotate_variation.pl \-downdb -webfrom annovar gnomad_genome  \--buildver hg19 /public/biosoft/ANNOVAR/annovar/humandb/ >down.log 2>&1 &## 16G Mar 12  2017 /public/biosoft/ANNOVAR/annovar/humandb/hg19_gnomad_genome.txt

仍然是有 16 G,唉,人类遗传研究不容易啊, 简单查看文件内容如下:

#Chr    Start   End Ref Alt gnomAD_genome_ALL   gnomAD_genome_AFR   gnomAD_genome_AMR   gnomAD_genome_ASJ   gnomAD_genome_EAS   gnomAD_genome_FIN   gnomAD_genome_NFE   gnomAD_genome_OTH1       13538   13538   G       A       3.378e-05       0       0       0       0       0       7.225e-05       01       13540   13540   T       C       0       0       0       0       0       0       0       01       10327   10327   T       C       0.256   0.2609  0.5     0.25    0.5     0.2759  0.1981  0.2143

虽然不是vcf格式了,但是该有的信息都还在,很容易去gnomAD数据库查询情况,比如:http://gnomad.broadinstitute.org/variant/12-121437382-A-G 相信正常人都可以看出这样的url是有规律的,自己感兴趣的变异位点,可以链接到网站里面查看下详细的信息。

http://gnomad.broadinstitute.org/variant/12-121435427-G-Ahttp://gnomad.broadinstitute.org/variant/12-121437382-A-Ghttp://gnomad.broadinstitute.org/variant/1-13569-C-T

比较重要的信息,就是变异的基因组坐标以及其在不同人群发生的频率咯:

#ChrStartEndRefAltgnomAD_genome_ALLgnomAD_genome_AFRgnomAD_genome_AMRgnomAD_genome_ASJgnomAD_genome_EASgnomAD_genome_FINgnomAD_genome_NFEgnomAD_genome_OTH

人群的全称是:ALL, AFR (African), AMR (Admixed American), EAS (East Asian), FIN (Finnish), NFE (Non-finnish European), OTH (other), SAS (South Asian).

这里值得一提的是,ANNOVAR这个软件提供的 hg19_gnomad_genome.txt文件,有3亿行,意味着人类几乎10%的位点都被囊括了,而大家看到上面截取的文件内容里面有很多位点,在任何人群里面的发生频率都是0,理论上这样的位点是不需要列出的。在gnomAD数据库里面有12,288,392个位点都是人群频率大于5%,有 281,634,375是小于5%的。 也就是说人群频率大于5%的位点是少数派,人类这个整体,差异没有我们想象的那么大。

根据人群频率来进行过滤

/public/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old  snp.vcf  >snp_input/public/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old  indel.vcf  >indel_input/public/biosoft/ANNOVAR/annovar/annotate_variation.pl \-filter -dbtype gnomad_genome -buildver hg19 -out snp_filter snp_input \/public/biosoft/ANNOVAR/annovar/humandb/ -score_threshold  0.05/public/biosoft/ANNOVAR/annovar/annotate_variation.pl \-filter -dbtype gnomad_genome -buildver hg19 -out indel_filter indel_input \/public/biosoft/ANNOVAR/annovar/humandb/ -score_threshold  0.05

这种需要进行格式转换的软件我其实不太喜欢用,把标准的vcf文件给转换了,到时候其它下游分析,可能还得转回来,太麻烦了。还是简单给大家看看日志吧,这个也很重要:

NOTICE: Read 3784343 lines and wrote 3784225 different variants at 3784225 genomic positions (3784225 SNPs and 0 indels)NOTICE: Among 3784225 different variants at 3784225 positions, 2200297 are heterozygotes, 1583928 are homozygotesNOTICE: Among 3784225 SNPs, 2535708 are transitions, 1246859 are transversions (ratio=2.03), 1658 have more than 2 allelesNOTICE: Read 754755 lines and wrote 723138 different variants at 754637 genomic positions (0 SNPs and 754637 indels)NOTICE: Among 754637 different variants at 754637 positions, 410912 are heterozygotes, 312226 are homozygotesNOTICE: Among 0 SNPs, 0 are transitions, 0 are transversions (ratio=NA)

对3784343个的SNP位点来说,3353921个因为人群频率大于了0.05会被过滤掉,还剩下430304值得我好好研究一下。

但是,430304个变异位点还是有点多啊!!!!

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-12-31

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏机器人网

史上最全汽车图解,老司机也未必认得完 (收藏版)

这份史上最全汽车各部件图解,非常值得收藏!就算是老师傅,有很多部件的名字你肯定听说过但不一定都知道在哪个位置吧。 ? 打开发动机盖,就是这个样子了,4A13...

2758
来自专栏FreeBuf

如何自己动手制作一个靠谱的PM2.5检测仪

气象专家和医学专家认为,由细颗粒物造成的灰霾天气对人体健康的危害甚至要比沙尘暴更大。粒径10微米以上的颗粒物,会被挡在人的鼻子外面;粒径在2.5微米至10微米之...

3408
来自专栏用户2442861的专栏

Pyqt phonon的使用

Qt phonon地址:http://wenku.baidu.com/link?url=nH_dZ8lZbXHy8N5__8jAWLXcuMYf4yRjdCK...

401
来自专栏大数据文摘

量子生命|生物信息是如何编码的, 又是如何遗传的?

1021
来自专栏吉浦迅科技

辣眼睛,看别人如何用Jetson TX1花式逗猫!

如果你最近在阅读人工智能(artificial intelligence,简称AI)、深度学习(deep learning)和神经网络(neural netwo...

2705
来自专栏自然语言处理

微博话题爬取与存储分析(下)

本文基于python以新浪微博为数据平台,从数据采集、关键字提取、数据存储三个角度,用最简单的策略来挖掘我们的“黄金”。

752
来自专栏FreeBuf

RSA大会播报 – 2014最佳安全博客提名(国外篇)

今年美国RSA大会将在这个月的23-28号举行,每年大会上都会评出过去一年来业内最佳安全博客(Security Bloggers Network Social ...

1545
来自专栏大数据文摘

数学之美:图论和网络爬虫

2544
来自专栏生信技能树

TCGA的28篇教程- 使用R语言的RTCGAToolbox包获取TCGA数据

前些天被TCGA的终结新闻刷屏,但是一直比较忙,还没来得及仔细研读,但是笔记本躺着的一些TCGA教程快发霉了,借此契机好好整理一下吧,预计28篇教程!

1833
来自专栏CDA数据分析师

骑行在华盛顿 针对320万次共享单车骑行数据的分析

原作者 Alex Freeman 编译 Mika 本文为 CDA 数据分析师原创作品,转载需授权 共享单车的蓬勃发展,让人们的短途出行更加方便,这种绿色低...

21110

扫码关注云+社区