用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 条评论
登录 后参与评论

相关文章

来自专栏CSDN技术头条

你被复制了吗?大数据带你看中国人的名字 重名TOP榜

名字是一个人的代号,起名字也是一门学问。你想知道中国人重名最多的名字是哪个吗?你想知道中国人重名最多的名都有哪些吗?再看看身边“10后”的娃们,有没有叫子涵、子...

1878
来自专栏生信宝典

UCSC XENA - 集大成者(TCGA, ICGC)

TCGA有自己的一批工具,ICGC也有自己的网站,但好的资源都是要整合起来,整合越多越好(虽然事实不一定如此,但有这个想法的人不少),用着才更方便。这就靠今天介...

6333
来自专栏Alan's Lab

windows 10 私人不负责任评测(多图预警)

最牛的一个: Mac 下 Flash 一直是个痛,只有用虚拟机的 Win 7 才勉强敢看,但拿 Win 10 试了一盘三国杀+两部腾讯视频(拿 Mac 看过腾讯...

863
来自专栏Python爬虫与算法进阶

大佬,我代码哪错了?

问题无处不在 我: “大佬,帮我看看这个问题错在哪了?” 大佬: “你的代码呢、你的错误提示呢?” 我: “好的,我发给你” 大佬: “。。。 再见” 留下一脸...

3488
来自专栏FreeBuf

网络蜜罐技术探讨

*原创作者bt0sea,本文属FreeBuf原创奖励计划,未经许可禁止转载 之前在FreeBuf发表的第一篇有关蜜罐文章,引起了业界不小的轰动,但是...

3789
来自专栏生信宝典

生信宝典之傻瓜式 (二) 如何快速查找指定基因的调控网络

我是谁?我在哪儿?我在查什么? 在信息爆炸的时代,相信很多小伙伴在查文章时会因信息量太大而抓狂。今天带来一款设计简洁、功能全面的基因功能查询工具,助你事半功倍,...

2106
来自专栏生信技能树

【直播】我的基因组63:wegene芯片跟二代测序的简单比较

首先要说的是,wegene公司效率很赞。 从我送样品(2017-2-20)到拿到测序数据(2017-3-10),总共耗时不到20天~~~ ? 前面我们已经介绍过...

3748
来自专栏Golang语言社区

【Golang语言社区】--无线视频监控传输原理

监视前端一体摄像机安装在圆形防护罩内,内带全方位云台。摄制的图像转换成视频信号传输到微波发射机的调制端,微波发射机将其加载到载波上,经微波天线定向辐射到监控中心...

4599
来自专栏FreeBuf

谷歌网络服务宕机,中国电信背锅

上周日,谷歌旗下的云服务、YouTube等网络服务在全球范围内均发生了数小时的宕机,外媒称因遭到来自中国电信IP的BGP劫持导致故障发生。虽然这次事件为中国电信...

1022
来自专栏SDNLAB

SDN实战团分享(三十一):Segment Routing meet SDN

一、介绍 ? 在1990年代Yakov, Eric Rosen, Kompella很多业界先驱(仅列举了Juniper公司的MPLS业界领袖,其他公司也有 很多...

63616

扫码关注云+社区