gnomAD数据库简介（二）

在gnomAD数据库简介（一）中，我们简单介绍了基因组学遗传分析中人群变异频率的重要性，以及gnomAD数据库的一些背景。

本篇主要侧重gnomAD的后台数据下载和简单评估。

gnomAD后台数据下载

gnomAD数据下载的几个方式：

测试一下gsutil命令：

pip install gsutil

cd /home/shw/public/gnomAD

gsutil ls gs://gcp-public-data--gnomad/release/
gsutil ls gs://gcp-public-data--gnomad/release/2.1.1/liftover_grch38/vcf/exomes

为了简便一些，我们还是使用熟悉的wget命令下载：

nohup wget -c https://storage.googleapis.com/gcp-public-data--gnomad/release/2.1.1/liftover_grch38/vcf/exomes/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz &
head nohup.out

后台数据简单测试

查看上述获取的gnomAD（exomes, v2.1.1, LiftOver）VCF文件记录的变异位点个数：

zcat gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz | wc -l
# 17205543

gnomAD的这个外显子组数据共收录了约1,720万个变异位点！要知道人类总的外显子组位点数约为3,000万。这个比例依然很难得。随便找个基因的外显子序列，其中一半以上的核苷酸都能在gnomAD查到人群变异频率！

在该VCF文件中随机选择一个位点进行比较和测试，例如：rs1479269360

gnomAD后台数据（VCF文件的第5000行）

# 查看VCF文件的表头：
zcat gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz | head -n 5000 | grep -v '##' | head -n 1

# 查看VCF文件某一个变异位点的人群频率：
zcat gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz | head -n 5000 | tail -n 1

（人群变异频率）AF=7.44679e-06

另外注释有：转录本ID、密码子变化、（反式）调控位点注释等信息

gnomAD在线检索（AF完全匹配）

另有人群的亚群频率、年龄分布、基因型质量、测序深度、IGV等展示信息

dbSNP在线检索（发现居然没有该位点的AF）

另有临床意义等其它信息：

提取gnomAD的人群变异频率

从刚才的gnomAD（exomes, v2.1.1, LiftOver）VCF文件中提取AF信息：

nohup zcat gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz | sed 's/AF=/\t/g' | cut -f 9 | sed 's/;/\t/g' | cut -f 1 > gnomad.exomes.r2.1.1.sites.liftover_grch38.AF.txt
zcat gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz | cut -f 1-7 > gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.1-7col.txt &
# 按列合并：
paste  gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.1-7col.txt gnomad.exomes.r2.1.1.sites.liftover_grch38.AF.txt | grep -v '##' >  gnomad.exomes.r2.1.1.sites.liftover_grch38.7col_AF.txt

测试镰刀型贫血症的致病HBB的致病变异位点：rs334

grep -w rs334 gnomad.exomes.r2.1.1.sites.liftover_grch38.7col_AF.txt

chr11   5227002 rs334   T       A       2136270.15      PASS    3.47958e-03

完全匹配dbSNP网站上Frequency中的GnomAD_exome，且后者具有最大的人群基数：

https://www.ncbi.nlm.nih.gov/snp/rs334

使用gnomAD（v2.1.1）在线检索：

令人惊喜的是，gnomAD在线检索结果也提供了SIFT, Polyphen等in-sillico有害性预测，以及ClinVar相关注释信息：

关于ClinVar的详细介绍，及其对rs334注释，请查看：ClinVar数据库详解。

继续使用gnomAD（v3.1.1）在线检索：rs334（大小写敏感！）。结果中居然还有CADD和REVEL（In Silico Predictors）打分：

关于gnomAD的总的变异位点数

上述操作中，从gnomAD（exomes, v2.1.1, LiftOver）的VCF文件提取了AF（等位基因人群频率）信息，下面是其总的位点数：

wc -l   gnomad.exomes.r2.1.1.sites.liftover_grch38.7col_AF.txt
# 17,201,297 gnomad.exomes.r2.1.1.sites.liftover_grch38.7col_AF.txt

当然，我们更想了解所有3,000万个位点的变异频率。因为说不准哪天我们自己的外显子组测序数据就测到了一个导致氨基酸变异的位点，但恰好未被gnomAD收录（这种情况是存在的），此时由于不知道其AF，按照通常的思路只好考虑将其舍弃：只保留gnomAD中收录的、且AF<5%的位点。

那么gnomAD未收录的位点均被舍弃。也就是说，最终致病位点只能限制在gnomAD所收录的位点中（这依赖于gnomAD，是比较被动的）。此为“过分的舍弃”。

另一个思路，只过滤掉gnomAD中收录的、且AF>10%变异的位点，但保留下来的某些位点仍然可能在人群中存在高频变异（AF>10%），而这些位点有可能是耐受的、良性的或非致病的位点。此为“过多的保留”。

因此一些研究或高水平文献中不止参考了gnomAD，也参考了1000 Genomes和Bale database等数据库中收录的位点，目的就是尽量减少“过分的舍弃”和“过多的保留”。

因此我们还是希望gnomAD能覆盖到全部外显子序列（~3,000万个位点），这无疑是一个巨大挑战。

撰写：宋红卫

校对：宋红卫，叶明皓

更多人类遗传学知识、文献和分析技术

请关注和星标聊生信