用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个变异位点还是有点多啊！！！！