Add assemble_constraint_context_ht function to create a fully annotated context HT for computing constraint on

[jkgoodrich]

Also adds some functions to help transform the methylation data for annotation onto the context HT:

• transform_methylation_level
• transform_grch37_methylation
• transform_grch38_methylation

Tested with:

from gnomad.resources.grch38.gnomad import coverage, all_sites_an, public_release
from gnomad.resources.grch38.reference_data import methylation_sites, vep_context

context_ht = vep_context.ht()
coverage_hts = {
    "exomes": coverage("exomes").ht(),
    "genomes": coverage("genomes").ht(),
}
an_hts = {
    "exomes": all_sites_an("exomes").ht(),
    "genomes": all_sites_an("genomes").ht(),
}
exome_ht = public_release("exomes").ht()
genomes_ht = public_release("genomes").ht()
freq_hts = {
    "exomes": exome_ht.select("freq"),
    "genomes": genomes_ht.select("freq"),
}
filter_hts = {
    "exomes": exome_ht.select("filters"),
    "genomes": genomes_ht.select("filters"),
}
methylation_ht = methylation_sites.ht()
gerp_ht = hl.experimental.load_dataset(name="gerp_scores", version="hg19", reference_genome="GRCh38")

annotated_context_ht = assemble_constraint_context_ht(
    context_ht,
    coverage_hts=coverage_hts,
    an_hts=an_hts,
    freq_hts=freq_hts,
    filter_hts=filter_hts,
    methylation_ht=methylation_ht,
    gerp_ht=gerp_ht,
    transformation_funcs=None,
)

annotated_context_ht.describe()
annotated_context_ht.show()

and

from gnomad.resources.grch37.gnomad import coverage, public_release
from gnomad.resources.grch37.reference_data import methylation_sites, vep_context

context_ht = vep_context.ht()
coverage_hts = {
    "exomes": coverage("exomes").ht(),
    "genomes": coverage("genomes").ht(),
}
an_hts = None
exome_ht = public_release("exomes").ht()
genomes_ht = public_release("genomes").ht()
freq_hts = {
    "exomes": exome_ht.select("freq"),
    "genomes": genomes_ht.select("freq"),
}
filter_hts = {
    "exomes": exome_ht.select("filters"),
    "genomes": genomes_ht.select("filters"),
}
methylation_ht = methylation_sites.ht()
gerp_ht = hl.experimental.load_dataset(name="gerp_scores", version="hg19", reference_genome="GRCh37")

annotated_context_ht = assemble_constraint_context_ht(
    context_ht,
    coverage_hts=coverage_hts,
    an_hts=an_hts,
    freq_hts=freq_hts,
    filter_hts=filter_hts,
    methylation_ht=methylation_ht,
    gerp_ht=gerp_ht,
    transformation_funcs=None,
)

annotated_context_ht.describe()
annotated_context_ht.show()

[KoalaQin]

I'm sending the first round of my review back, I think it's very clear to combine these tables to context HT with dictionaries, I have a few questions to help me imagine the abstract part.

[KoalaQin]

Why do you need 'most_severe_consequence' under both 'vep' and 'transcript_consequences'?

[jkgoodrich]

transcript_consequences is a list, so the most_severe_consequence in each element of that list is the most severe consequence for the specific transcript. The higher level most_severe_consequence is the most severe consequence across all transcripts at the variant

[KoalaQin]

oh, I didn't know that, I thought one variant would only have one consequence in one transcript, it would also choose between "missense" and "splice region variant" to be most severe for one transcript.

[jkgoodrich]

Yeah, it's whatever the most significant consequence is (based on csq_order) in the consequence_terms annotation of the element in transcript_consequences.

[KoalaQin]

I don't see coverage_stats in the coverage_hts in the ones you tested, which coverage_hts are you referring to?

[jkgoodrich]

The comment specifies that this is relevant to v4...

It's added here:

gnomad_methods/gnomad/utils/sparse_mt.py

Line 1334 in d29acb3

"coverage_stats": (

[KoalaQin]

I meant when I tried to get coverage('exomes').ht() for v4, I don't see coverage_stats as you mentioned in the comment, but the nested rows are there. Are there other versions of coverage table?

[KoalaQin]

I see there's a difference between GRCh37 and GRCh38 coverage table, the columns changed from:

Row fields:
    'row_id': int64 
    'locus': locus<GRCh37> 
    'mean': float64 
    'median': int32 

to:

Row fields:
    'locus': locus<GRCh38> 
    'mean': float64 
    'median_approx': int32 
    'total_DP': int64 

[jkgoodrich]

It's in this one: "gs://gnomad/release/4.0/ht/exomes/gnomad.exomes.v4.0.coverage.ht"

from gnomad_qc.v4.resources.release import release_coverage
ht = release_coverage("exomes", public=False).ht()
ht.describe()

----------------------------------------
Global fields:
    'coverage_stats_meta': array<dict<str, str>> 
    'coverage_stats_meta_sample_count': array<int32> 
----------------------------------------
Row fields:
    'locus': locus<GRCh38> 
    'coverage_stats': array<struct {
        mean: float64, 
        median_approx: int32, 
        total_DP: int64, 
        over_1: float64, 
        over_5: float64, 
        over_10: float64, 
        over_15: float64, 
        over_20: float64, 
        over_25: float64, 
        over_30: float64, 
        over_50: float64, 
        over_100: float64
    }> 
----------------------------------------
Key: ['locus']
----------------------------------------

The point is that it can work with either possible schema

[KoalaQin]

Okay, that's what I want to know, thanks.

[KoalaQin]

A few comments, we're close.
Sorry, it's been a while and I forgot.

[KoalaQin]

oh, I didn't know that, I thought one variant would only have one consequence in one transcript, it would also choose between "missense" and "splice region variant" to be most severe for one transcript.

[KoalaQin]

I meant when I tried to get coverage('exomes').ht() for v4, I don't see coverage_stats as you mentioned in the comment, but the nested rows are there. Are there other versions of coverage table?

[KoalaQin]

LGTM!