add gnomad exome allele frequencies.
MODIFY
docs/content/database_schema.rst
docs/content/database_schema.rst
@@ -169,6 +169,20 @@ max_aaf_all FLOAT the maximum of aaf_esp_ea, aaf_esp_aa, a
exac_num_het INTEGER The number of heterozygote genotypes observed in ExAC. Pulled from the ExAC ``AC_Het`` INFO field.
exac_num_hom_alt INTEGER The number of homozygous alt. genotypes observed in ExAC. Pulled from the ExAC ``AC_Het`` INFO field.
exac_num_chroms INTEGER The number of chromosomes underlying the ExAC variant call. Pulled from the ExAC ``AN_Adj`` INFO field.
+
+aaf_gnomad_all FLOAT Allele frequency (population independent) of the variant in gnomad,
+aaf_gnomad_afr FLOAT Allele frequency (AFR population) of the variant in gnomad
+aaf_gnomad_amr FLOAT Allele frequency (AMR population) of the variant in gnomad
+aaf_gnomad_asj FLOAT Allele frequency (ASJ population) of the variant in gnomad
+aaf_gnomad_eas FLOAT Allele frequency (EAS population) of the variant in gnomad
+aaf_gnomad_fin FLOAT Allele frequency (FIN population) of the variant in gnomad
+aaf_gnomad_nfe FLOAT Allele frequency (NFE population) of the variant in gnomad
+aaf_gnomad_oth FLOAT Allele frequency (OTH population) of the variant in gnomad
+aaf_gnomad_sas FLOAT Allele frequency (SAS population) of the variant in gnomad
+gnomad_num_het INTEGER Number of het genotypes observed in gnomad
+gnomad_num_hom_alt INTEGER Number of hom_alt genotypes observed in gnomad
+gnomad_num_chroms INTEGER Number of chromosomes genotyped in gnomad
+
======================== ========= =================================================================================================

MODIFY
docs/content/history.rst
docs/content/history.rst
@@ -8,6 +8,7 @@ Release History
0.20.0
======
+#. Include gnomad exome allele frequencies for all populations and include these in calc of max_aaf_all.
#. Initial support for python 3.
#. Support annotating with bcftools BCSQ (http://biorxiv.org/content/early/2016/12/01/090811) in addition to VEP and SnpEff.
#. Fix bug that resulted in the message: `UnsupportedOperation: IOStream has no fileno.` when loading. (Thanks Andrew O. and Ryan R.)

ADD
gemini/annotation_provenance/make-gnomad.sh
+wget -O - https://data.broadinstitute.org/gnomAD/release-170228/exomes/vcf/gnomad.exomes.r2.0.1.sites.vcf.gz | zcat - | perl -pe 's/;CSQ.+//g' | less | bgzip -c > gnomad.exomes.r2.0.1.sites.no-VEP.vcf.gz
+
+zcat gnomad.exomes.r2.0.1.sites.no-VEP.vcf.gz \
+ | perl -pe 's/(;[^_\s]+_HIST(_AL[LT])?=[^;]+)//g' \
+ | vt decompose -s - \
+ | vt normalize -r ~/bcbio/genomes/Hsapiens/g1k_v37_decoy/seq/g1k_v37_decoy.fa - \
+ | bgzip -c > gnomad.exomes.r2.0.1.sites.no-VEP.nohist.tidy.vcf.gz
+
+tabix -f gnomad.exomes.r2.0.1.sites.no-VEP.nohist.vcf.gz

MODIFY
gemini/annotations.py
gemini/annotations.py
@@ -53,6 +53,7 @@ def get_anno_files(args):
'fitcons': os.path.join(anno_dirname, "hg19_fitcons_fc-i6-0_V1-01.bed.gz"),
'cosmic': os.path.join(anno_dirname, 'cosmic-v68-GRCh37.tidy.vcf.gz'),
'exac': os.path.join(anno_dirname, 'ExAC.r0.3.sites.vep.tidy.vcf.gz'),
+ 'gnomad': os.path.join(anno_dirname, 'gnomad.exomes.r2.0.1.sites.no-VEP.nohist.tidy.vcf.gz'),
'geno2mp': os.path.join(anno_dirname, 'geno2mp.variants.tidy.vcf.gz'),
}
# optional annotations
@@ -175,6 +176,23 @@ ExacInfo = collections.namedtuple("ExacInfo",
EXAC_EMPTY = ExacInfo(False, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1)
+GnomadInfo = collections.namedtuple('GnomadInfo',
+ "aaf_ALL \
+ aaf_AFR \
+ aaf_AMR \
+ aaf_ASJ \
+ aaf_EAS \
+ aaf_FIN \
+ aaf_NFE \
+ aaf_OTH \
+ aaf_SAS \
+ num_het \
+ num_hom_alt \
+ num_chroms")
+
+GNOMAD_EMPTY = GnomadInfo(-1, -1, -1, -1, -1, -1,
+ -1, -1, -1, -1, -1, -1)
+
def load_annos(args):
"""
Populate a dictionary of Tabixfile handles for
@@ -776,6 +794,68 @@ def get_geno2mp_ct(var):
# missing is -1
return -1
+def get_gnomad_info(var, empty=GNOMAD_EMPTY):
+ info_map = {}
+ afs = {}
+ for hit in annotations_in_vcf(var, "gnomad", "vcf", "grch37"):
+ # Does not handle anything beyond var.ALT[0] in the VCF (in case of multi-allelic variants)
+ # var.start is used since the chromosomal pos in pysam.asVCF is zero based (hit.pos)
+ # and would be equivalent to (POS-1) i.e var.start
+ if var.start != hit.pos or var.REF != hit.ref:
+ continue
+
+ # This would look for var.ALT[0] matches to
+ # any of the multiple alt alleles represented in the EXAC file
+ ALT = hit.alt.split(",")
+ for allele_num, each in enumerate(ALT):
+ if each != var.ALT[0]:
+ continue
+
+ # Store the allele index of the match to retrieve the right frequencies
+ for info in hit.info.split(";"):
+ if "=" in info:
+ (key, value) = info.split("=", 1)
+ info_map[key] = value
+
+ # Population independent raw (non-adjusted) allele frequencies given by AF
+ if info_map.get('AF', '.') != '.':
+ aaf_ALL = float(info_map['AF'].split(",")[allele_num])
+ else:
+ aaf_ALL = -1
+
+ for grp in ('_AFR', '_AMR', '_ASJ', '_EAS', '_FIN', '_NFE', '_OTH', '_SAS'):
+ ac = info_map.get('AC' + grp)
+ if ac is None: continue
+
+ an = info_map.get('AN' + grp)
+ if an is None: continue
+
+ if an == '0':
+ afs[grp] = 0
+ continue
+
+ ac_list = ac.split(",")
+ afs[grp] = float(ac_list[allele_num]) / float(an)
+
+ nhm = sum(map(int, info_map.get('GC_Male', '').split(",")[1:-1]))
+ nhf = sum(map(int, info_map.get('GC_Female', '').split(",")[1:-1]))
+
+ num_hets = nhm + nhf
+ num_homs = int(info_map.get("Hom", -1))
+
+ called_chroms = int(info_map.get('AN', -1))
+
+ return GnomadInfo(aaf_ALL, float(afs['_AFR']),
+ float(afs['_AMR']), float(afs['_ASJ']),
+ float(afs['_EAS']), float(afs['_FIN']),
+ float(afs['_NFE']), float(afs['_OTH']),
+ float(afs['_SAS']), num_hets, num_homs,
+ called_chroms)
+
+ return empty
+
+
+
def get_exac_info(var, empty=EXAC_EMPTY):
"""
Returns the allele frequencies from the Exac data (Broad)

MODIFY
gemini/database.py
gemini/database.py
@@ -259,6 +259,20 @@ def create_tables(path, effect_fields=None):
exac_num_het int,
exac_num_hom_alt int,
exac_num_chroms int,
+
+ aaf_gnomad_all decimal(2,7),
+ aaf_gnomad_afr decimal(2,7),
+ aaf_gnomad_amr decimal(2,7),
+ aaf_gnomad_asj decimal(2,7),
+ aaf_gnomad_eas decimal(2,7),
+ aaf_gnomad_fin decimal(2,7),
+ aaf_gnomad_nfe decimal(2,7),
+ aaf_gnomad_oth decimal(2,7),
+ aaf_gnomad_sas decimal(2,7),
+
+ gnomad_num_het int,
+ gnomad_num_hom_alt int,
+ gnomad_num_chroms int,
%s""" % effect_string.rstrip(","),
variant_impacts="""

MODIFY
gemini/gemini_load_chunk.py
gemini/gemini_load_chunk.py
@@ -422,6 +422,7 @@ class GeminiLoader(object):
cosmic_ids = annotations.get_cosmic_info(var)
fitcons = annotations.get_fitcons(var)
Exac = annotations.get_exac_info(var)
+ gnomad = annotations.get_gnomad_info(var)
#load CADD scores by default
if self.args.skip_cadd is False:
@@ -447,6 +448,7 @@ class GeminiLoader(object):
esp = annotations.ESPInfo(False, -1, -1, -1, 0)
thousandG = annotations.EMPTY_1000G
Exac = annotations.EXAC_EMPTY
+ gnomad = annotations.GNOMAD_EMPTY
recomb_rate = None
gms = annotations.GmsTechs(None, None, None)
grc = None
@@ -704,7 +706,24 @@ class GeminiLoader(object):
aaf_adj_exac_sas=Exac.aaf_SAS,
exac_num_het=Exac.num_het,
exac_num_hom_alt=Exac.num_hom_alt,
- exac_num_chroms=Exac.num_chroms)
+ exac_num_chroms=Exac.num_chroms,
+
+ aaf_gnomad_all=gnomad.aaf_ALL,
+ aaf_gnomad_afr=gnomad.aaf_AFR,
+ aaf_gnomad_amr=gnomad.aaf_AMR,
+ aaf_gnomad_asj=gnomad.aaf_ASJ,
+ aaf_gnomad_eas=gnomad.aaf_EAS,
+ aaf_gnomad_fin=gnomad.aaf_FIN,
+ aaf_gnomad_nfe=gnomad.aaf_NFE,
+ aaf_gnomad_oth=gnomad.aaf_OTH,
+ aaf_gnomad_sas=gnomad.aaf_SAS,
+
+ gnomad_num_het=gnomad.num_het,
+ gnomad_num_hom_alt=gnomad.num_hom_alt,
+ gnomad_num_chroms=gnomad.num_chroms,
+
+
+ )
variant['max_aaf_all'] = max(-1,
variant['aaf_esp_ea'],

MODIFY
gemini/install-data.py
gemini/install-data.py
@@ -52,6 +52,7 @@ anno_files = \
'summary_gene_table_v75',
'cancer_gene_census.20140120.tsv',
'ExAC.r0.3.sites.vep.tidy.vcf.gz',
+'gnomad.exomes.r2.0.1.sites.no-VEP.nohist.tidy.vcf.gz',
'geno2mp.variants.tidy.vcf.gz',
]
extra_anno_files = {"gerp_bp": "hg19.gerp.bw", "cadd_score": "whole_genome_SNVs.tsv.compressed.gz"}

MODIFY
master-test.sh
master-test.sh
@@ -70,6 +70,7 @@ bash test-clinvar.sh
# Test Exac
bash test-exac.sh
+bash test-gnomad.sh
# Test population_gen metrics
bash test-pop.sh

ADD
test/test-gnomad.sh
+#############################################
+# 1. Test ExAC allele frequencies
+#############################################
+check()
+{
+ if diff $1 $2; then
+ echo ok
+ else
+ echo fail
+ fi
+}
+export -f check
+
+
+echo " gnomad.t01...\c"
+echo "chrom start end ref alt aaf_gnomad_all aaf_gnomad_all aaf_gnomad_afr aaf_gnomad_amr aaf_gnomad_eas aaf_gnomad_fin aaf_gnomad_nfe aaf_gnomad_oth aaf_gnomad_sas
+chr1 985954 985955 G C -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0
+chr1 1199488 1199489 G A -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0
+chr1 1959698 1959699 G A 0.0174018 0.0174018 0.0036564286664 0.00955851604038 0.000116387337058 0.0382032667877 0.0233466915698 0.0213392200147 0.00599856556041
+chr1 161276552 161276553 G T -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0
+chr1 247587092 247587093 C T -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0
+chr9 112185055 112185056 C G 8.12315e-06 8.12315e-06 0.000130701869037 0.0 0.0 0.0 0.0 0.0 0.0
+chr12 121432116 121432118 GC G 6.27201e-05 6.27201e-05 0.000154202004626 3.19570497252e-05 6.25312656328e-05 5.06790999392e-05 8.0064051241e-05 0.000198098256735 0.0
+chr14 105420589 105420590 C T 0.00166166 0.00166166 0.000261574679571 0.000863609291245 0.0 0.000183941874368 0.00187292768169 0.00164594001463 0.00162453700695" > exp.$$
+gemini query --header -q "select chrom, start, end, ref, alt, aaf_gnomad_all, aaf_gnomad_all, \
+ aaf_gnomad_afr, aaf_gnomad_amr, aaf_gnomad_eas, aaf_gnomad_fin, \
+ aaf_gnomad_nfe, aaf_gnomad_oth, aaf_gnomad_sas from variants" test.exac.db > obs.$$
+
+
+check obs.$$ exp.$$
+#rm obs.$$ exp.$$
+
+#############################################

add gnomad exome allele frequencies.