人群频率 | gnomAD数据库 (二) 后台数据的获取及质量评估

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

本文主要侧重于gnomAD下载和简单评估背景数据。

gnomAD下载后台数据

gnomAD下载数据的几种方式：

测试一下gsutil命令：

pip install gsutil   cd /home/shw/public/gnomAD   gsutillsgs://gcp-public-data--gnomad/release/ gsutillsgs://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 & headnohup.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信息：

nohupzcatgnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz|sed's/AF=/\t/g'|cut-f9|sed's/;/\t/g'|cut-f1>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-wrs334gnomad.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详细介绍，及其对rs请查看334注释：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

当然，我们想知道所有3000万个位点的变异频率。因为不确定哪一天我们自己的外显子组测序数据测量了一个导致氨基酸变异的位点，但恰好没有gnomAD收录(这种情况是存在的)，此时不知道AF，按照通常的思路只好考虑放弃：只保留gnomAD中收录的、且AF<5%的位点。

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

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

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

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

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

请关注和星标聊生信