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

相关文章

来自专栏CDA数据分析师

深入对比数据科学工具箱:Python和R之争

概述 在真实的数据科学世界里,我们会有两个极端,一个是业务,一个是工程。偏向业务的数据科学被称为数据分析(Data Analysis),也就是A型数据科学。偏向...

1937
来自专栏编程一生

IO和socket编程

973
来自专栏程序人生

golang五周岁

时间过得真是太快,一晃眼golang都五周岁了。五年,对于人生来讲也许有些漫长,但对于一门语言的发展来说,仅仅是沧海一粟。相对于它试图挑战的「系统级」语言:c,...

3366
来自专栏PHP技术

8 个不得不说的 MySQL 陷阱

Mysql安装简单,速度较快,功能丰富。另外它还是开源运动的标杆,它的伟大成就向我们展示了一个成功的公司是可以建立在开源代码之上的。 然而用过mysql的人都曾...

3235
来自专栏坚毅的PHP

【读书笔记】linux编程艺术

一、哲学 1、Rule of Modularity: Write simple parts connected by clean interfaces. 1、模...

3219
来自专栏Python中文社区

Python=R+SQL/Hive?在数据分析与挖掘方面该选哪一个?

什么是R语言? R语言,一种自由软件编程语言与操作环境,主要用于统计分析、绘图、数据挖掘。R本来是由来自新西兰奥克兰大学的罗斯·伊哈卡和罗伯特·杰特曼开发(也因...

19310
来自专栏祝威廉

这些年,我工作上走过的路

我走过了毕业季,创业征途,踏进开源之路,转型进入大数据,到最后有缘接触机器学习。每个章节,我都会提及对应那个阶段对技术的感悟,自己做的一些具体事情。

522
来自专栏数据小魔方

大连市2016年空气质量数据可视化~

前几天发现了一个很有趣的包——openair,可以将年度时间序列刻画成周年日历热图,感觉这种形式非常适合用于呈现年度空气质量可视化,所以抓空爬了一些大连市201...

2766
来自专栏更流畅、简洁的软件开发方式

根据实例说数据库设计(一)—— 人员管理

     说在前面      可能您会问,树的系列还差第三篇没有写呢,怎么就又说数据库设计了?因为如果写第三篇的话,那么就涉及到了权限,而权限里面又涉及到了人员...

2487
来自专栏csxiaoyao

USTC高级软件工程课程学习心得

34710

扫描关注云+社区