gnomad.utils.constraint

Script containing generic constraint functions that may be used in the constraint pipeline.

`gnomad.utils.constraint.COVERAGE_CUTOFF`	Minimum median exome coverage differentiating high coverage sites from low coverage sites.
`gnomad.utils.constraint.CLASSIC_LOF_ANNOTATIONS`	Classic loss-of-function VEP annotations.
`gnomad.utils.constraint.DEFAULT_FIELDS_TO_SUM`	Default fields summed when aggregating expected variants by constraint group.
`gnomad.utils.constraint.DEFAULT_GENCODE_ANNOTATIONS`	Default GENCODE annotations added to a Table by transcript id.
`gnomad.utils.constraint.get_mu_annotation_expr`(ht, ...)	Get mutation rate annotation expression from the mutation rate Table.
`gnomad.utils.constraint.annotate_with_mu`(ht, ...)	Annotate SNP mutation rate for the input Table.
`gnomad.utils.constraint.variant_observed_expr`([...])	Return 1 if a variant meets frequency criteria, 0 otherwise.
`gnomad.utils.constraint.variant_observed_and_possible_expr`(...)	Return per-variant observed and possible count expressions.
`gnomad.utils.constraint.weighted_sum_agg_expr`(...)	Return the weighted aggregate sum of `expr` weighted by `weight_expr`.
`gnomad.utils.constraint.counts_agg_expr`([...])	Return an aggregation expression for variant and singleton counts.
`gnomad.utils.constraint.build_constraint_consequence_groups`(...)	Build constraint consequence groups.
`gnomad.utils.constraint.count_observed_and_possible_by_group`(ht, ...)	Aggregate observed and possible variant counts by substitution context group.
`gnomad.utils.constraint.explode_downsamplings_oe`(ht, ...)	Explode observed and expected downsampling counts for each genetic ancestry group and metric.
`gnomad.utils.constraint.annotate_mutation_type`(t)	Annotate mutation types.
`gnomad.utils.constraint.trimer_from_heptamer`(t)	Trim heptamer context to create trimer context.
`gnomad.utils.constraint.collapse_strand`(t)	Return the deduplicated context by collapsing DNA strands.
`gnomad.utils.constraint.transform_methylation_level`([...])	Transform methylation level into a 0 to 2 scale.
`gnomad.utils.constraint.assemble_constraint_context_ht`(ht)	Assemble context Table with necessary annotations for constraint calculations.
`gnomad.utils.constraint.calculate_gerp_cutoffs`(ht)	Find GERP cutoffs at the given percentile thresholds.
`gnomad.utils.constraint.calibration_model_group_expr`(...)	Get the calibration model grouping annotation for a variant.
`gnomad.utils.constraint.build_models`(ht, ...)	Build coverage and plateau models.
`gnomad.utils.constraint.build_plateau_models`(...)	Build plateau models to calibrate mutation rate against proportion observed.
`gnomad.utils.constraint.build_coverage_model`(...)	Build coverage model.
`gnomad.utils.constraint.get_all_gen_anc_lengths`(ht, ...)	Get the minimum length of observed variant counts array for each genetic ancestry group downsampling.
`gnomad.utils.constraint.get_constraint_grouping_expr`(...)	Collect annotations used for constraint groupings.
`gnomad.utils.constraint.annotate_exploded_vep_for_constraint_groupings`(ht)	Explode a VEP annotation and add constraint grouping fields.
`gnomad.utils.constraint.apply_plateau_models`(...)	Compute the predicted probability observed.
`gnomad.utils.constraint.coverage_correction_expr`(...)	Compute the coverage correction expression.
`gnomad.utils.constraint.apply_models`(...[, ...])	Apply calibration models to compute expected variant counts.
`gnomad.utils.constraint.aggregate_constraint_metrics_expr`(t)	Get an aggregation expression for the sum of specified fields and expressions.
`gnomad.utils.constraint.compute_pli`(ht, ...)	Compute the pLI score using the observed and expected variant counts.
`gnomad.utils.constraint.oe_confidence_interval`(...)	Compute a confidence interval around the observed/expected ratio.
`gnomad.utils.constraint.calculate_raw_z_score`(...)	Compute the signed raw z-score using observed and expected variant counts.
`gnomad.utils.constraint.get_constraint_flags`(...)	Determine the constraint flags that define why constraint will not be calculated.
`gnomad.utils.constraint.calculate_raw_z_score_sd`(...)	Calculate the standard deviation of the raw z-score.
`gnomad.utils.constraint.add_gencode_transcript_annotations`(ht, ...)	Add GENCODE annotations to Table based on transcript id.
`gnomad.utils.constraint.rank_and_assign_bins`(...)	Rank rows by a numeric expression and assign bin labels.
`gnomad.utils.constraint.compute_percentile_thresholds`(ht, ...)	Compute approximate percentile thresholds for a metric expression.
`gnomad.utils.constraint.annotate_bins_by_threshold`(ht, ...)	Annotate rows with bin assignments using pre-computed thresholds.
`gnomad.utils.constraint.rank_array_element_metrics`(ht, ...)	Rank metrics within array elements and annotate rank structs back.

Script containing generic constraint functions that may be used in the constraint pipeline.

gnomad.utils.constraint.COVERAGE_CUTOFF = 30

Minimum median exome coverage differentiating high coverage sites from low coverage sites.

Low coverage sites require an extra calibration when computing the proportion of expected variation.

gnomad.utils.constraint.CLASSIC_LOF_ANNOTATIONS = ('stop_gained', 'splice_donor_variant', 'splice_acceptor_variant'): Classic loss-of-function VEP annotations.

gnomad.utils.constraint.DEFAULT_FIELDS_TO_SUM = ('mu_snp', 'mu', 'observed_variants', 'possible_variants', 'predicted_proportion_observed', 'coverage_correction', 'expected_variants'): Default fields summed when aggregating expected variants by constraint group.

gnomad.utils.constraint.DEFAULT_GENCODE_ANNOTATIONS = ('transcript_id_version', 'gene_id_version', 'level', 'transcript_type', 'start_position', 'end_position'): Default GENCODE annotations added to a Table by transcript id.

gnomad.utils.constraint.get_mu_annotation_expr(ht, mutation_ht, mu_expr='mu_snp')[source]

Get mutation rate annotation expression from the mutation rate Table.

Note

Function expects that ht includes mutation_ht’s key fields. Note that these annotations don’t need to be the keys of ht. The standard keys used are: ‘context’, ‘ref’, ‘alt’, and ‘methylation_level’.

This function converts the mutation rate Table into a dictionary and uses the mutation rate HT key fields in the input Table to get the mutation rate annotation expression. The dictionary is used instead of joining the mutation rate Table to the input Table to avoid unnecessary shuffling given the small size of the mutation rate Table.

Parameters:

ht (Table) – Input Table that will be annotated with the returned mutation rate expression.
mutation_ht (Table) – Mutation rate Table.
mu_expr (Union[str, Float64Expression]) – Mutation rate expression or annotation name in mutation_ht. Default is ‘mu_snp’.

Return type:

Float64Expression

Returns:

Mutation rate annotation expression.

gnomad.utils.constraint.annotate_with_mu(ht, mutation_ht, mu_annotation='mu_snp')[source]

Annotate SNP mutation rate for the input Table.

Note

Function expects that ht includes mutation_ht’s key fields. Note that these annotations don’t need to be the keys of ht. The standard keys used are: ‘context’, ‘ref’, ‘alt’, and ‘methylation_level’.

Parameters:

ht (Table) – Input Table to annotate.
mutation_ht (Table) – Mutation rate Table.
mu_annotation (str) – The name of mutation rate annotation in mutation_ht. Default is ‘mu_snp’.

Return type:

Table

Returns:

Table with mutational rate annotation added.

gnomad.utils.constraint.variant_observed_expr(freq_expr=None, ht=None, singleton=False, max_af=None, count_missing=False)[source]

Return 1 if a variant meets frequency criteria, 0 otherwise.

One of ht or freq_expr must be provided. If freq_expr is not provided, ht.freq is used. When freq_expr is an ArrayExpression, the first element is used.

The returned count is 1 when:

singleton is True and freq_expr.AC == 1.

max_af is not None and freq_expr.AC > 0 and freq_expr.AF <= max_af.

Neither singleton nor max_af is set but freq_expr is available: freq_expr.AC > 0.

Neither singleton nor max_af is set and no freq_expr is available: unconditionally 1.

Parameters:

ht (Optional[Table]) – Input Table. Used to look up freq when freq_expr is None.
freq_expr (Union[StructExpression, ArrayExpression, None]) – StructExpression (or ArrayExpression of structs) with AC and AF fields. When an array, freq_expr[0] is used.
singleton (bool) – Count only singletons (AC == 1). When True, max_af is ignored. Default is False.
max_af (Optional[float]) – Maximum allele frequency threshold. Default is None (no cutoff).
count_missing (bool) – Whether to count a variant whose count expression is missing (e.g. frequency is None). True counts it as 1, False as 0. Default is False.

Return type:

Int32Expression

Returns:

Int32Expression equal to 0 or 1.

gnomad.utils.constraint.variant_observed_and_possible_expr(freq_expr, max_af=None, use_possible_adj=True)[source]

Return per-variant observed and possible count expressions.

For each element of the frequency array, observed_variants is 1 when the variant meets the frequency criteria (AC > 0, and AF <= max_af if set) and 0 otherwise. possible_variants uses the same criteria, except that a missing frequency counts as possible: a variant absent from the dataset could still have been observed. A variant with a defined frequency that fails the criteria (AC == 0, or AF > max_af) is therefore 0 in both, so it is excluded from the denominator as well as the numerator. With max_af set, this removes common variants from the possible counts.

When freq_expr itself is missing (e.g. the variant is not in the dataset), observed_variants is missing and possible_variants is 1, or missing when use_possible_adj is False.

When use_possible_adj is True (default), the possible count is a scalar derived from the first (adj) element of the frequency array. When False, it is an array with one entry per array element, matching the shape of observed_variants.

Parameters:

freq_expr (ArrayExpression) – Array of frequency structs with AC and AF fields.
max_af (Optional[float]) – Maximum allele frequency threshold. Variants above it are excluded from both counts. Default is None (no cutoff).
use_possible_adj (bool) – If True, compute possible count from only the adj (first) frequency element. If False, compute per-downsampling. Default is True.

Return type:

StructExpression

Returns:

Struct with observed_variants (array of int) and possible_variants (int if use_possible_adj, else array of int).

gnomad.utils.constraint.weighted_sum_agg_expr(expr, weight_expr)[source]

Return the weighted aggregate sum of expr weighted by weight_expr.

Both parameters may be scalar or array numeric expressions:

Both arrays: elements are multiplied pairwise and summed per-element across rows (hl.agg.array_sum).

Both scalars: standard hl.agg.sum(expr * weight_expr).

Mixed: the scalar is broadcast across the array elements and summed per-element across rows.

Parameters:

expr (Union[ArrayNumericExpression, NumericExpression]) – Numeric expression (scalar or array) to be weighted.
weight_expr (Union[ArrayNumericExpression, NumericExpression]) – Numeric expression (scalar or array) to weight by.

Return type:

Union[ArrayExpression, NumericExpression]

Returns:

Weighted aggregate sum expression.

gnomad.utils.constraint.counts_agg_expr(freq_expr=None, ht=None, count_singletons=False, max_af=None, count_missing=False)[source]

Return an aggregation expression for variant and singleton counts.

Aggregates per-variant counts (via variant_observed_expr) across rows. Each variant contributes 0 or 1 to the count based on its frequency metadata (AC, AF). The result is a struct with variant_count and, when count_singletons is True, singleton_count.

One of freq_expr or ht must be provided. If freq_expr is not provided, ht.freq is used as the fallback.

The shape of freq_expr controls whether the output counts are scalars or arrays:

StructExpression (single frequency entry with AC / AF): returns scalar variant_count and singleton_count.

ArrayExpression (array of frequency structs, e.g. one per downsampling): returns array-valued counts where each element corresponds to a position in the input array. Internally the array is mapped through variant_observed_expr and summed element-wise with hl.agg.array_sum.

When max_af is set, only variants with AF <= max_af and AC > 0 are counted. Singleton counting (count_singletons=True) counts only variants with AC == 1, independent of the max_af filter.

Parameters:

freq_expr (Union[StructExpression, ArrayExpression, None]) – Frequency expression — an ArrayExpression of structs or a single StructExpression with AC and AF fields. If None, falls back to ht.freq.
ht (Optional[Table]) – Input Table. Used to look up freq when freq_expr is None.
count_singletons (bool) – If True, include a singleton_count field in the result. Default is False.
max_af (Optional[float]) – Maximum allele frequency threshold. Variants with AF > max_af are excluded from variant_count. Does not affect singleton_count. Default is None (no cutoff).
count_missing (bool) – Whether to count a variant whose frequency is missing. True counts it as 1, False as 0. Default is False.

Return type:

StructExpression

Returns:

Aggregation StructExpression with variant_count (and optionally singleton_count). Values are scalars when freq_expr is a StructExpression, or arrays when it is an ArrayExpression.

gnomad.utils.constraint.build_constraint_consequence_groups(csq_expr, lof_modifier_expr, classic_lof_annotations=('stop_gained', 'splice_donor_variant', 'splice_acceptor_variant'), additional_groupings=None, additional_grouping_combinations=None)[source]

Build constraint consequence groups.

Builds constraint groups based on the consequence expression and LoF modifier expression. By default, the following groups are built:

csq_set: synonymous_variant, missense_variant

lof: classic, hc_lc, hc

The resulting meta and corresponding constraint group filters are:

{“csq_set”: “syn”}: synonymous_variant

{“csq_set”: “mis”}: missense_variant

{“lof”: “classic”}: classic LoF annotations (stop_gained, splice_donor_variant, splice_acceptor_variant)

{“lof”: “hc_lc”}: LOFTEE HC or LC modifier

{“lof”: “hc”}: LOFTEE HC modifier only

Additional groupings can be added via additional_groupings, and grouping combinations via additional_grouping_combinations.

Note

A combination of more than one grouping builds the union of the combined groups, not their intersection. For example, with an additional_groupings entry {"am": {"am_per_99": ...}}, the combination ["lof", "am"] produces groups such as {"lof": "hc", "am": "am_per_99"} that contain variants that are HC LoF or in the top 1% of AlphaMissense scores. This differs from frequency metadata (e.g. freq_meta), where multiple keys mean all of them apply.

Parameters:

csq_expr (StringExpression) – VEP most severe consequence expression (e.g., ht.most_severe_consequence).
lof_modifier_expr (StringExpression) – LOFTEE modifier expression (e.g., ht.modifier).
classic_lof_annotations (Tuple[str, ...]) – Classic LoF annotations. Default is ("stop_gained", "splice_donor_variant", "splice_acceptor_variant").
additional_groupings (Optional[Dict[str, Dict[str, BooleanExpression]]]) – Additional groupings to add to the constraint groups. Default is None.
additional_grouping_combinations (Optional[List[List[str]]]) – Additional grouping combinations to add to the constraint groups. Groupings combined in one entry are joined as a union (see the note above). Default is None.

Return type:

Tuple[List[BooleanExpression], List[Dict[str, str]]]

Returns:

Tuple of (constraint group filter expressions, meta dicts).

gnomad.utils.constraint.count_observed_and_possible_by_group(ht, possible_expr, observed_expr, additional_grouping=('methylation_level',), partition_hint=100, weight_exprs=None, additional_agg_sum_exprs=None)[source]

Aggregate observed and possible variant counts by substitution context group.

Groups rows by context, ref, alt, and any fields in additional_grouping, then sums observed and possible counts within each group.

Parameters:

ht (Table) – Input Table with context, ref, alt fields and any fields named in additional_grouping.
possible_expr (Int32Expression) – Per-variant possible count (scalar).
observed_expr (ArrayExpression) – Per-variant observed count (array, one element per array element of the frequency array).
additional_grouping (Union[List[str], Tuple[str, ...]]) – Field names to append to the base (context, ref, alt) grouping. Default is ("methylation_level",).
partition_hint (int) – Target number of partitions for the group_by. Default is 100.
weight_exprs (Union[List[str], Dict[str, Union[ArrayExpression, NumericExpression]], None]) – Weighted sums of possible_expr to include. Pass field names (looked up on ht) or a dict mapping output names to weight expressions. Each produces weighted_sum_agg_expr(possible_expr, weight).
additional_agg_sum_exprs (Union[List[str], Dict[str, Union[ArrayExpression, NumericExpression]], None]) – Extra fields to sum alongside observed/possible. Pass field names (looked up on ht) or a dict mapping output names to expressions. Arrays use hl.agg.array_sum; scalars use hl.agg.sum.

Return type:

Table

Returns:

Grouped Table with observed_variants, possible_variants, and any weighted/additional sum fields.

gnomad.utils.constraint.explode_downsamplings_oe(ht, downsampling_meta, metrics=['syn', 'lof', 'mis'])[source]

Explode observed and expected downsampling counts for each genetic ancestry group and metric.

The input ht must contain struct of downsampling information for genetic ancestry groups under each metric name. For example: ‘lof’: struct {gen_anc_exp: struct {global: array<float64>}.

Parameters:

ht (Table) – Input Table.
metrics (List[str]) – List of metrics to explode. Default is ‘[‘syn’, ‘lof’, ‘mis’]’.
downsampling_meta (Dict[str, List[str]]) – Dictionary containing downsampling metadata. Keys are the genetic ancestry group names and values are the list of downsamplings for that genetic ancestry group. Example: {‘global’: [‘5000’, ‘10000’], ‘afr’: [‘5000’, ‘10000’]}.

Return type:

Table

Returns:

Table with downsampling counts exploded so that observed and expected metric counts for each pair of genetic ancestry groups and downsamplings is a separate row.

gnomad.utils.constraint.annotate_mutation_type(t, context_length=None, num_scan_context_length=100)[source]

Annotate mutation types.

The following annotations are added to the output Table:

cpg
transition
mutation_type - one of “CpG”, “non-CpG transition”, or “transversion”
mutation_type_model

..note:

This function uses the term ‘mutation_type’ because ‘variant_type’ is already used in this repo to indicate a variant’s multiallelic and SNP/indel status.

Parameters:

t (Union[MatrixTable, Table]) – Input Table or MatrixTable.
context_length (Optional[int]) – Length of the ‘context’ annotation in ‘t’. If this is not specified, the value will be determined by examining the first num_scan_context_length values of the ‘context’ annotation. Default is None.
num_scan_context_length (Optional[int]) – Number of values in the ‘context’ annotation to use for determining context_length if it is not specified. If set to None, all values in ‘context’ will be used. Default is 100.

Return type:

Union[MatrixTable, Table]

Returns:

Table with mutation type annotations added.

gnomad.utils.constraint.trimer_from_heptamer(t)[source]

Trim heptamer context to create trimer context.

Parameters:: t (Union[MatrixTable, Table]) – Input MatrixTable or Table with context annotation.
Return type:: Union[MatrixTable, Table]
Returns:: MatrixTable or Table with trimer context annotated.

gnomad.utils.constraint.collapse_strand(t)[source]

Return the deduplicated context by collapsing DNA strands.

Function returns the reverse complement for ‘ref, ‘alt’, and ‘context’ if the reference allele is either ‘G’ or ‘T’.

The following annotations are added to the output Table:

was_flipped - whether the ‘ref, ‘alt’, and ‘context’ were flipped (reverse complement taken)

Parameters:

ht – Input Table.
t (Union[Table, MatrixTable]) –

Return type:

Union[Table, MatrixTable]

Returns:

Table with deduplicated context annotation (ref, alt, context, was_flipped).

gnomad.utils.constraint.transform_methylation_level(methylation_expr='methylation_level', methylation_cutoffs=(0, 5), ht=None)[source]

Transform methylation level into a 0 to 2 scale.

The methylation level is transformed to a 0-2 scale based on the provided cutoffs. The methylation level is assigned a value of 2 if it is greater than the second cutoff, 1 if it is greater than the first cutoff (but less than or equal to the second), and 0 otherwise.

Parameters:

methylation_expr (Union[str, NumericExpression]) – Methylation level expression or annotation name in ht. If methylation_expr is a string, ht must be provided.
methylation_cutoffs (Tuple[Union[int, float], Union[int, float]]) – Tuple of two integers/floats representing the cutoffs for the methylation level transformation. Default is (0, 5).
ht (Optional[Table]) – Input Table. Default is None.

Return type:

Union[Table, NumericExpression]

Returns:

Table with methylation level annotation added or methylation level expression.

gnomad.utils.constraint.assemble_constraint_context_ht(ht, methylation_ht=None, gerp_ht=None, coverage_hts=None, an_hts=None, freq_hts=None, filter_hts=None, transformation_funcs=None)[source]

Assemble context Table with necessary annotations for constraint calculations.

Note

Checks for ‘was_split’ annotation in Table. If not present, splits multiallelic sites.

The following annotations are added to the output Table:

ref: Reference allele.

alt: Alternate allele.

context: Trimer context.

annotations added by annotate_mutation_type(), collapse_strand(), and add_most_severe_csq_to_tc_within_vep_root().

Depending on the annotations provided, the following annotations may also be added:

‘methylation_level’: Methylation level annotation will be added if methylation_ht is provided.

‘gerp’: GERP annotation will be added if gerp_ht is provided.

‘coverage’: Coverage annotations will be added if coverage_hts is provided.

‘AN’: Allele number annotations will be added if an_hts is provided.

‘freq’: Frequency annotations will be added if freq_hts is provided.

‘filters’: Filter annotations will be added if filter_hts is provided.

The transformation_funcs parameter can be used to transform the HTs before adding annotations. The keys should be the annotation names (‘coverage’, ‘gerp’, etc.) and the values should be functions that take the annotation Table and the ht that is being annotated as input and return the transformed and keyed annotation. If not provided, the following default transformations are used:

‘methylation_level’: Uses the ‘methylation_level’ annotation in the Methylation Table after transforming the methylation level to a 0-2 scale using transform_grch37_methylation() or transform_grch38_methylation().

‘gerp’: Uses the ‘S’ annotation in the GERP Table. If ‘S’ is missing, it defaults to 0.

‘coverage’: If necessary, pulls out the first element of coverage statistics (which includes all samples). Relevant to v4, where coverage stats include additional elements to stratify by UK Biobank subset and platforms.

‘AN’: Uses the ‘AN’ annotation in the allele number Table.

‘freq’: Uses the ‘freq’ annotation in the frequency Table.

‘filters’: Uses the ‘filters’ annotation in the filter Table.

The following global annotations may also be added to the output Table:

‘an_globals’: Global allele number annotations ‘strata_sample_count’ and ‘strata_meta’ will be added if an_hts is provided.

‘freq_globals’: Global frequency annotations ‘freq_meta_sample_count’ and ‘freq_meta’ will be added if freq_hts is provided.

Parameters:

ht (Table) – Input context Table with VEP annotation.
methylation_ht (Table) – Optional Table with methylation annotation. Default is None.
gerp_ht (Table) – Optional Table with GERP annotation. Default is None.
coverage_hts (Dict[str, Table]) – An optional Dictionary with key as one of ‘exomes’ or ‘genomes’ and values as corresponding coverage Tables. Default is None.
an_hts (Dict[str, Table]) – A Dictionary with key as one of ‘exomes’ or ‘genomes’ and values as corresponding allele number Tables. Default is None.
transformation_funcs (Optional[Dict[str, Callable]]) – An optional Dictionary to transform the HTs before adding annotations. Default is None, which will use some default transformations as described above.
freq_hts (Dict[str, Table]) –
filter_hts (Dict[str, Table]) –

Return type:

Table

Returns:

Table with sites split and necessary annotations.

gnomad.utils.constraint.calculate_gerp_cutoffs(ht, gerp_expr=None, lower_percentile=0.05, upper_percentile=0.95)[source]

Find GERP cutoffs at the given percentile thresholds.

Note

Uses hl.agg.approx_quantiles, so results are approximate.

Parameters:

ht (Table) – Input Table.
gerp_expr (Optional[Float64Expression]) – GERP score expression. Default is ht.gerp.
lower_percentile (float) – Lower percentile threshold (0-1). Default is 0.05.
upper_percentile (float) – Upper percentile threshold (0-1). Default is 0.95.

Return type:

Tuple[float, float]

Returns:

Tuple of (lower cutoff, upper cutoff) GERP scores.

gnomad.utils.constraint.calibration_model_group_expr(exomes_coverage_expr, cpg_expr, low_cov_cutoff=None, high_cov_cutoff=30, upper_cov_cutoff=None, skip_coverage_model=False, additional_grouping_exprs=None, cpg_in_high_only=False)[source]

Get the calibration model grouping annotation for a variant.

The calibration model expression is a struct with the following fields:

high_or_low_coverage: Whether the variant belongs to the high or low coverage calibration model. The variant is assigned to the high coverage model if the exome coverage is greater than or equal to ‘high_cov_cutoff’ and less than or equal to ‘upper_cov_cutoff’ (if provided). The variant is assigned to the low coverage model if skip_coverage_model is False and the exome coverage is greater than ‘low_cov_cutoff’ (if provided) and less than ‘high_cov_cutoff’.

model_group: Struct of the fields the models are grouped by:

cpg: Whether the variant is a CpG (cpg_expr).

Any fields in additional_grouping_exprs (e.g. genomic_region).

The expression is missing for variants that belong to neither model.

Parameters:

exomes_coverage_expr (Int32Expression) – Exome coverage expression.
cpg_expr (BooleanExpression) – CpG expression.
low_cov_cutoff (Optional[int]) – Low coverage cutoff. Default is None.
high_cov_cutoff (int) – High coverage cutoff. Default is COVERAGE_CUTOFF.
upper_cov_cutoff (Optional[int]) – Upper coverage cutoff. Default is None.
skip_coverage_model (bool) – Whether to skip the coverage model. Default is False.
additional_grouping_exprs (Optional[Dict[str, StringExpression]]) – Optional Dictionary of additional expressions to group by. Default is None.
cpg_in_high_only (bool) – Whether to only include cpg_expr in the model grouping for the high coverage model. Default is False.

Return type:

StructExpression

Returns:

StructExpression with high_or_low_coverage and model_group fields.

gnomad.utils.constraint.build_models(ht, coverage_expr, weighted=False, keys=('context', 'ref', 'alt', 'methylation_level'), model_group_expr=None, high_cov_definition=30, upper_cov_cutoff=None, skip_coverage_model=False, log10_coverage=True, additional_grouping=())[source]

Build coverage and plateau models.

This function builds models (plateau_models) using linear regression to calibrate mutation rate estimates against the proportion observed of each substitution, context, and methylation level in ht.

Two plateau models are fit, one for CpG transitions, and one for the remainder of sites (transversions and non CpG transitions).

The plateau models only consider high coverage sites, or sites above a median coverage of high_cov_definition and median coverage below upper_cov_cutoff.

Plateau model: adjusts proportion of expected variation based on location in the genome and CpG status.

The x and y of the plateau models:

x: mu_snp - mutation rate
y: proportion observed (‘observed_variants’ / ‘possible_variants’)

This function also builds models (coverage models) to calibrate the proportion of expected variation at low coverage sites (sites below high_cov_definition).

The coverage models are built by creating a scaling factor across all high coverage sites, applying this ratio to the low coverage sites, and running a linear regression.

Coverage model: corrects proportion of expected variation at low coverage sites. Low coverage sites are defined as sites with median coverage < high_cov_definition.

The x and y of the coverage model:

x: groupings of exome coverage at low coverage sites (log10 transformed if requested)

y: sum(‘observed_variants’) / (high_coverage_scale_factor * sum(‘possible_variants’ * ‘mu_snp’)) at low coverage sites

high_coverage_scale_factor = sum(‘observed_variants’) / sum(‘possible_variants’ * ‘mu_snp’) at high coverage sites

Note

This function expects that the input Table (ht) contains only high quality synonymous variants below 0.1% frequency.

The following fields are expected in ht:

context: trinucleotide genomic context.

ref: the reference allele.

alt: the alternate allele.

methylation_level: methylation level.

mu_snp: mutation rate.

cpg: whether the variant is a CpG.

observed_variants: the number of observed variants in the dataset for each variant. Note that the term “variant” here refers to a specific substitution, context, methylation level, and coverage combination.

possible_variants: the number of possible variants in the dataset for each variant.

Parameters:

ht (Table) – Input Table.
coverage_expr (Int32Expression) – Expression that defines the coverage metric.
weighted (bool) – Whether to weight the plateau models (a linear regression model) by ‘possible_variants’. Default is False.
keys (Tuple[str, ...]) – Annotations used to group observed and possible variant counts. Default is (“context”, “ref”, “alt”, “methylation_level”).
model_group_expr (Optional[StructExpression]) – Expression with high_or_low_coverage and model_group annotations to group variants into high or low coverage models. If model_group includes a genomic_region field, the coverage model is built only on sites where it is “autosome_or_par”; otherwise all sites are used. If not provided, the calibration_model_group_expr function is used to define the grouping, which does not include genomic_region.
high_cov_definition (int) – Coverage threshold for high/low classification. Sites with coverage above this cutoff are considered well covered. Default is COVERAGE_CUTOFF.
upper_cov_cutoff (Optional[int]) – Upper coverage cutoff. Sites with coverage above this cutoff are excluded from the high coverage Table. Default is None.
skip_coverage_model (bool) – Whether to skip generating the coverage model. If set to True, None is returned instead of the coverage model. Default is False.
log10_coverage (bool) – Whether to convert coverage sites with log10 when building the coverage model. Default is True.
additional_grouping (Tuple[str, ...]) – Additional annotations to group by before counting the observed and possible variants. Default is ().

Return type:

Tuple[Optional[ArrayNumericExpression], DictExpression]

Returns:

Tuple of (coverage model, plateau models). The coverage model is an [intercept, slope] ArrayNumericExpression, or None when skip_coverage_model is True.

gnomad.utils.constraint.build_plateau_models(mu_snp_expr, observed_variants_expr, possible_variants_expr, weighted=False, cpg_expr=None, model_group_expr=None)[source]

Build plateau models to calibrate mutation rate against proportion observed.

Fits a linear regression of observed_variants_expr / possible_variants_expr on mu_snp_expr. When either observed or possible expressions are arrays (e.g., one model per downsampling), the regression is applied element-wise via hl.agg.array_agg.

When model_group_expr or cpg_expr is provided, the result is a DictExpression keyed by the grouping struct. Otherwise the result is the regression beta directly.

Parameters:

mu_snp_expr (Float64Expression) – Mutation rate expression.
observed_variants_expr (Union[ArrayExpression, Int64Expression]) – Observed variant counts (scalar or array).
possible_variants_expr (Union[ArrayExpression, Int64Expression]) – Possible variant counts (scalar or array).
weighted (bool) – If True, use weighted least squares with possible_variants_expr as weights. Default is False.
cpg_expr (Optional[BooleanExpression]) – Boolean expression indicating CpG sites. When provided, adds a cpg field to the grouping struct.
model_group_expr (Optional[StructExpression]) – Struct expression to group by in the aggregation.

Return type:

Union[DictExpression, ArrayExpression, StructExpression]

Returns:

Regression betas, optionally grouped by model_group_expr (and/or cpg_expr).

gnomad.utils.constraint.build_coverage_model(low_coverage_oe_expr, coverage_expr, log10_coverage=False)[source]

Build coverage model.

This function uses linear regression to build a model of coverage to correct proportion of expected variation at low coverage sites.

The x and y of the coverage model:

x: coverage_expr

y: low_coverage_oe_expr

Parameters:

low_coverage_oe_expr (Float64Expression) – The Float64Expression of observed:expected ratio for a given coverage level.
coverage_expr (Float64Expression) – The Float64Expression of the coverage expression.
log10_coverage (bool) – Whether to convert coverage sites by log10 when building the coverage model. Default is False.

Return type:

StructExpression

Returns:

StructExpression with intercept and slope of the model.

gnomad.utils.constraint.get_all_gen_anc_lengths(ht, gen_ancs, obs_expr)[source]

Get the minimum length of observed variant counts array for each genetic ancestry group downsampling.

The observed variant counts for each genetic ancestry group in gen_ancs are specified by annotations on the obs_expr expression.

The function also performs a check that arrays of variant counts within genetic ancestry group downsamplings all have the same lengths.

Parameters:

ht (Table) – Input Table containing obs_expr.
gen_ancs (Tuple[str, ...]) – Genetic ancestry groups used to categorize observed variant counts in downsamplings.
obs_expr (StructExpression) – Expression for the genetic ancestry group observed variant counts. Should be a struct containing an array for each gen_anc in gen_ancs.

Return type:

List[Tuple[int, str]]

Returns:

List of (length, gen_anc) tuples, one per group in gen_ancs, where length is the minimum array length for that group.

gnomad.utils.constraint.get_constraint_grouping_expr(vep_annotation_expr, coverage_expr=None, include_transcript_group=True, include_canonical_group=True, include_mane_select_group=False)[source]

Collect annotations used for constraint groupings.

Function collects the following annotations:

annotation - most_severe_consequence from vep_annotation_expr

modifier - first non-missing of lof or polyphen_prediction from vep_annotation_expr, or the literal “None”

gene - gene_symbol from vep_annotation_expr

gene_id - gene_id from vep_annotation_expr

coverage - exome coverage if coverage_expr is specified

transcript - transcript_id from vep_annotation_expr (added when include_transcript_group is True)

canonical - from vep_annotation_expr (added when include_canonical_group is True)

mane_select - from vep_annotation_expr (added when include_mane_select_group is True)

Note

This function expects that the following fields are present in vep_annotation_expr:

lof

most_severe_consequence

gene_symbol

gene_id

polyphen_prediction (optional; missing used if absent)

transcript_id (if include_transcript_group is True)

canonical (if include_canonical_group is True)

mane_select (if include_mane_select_group is True)

Parameters:

vep_annotation_expr (StructExpression) – StructExpression of VEP annotation.
coverage_expr (Optional[Int32Expression]) – Int32Expression of exome coverage. Default is None.
include_transcript_group (bool) – Whether to include the transcript annotation in the groupings. Default is True.
include_canonical_group (bool) – Whether to include canonical annotation in the groupings. Default is True.
include_mane_select_group (bool) – Whether to include mane_select annotation in the groupings. Default is False.

Return type:

Dict[str, Union[StringExpression, Int32Expression, BooleanExpression]]

Returns:

Dict mapping annotation names to Hail expressions.

gnomad.utils.constraint.annotate_exploded_vep_for_constraint_groupings(ht, coverage_expr=None, vep_annotation='transcript_consequences', include_canonical_group=True, include_mane_select_group=False)[source]

Explode a VEP annotation and add constraint grouping fields.

Explodes the specified VEP annotation (vep_annotation) and adds the following annotations via get_constraint_grouping_expr:

annotation - most_severe_consequence from vep_annotation

modifier - first non-missing of lof or polyphen_prediction from vep_annotation, or the literal “None”

gene - gene_symbol from vep_annotation

gene_id - gene_id from vep_annotation

coverage - exome coverage if coverage_expr is specified

transcript - transcript_id from vep_annotation (added when vep_annotation is “transcript_consequences”)

canonical - from vep_annotation (added when include_canonical_group is True)

mane_select - from vep_annotation (added when include_mane_select_group is True)

Note

This function expects that a vep annotation is present in ht.

Note

Only the grouping fields are kept from the VEP annotation. The returned Table no longer has the vep annotation or the exploded vep_annotation struct, which keeps per-variant tables built from it small. All other fields of ht are kept.

Parameters:

ht (Table) – Input Table.
coverage_expr (Optional[Int32Expression]) – Expression that defines the coverage metric. Default is None.
vep_annotation (str) – Name of annotation in vep (one of "transcript_consequences" and "worst_csq_by_gene") that will be used for obtaining constraint annotations. Default is "transcript_consequences".
include_canonical_group (bool) – Whether to include canonical annotation in the groupings. Default is True. Ignored unless vep_annotation is "transcript_consequences".
include_mane_select_group (bool) – Whether to include mane_select annotation in the groupings. Default is False. Ignored unless vep_annotation is "transcript_consequences".

Return type:

Tuple[Table, Tuple[str, ...]]

Returns:

Tuple of (Table with one row per element of vep_annotation and the grouping fields added in place of the VEP annotation, names of the added grouping fields).

gnomad.utils.constraint.apply_plateau_models(mu_expr, plateau_models_expr)[source]

Compute the predicted probability observed.

The predicted probability observed is computed as the mutation rate adjusted by the plateau model.

Parameters:

mu_expr (Float64Expression) – Mutation rate expression.
plateau_models_expr (ArrayExpression) – This can be either a single plateau model, where the first element is the intercept and the second element is the slope, or an array of plateau models.

Return type:

Union[ArrayExpression, Float64Expression]

Returns:

Float64Expression of the predicted probability observed, or an ArrayExpression of them when plateau_models_expr is an array of models.

gnomad.utils.constraint.coverage_correction_expr(coverage_expr, coverage_model, low_coverage_expr=None, coverage_cutoff=None, log10_coverage=False)[source]

Compute the coverage correction expression.

Note

One and only one of low_coverage_expr or coverage_cutoff must be specified.

The coverage correction expression is computed as follows:

If the coverage is 0, the coverage correction is 0.

If the low coverage expression (low_coverage_expr) is True, or the coverage (coverage_expr) is below the coverage cutoff (coverage_cutoff), the coverage model is applied to the coverage.

Otherwise, the coverage correction is 1.

Parameters:

coverage_expr (Float64Expression) – Float64Expression of the coverage.
coverage_model (Union[Tuple[float, float], ArrayNumericExpression]) – Coverage model as [intercept, slope], either a tuple of floats or an ArrayNumericExpression (as returned by build_models).
low_coverage_expr (Optional[BooleanExpression]) – Optional BooleanExpression indicating whether the site is a low coverage site, and the coverage model should be applied. Default is None.
coverage_cutoff (Optional[int]) – Optional coverage cutoff. If specified, the coverage model is applied to sites with coverage below this cutoff. Default is None.
log10_coverage (bool) – Whether to convert coverage sites by log10 when applying the coverage model. Default is False.

Return type:

Float64Expression

Returns:

Float64Expression of the coverage correction.

gnomad.utils.constraint.apply_models(mu_expr, plateau_models_expr, possible_variants_expr, coverage_model=None, coverage_expr=None, cpg_expr=None, model_group_expr=None, high_cov_definition=30, log10_coverage=True)[source]

Apply calibration models to compute expected variant counts.

Applies plateau and (optionally) coverage models to produce a struct with mu, predicted_proportion_observed, expected_variants, and (when a coverage model is provided) coverage_correction.

Note

When a coverage model is applied, a row whose model_group_expr is missing gets missing mu, expected_variants, and coverage_correction (unless its coverage is 0, which always gives a correction of 0). calibration_model_group_expr is missing for sites in neither model, such as sites above upper_cov_cutoff or at or below low_cov_cutoff. hl.agg.sum skips missing values, so filter these rows out before applying the models if they should not silently drop out of aggregated expected counts.

Parameters:

mu_expr (Float64Expression) – Mutation rate expression.
plateau_models_expr (ArrayExpression) – Single plateau model (array of [intercept, slope]) or an array of plateau models.
possible_variants_expr (Int64Expression) – Possible variant counts to multiply the predicted proportion observed by.
coverage_model (Union[Tuple[float, float], ArrayNumericExpression, None]) – Coverage model as [intercept, slope], either a tuple of floats or an ArrayNumericExpression (as returned by build_models). Default is None.
coverage_expr (Optional[Int32Expression]) – Int32Expression of the coverage. Required when coverage_model is provided or model_group_expr is None.
cpg_expr (Optional[BooleanExpression]) – BooleanExpression indicating CpG sites. Required when model_group_expr is None.
model_group_expr (Optional[StructExpression]) – Expression with high_or_low_coverage annotation to group variants into high or low coverage models. If not provided, calibration_model_group_expr is used. Default is None.
high_cov_definition (int) – Coverage threshold for high/low classification. Default is COVERAGE_CUTOFF.
log10_coverage (bool) – Whether to log10-transform coverage when applying the coverage model. Default is True.

Return type:

StructExpression

Returns:

StructExpression with “mu”, “predicted_proportion_observed”, “expected_variants”, and optionally “coverage_correction”.

gnomad.utils.constraint.aggregate_constraint_metrics_expr(t, fields_to_sum=('mu_snp', 'mu', 'observed_variants', 'possible_variants', 'predicted_proportion_observed', 'coverage_correction', 'expected_variants'), additional_exprs_to_sum=None)[source]

Get an aggregation expression for the sum of specified fields and expressions.

An aggregate sum or array sum is created for each field in fields_to_sum.

Parameters:

t (Union[Table, StructExpression]) – Input Table or StructExpression.
fields_to_sum (Union[List[str], Tuple[str, ...]]) – Fields in t to get an aggregate sum expression for. Default is DEFAULT_FIELDS_TO_SUM.
additional_exprs_to_sum (Optional[Dict[str, Expression]]) – Dictionary of additional expressions to get an aggregate sum expression for. Field names are the keys and expressions are the values. Default is None.

Return type:

StructExpression

Returns:

StructExpression with the sum of each field and expression.

gnomad.utils.constraint.compute_pli(ht, obs_expr, exp_expr, expected_values=None, min_diff_convergence=0.001)[source]

Compute the pLI score using the observed and expected variant counts.

Full details on pLI can be found in the ExAC paper: Lek, M., Karczewski, K., Minikel, E. et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature 536, 285–291 (2016).

pLI is the probability of being loss-of-function intolerant, and this function computes that probability using the expectation-maximization (EM) algorithm.

We assume a 3 state model, where each gene fits into one of three categories with respect loss-of-function variation sensitivity:

Null: where protein truncating variation is completely tolerated by natural selection.

Recessive (Rec): where heterozygous pLoFs are tolerated but homozygous pLoFs are not.

Haploinsufficient (LI): where heterozygous pLoFs are not tolerated.

The function requires the expected amount of loss-of-function depletion for each of these states. The default provided is based on the observed depletion of protein-truncating variation in the Blekhman autosomal recessive and ClinGen dosage sensitivity gene sets (Supplementary Information Table 12 of the above reference):

Null: 1.0, assume tolerant genes have the expected amount of truncating variation.

Rec: 0.463, derived from the empirical mean observed/expected rate of truncating variation for recessive disease genes (0.463).

LI: 0.089, derived from the empirical mean observed/expected rate of truncating variation for severe haploinsufficient genes.

The output StructExpression will include the following annotations:

pLI: Probability of loss-of-function intolerance; probability that transcript falls into distribution of haploinsufficient genes.

pNull: Probability that transcript falls into distribution of unconstrained genes.

pRec: Probability that transcript falls into distribution of recessive genes.

Parameters:

ht (Table) – Input Table containing obs_expr and exp_expr.
obs_expr (Int64Expression) – Expression for the number of observed variants on each gene or transcript in ht.
exp_expr (Float64Expression) – Expression for the number of expected variants on each gene or transcript in ht.
expected_values (Optional[Dict[str, float]]) – Dictionary containing the expected values for ‘Null’, ‘Rec’, and ‘LI’ to use as starting values.
min_diff_convergence (float) – Minimum iteration change in LI to consider the EM model convergence criteria as met. Default is 0.001.

Return type:

StructExpression

Returns:

StructExpression for pLI scores.

gnomad.utils.constraint.oe_confidence_interval(obs_expr, exp_expr, alpha=0.05, method='gamma')[source]

Compute a confidence interval around the observed/expected ratio.

Two methods are available:

"gamma" (default): uses hl.qgamma to compute exact quantiles of the Gamma posterior. Fast and precise.
"poisson": sweeps a discretized Poisson likelihood over the OE parameter space [0, 2). Retained for backwards compatibility.

Note

The default changed from the discretized Poisson computation to "gamma", so calls that do not pass method return different bounds than before. "gamma" requires Hail >= 0.2.137, which added hl.qgamma; on older Hail versions pass method="poisson".

Parameters:

obs_expr (Int64Expression) – Observed variant count expression.
exp_expr (Float64Expression) – Expected variant count expression.
alpha (float) – Significance level for the confidence interval. Default is 0.05 (90% CI).
method (str) – CI method — "gamma" or "poisson". Default is "gamma".

Return type:

StructExpression

Returns:

Struct with lower and upper bounds.

gnomad.utils.constraint.calculate_raw_z_score(obs_expr, exp_expr)[source]

Compute the signed raw z-score using observed and expected variant counts.

The raw z-scores are positive when the transcript had fewer variants than expected, and are negative when transcripts had more variants than expected.

Parameters:

obs_expr (Int64Expression) – Observed variant count expression.
exp_expr (Float64Expression) – Expected variant count expression.

Return type:

Float64Expression

Returns:

Raw z-score expression.

gnomad.utils.constraint.get_constraint_flags(exp_expr, raw_z_expr, raw_z_lower_threshold=-5.0, raw_z_upper_threshold=5.0, flag_postfix='')[source]

Determine the constraint flags that define why constraint will not be calculated.

Flags which are added:

“no_exp_{flag_postfix}” - for genes that have missing or zero expected variants.
“outlier_{flag_postfix}” - for genes that are raw z-score outliers: (raw_z_expr < raw_z_lower_threshold) or (raw_z_expr > raw_z_upper_threshold).

Parameters:

exp_expr (Float64Expression) – Expression for the expected variant counts of pLoF, missense, or synonymous variants.
raw_z_expr (Float64Expression) – Expression for the signed raw z-score of pLoF, missense, or synonymous variants.
raw_z_lower_threshold (Union[float, Float64Expression, None]) – Lower threshold for the raw z-score. When raw_z_expr is less than this threshold it is considered an ‘outlier’. Default is -5.0.
raw_z_upper_threshold (Union[float, Float64Expression, None]) – Upper threshold for the raw z-score. When raw_z_expr is greater than this threshold it is considered an ‘outlier’. Default is 5.0.
flag_postfix (str) – Postfix to add to the end of the constraint flag names.

Return type:

Dict[str, Expression]

Returns:

Dictionary containing expressions for constraint flags.

gnomad.utils.constraint.calculate_raw_z_score_sd(raw_z_expr, flag_expr, mirror_neg_raw_z=True)[source]

Calculate the standard deviation of the raw z-score.

When using mirror_neg_raw_z is True, all the negative raw z-scores (defined by raw_z_expr) are combined with those same z-scores multiplied by -1 (to create a mirrored distribution).

Parameters:

raw_z_expr (Float64Expression) – Expression for the raw z-score.
flag_expr (StringExpression) – Expression for the constraint flags. z-score will not be calculated if flags are present.
mirror_neg_raw_z (bool) – Whether the standard deviation should be computed using a mirrored distribution of negative raw_z_expr.

Return type:

Expression

Returns:

StructExpression containing standard deviation of the raw z-score and the z-score.

gnomad.utils.constraint.add_gencode_transcript_annotations(ht, gencode_ht, annotations=('transcript_id_version', 'gene_id_version', 'level', 'transcript_type', 'start_position', 'end_position'), remove_y_par=True)[source]

Add GENCODE annotations to Table based on transcript id.

In addition to the annotations specified by annotations, the following computed annotations are always added:

chromosome

cds_length

num_coding_exons

Parameters:

ht (Table) – Input Table.
gencode_ht (Table) – Table with GENCODE annotations.
annotations (Union[List[str], Tuple[str, ...]]) – GENCODE annotations to add. Default is DEFAULT_GENCODE_ANNOTATIONS.
remove_y_par (bool) – Whether to remove features for the Y chromosome PAR regions. Default is True because the Y chromosome PAR regions are typically not included in the constraint calculations and both chrX and chrY will have the same ‘transcript_id’ field for these regions. This parameter can only be True if gencode_ht includes a ‘transcript_id_version’ field because Y_PAR (or PAR_Y, depending on the GENCODE release) is included in the version of the transcript, which has been stripped from the ‘transcript_id’ field.

Return type:

Table

Returns:

Table with transcript annotations from GENCODE added.

gnomad.utils.constraint.rank_and_assign_bins(value_expr, bin_granularities=None, prefix='')[source]

Rank rows by a numeric expression and assign bin labels.

Rank-based binning: every row receives a unique position in the sorted order, and bins are derived from that position. This differs from threshold-based binning (see annotate_bins_by_threshold()), where bins are assigned by comparing values against pre-computed boundary values.

Rows are ordered ascending by value_expr. Each row is assigned a 0-based {prefix}rank and a {prefix}bin_{name} field for every entry in bin_granularities, computed as hl.int(rank * multiplier / n_rows).

Rows where value_expr is missing are excluded from the ranking and receive missing rank and bin annotations. n_rows is therefore the number of rows with a defined value_expr, so bins are evenly sized and the full range of bins is populated.

Used by rank_array_element_metrics() to rank metrics within array elements.

Note

This function evaluates its input when it is called. It needs the number of rows with a defined value_expr for the bin denominator, so it counts them, and it caches the filtered and the ranked rows so the filter and the sort are not recomputed. In Hail, Table.cache() ignores the storage level and writes a checkpoint to a persist_Table* path under the temporary directory set by hl.init(tmp_dir=...), so with a GCS tmp_dir these Tables are written to GCS. They hold only value_expr and the key, and are not removed when the function returns.

Parameters:

value_expr (Float64Expression) – Numeric expression to rank by (ascending).
bin_granularities (Optional[Dict[str, int]]) – Mapping of bin name to multiplier. Each entry produces a {prefix}bin_{name} field. Default is {"percentile": 100, "decile": 10, "sextile": 6}.
prefix (str) – String prepended to rank and bin_{name} field names. Default is "" (no prefix).

Return type:

StructExpression

Returns:

Struct with {prefix}rank and {prefix}bin_{name} fields for each entry in bin_granularities.

gnomad.utils.constraint.compute_percentile_thresholds(ht, metric_expr, outlier_expr=None, transcript_filter_expr=None, percentiles=(1, 5, 10, 15, 25, 50, 75), quantile_k=1000)[source]

Compute approximate percentile thresholds for a metric expression.

Threshold-based binning, step 1: computes the boundary values that define bin edges. Its values, keyed by (granularity, metric) instead of percentile, form the thresholds argument of annotate_bins_by_threshold() (step 2), which assigns each row to a bin; see that function for an example.

This two-step approach differs from rank-based binning (see rank_and_assign_bins()), where every row receives a unique position in the sorted order. Threshold-based binning allows thresholds to be computed on a filtered subset (e.g., representative transcripts) and then applied to all rows, so multiple rows can share the same bin.

Optionally filters to a subset of rows and excludes outliers, then computes approximate quantile thresholds at the requested percentiles in a single aggregation pass.

Note

Uses hl.agg.approx_quantiles, so results are approximate. Increase quantile_k for higher accuracy.

Parameters:

ht (Table) – Input Table.
metric_expr (Float64Expression) – Float expression to compute thresholds for. Must be defined on ht.
outlier_expr (Optional[BooleanExpression]) – Optional boolean expression that is True for rows to exclude. When None (default), no outlier filtering is applied.
transcript_filter_expr (Optional[BooleanExpression]) – Optional boolean expression that is True for rows to include. When None (default), all rows are included.
percentiles (Tuple[float, ...]) – Percentile values (0-100) at which to compute thresholds. Default is (1, 5, 10, 15, 25, 50, 75).
quantile_k (int) – Accuracy parameter for hail.expr.aggregators.approx_quantiles(). Default is 1000.

Return type:

Dict[float, float]

Returns:

Dict mapping each percentile to its threshold value.

gnomad.utils.constraint.annotate_bins_by_threshold(ht, metric_exprs, thresholds, granularities, field_name='constraint_bins')[source]

Annotate rows with bin assignments using pre-computed thresholds.

Threshold-based binning, step 2: assigns each row to a bin by comparing its metric value against boundary values produced by compute_percentile_thresholds() (step 1). For each (granularity, metric) pair, the bin equals the number of threshold boundaries the value exceeds. Bin 0 is the most constrained (below all thresholds); bin N means the value exceeds all N boundaries.

This differs from rank-based binning (see rank_and_assign_bins()), where every row gets a unique position. Here, multiple rows can share a bin, and the thresholds may have been derived from a different subset of rows than those being annotated.

compute_percentile_thresholds() returns a dict of percentile to value, so thresholds is built from it by selecting the percentiles for each granularity. Example use, with decile bins for two metrics:

deciles = [10, 20, 30, 40, 50, 60, 70, 80, 90]
metric_exprs = {"lof": ht.lof_oe_upper, "mis": ht.mis_oe_upper}
thresholds = {}
for metric, expr in metric_exprs.items():
    pct_to_value = compute_percentile_thresholds(ht, expr, percentiles=deciles)
    thresholds[("decile", metric)] = [pct_to_value[p] for p in deciles]

ht = annotate_bins_by_threshold(ht, metric_exprs, thresholds, ["decile"])

Parameters:

ht (Table) – Input Table.
metric_exprs (Dict[str, Float64Expression]) – Mapping of metric name to the Float64Expression to bin (e.g. {"lof": ht.lof_oe_upper, "mis": ht.mis_oe_upper}).
thresholds (Dict[Tuple[str, str], List[float]]) – Mapping of (granularity, metric) to a list of threshold values in ascending order, built from compute_percentile_thresholds() output (see the example above).
granularities (Union[List[str], Tuple[str, ...]]) – Granularity names to iterate over (e.g. ["decile", "ventile"]). Each must appear as the first element of at least one key in thresholds.
field_name (str) – Name of the struct field to annotate on ht. Default is "constraint_bins".

Return type:

Table

Returns:

Table with an added struct field containing per-granularity, per-metric bin assignments.

gnomad.utils.constraint.rank_array_element_metrics(ht, array_field, element_value_fn, filter_fn=None, bin_granularities=None, rank_field_prefix='')[source]

Rank metrics within array elements and annotate rank structs back.

Rank-based binning for array fields: applies rank_and_assign_bins() independently to each element of an array field. For each element, element_value_fn extracts named metric values, which are ranked across rows (optionally on a filtered subset via filter_fn). Each array element is then annotated with {metric_name}_rank structs containing rank and bin fields.

Rows not matching filter_fn get missing rank annotations.

This differs from threshold-based binning (see compute_percentile_thresholds() and annotate_bins_by_threshold()), where bins are assigned by comparing values against pre-computed boundaries rather than sorted position.

Note

This function evaluates its input when it is called. It caches the full input Table (with an added row index) and the extracted metric values, reads one row to get the array length, and each rank_and_assign_bins() call it makes caches its own ranked rows. Caching avoids recomputing ht and filter_fn for every metric, but in Hail Table.cache() writes a checkpoint under the temporary directory (see rank_and_assign_bins()), so this writes a copy of ht there. The copies are not removed when the function returns.

Parameters:

ht (Table) – Input Table.
array_field (str) – Name of the array field on ht.
element_value_fn (Callable[[StructExpression], Dict[str, Float64Expression]]) – Function that takes an array element (StructExpression) and returns a dict mapping metric names to Float64Expressions to rank. Applied identically to every element.
filter_fn (Optional[Callable[[Table], BooleanExpression]]) – Optional function that takes a Table and returns a BooleanExpression to filter rows before ranking. When None, all rows are ranked.
bin_granularities (Optional[Dict[str, int]]) – Bin granularities passed to rank_and_assign_bins().
rank_field_prefix (str) – Prefix for the rank and bin_{name} sub-fields within each rank struct. Passed through to rank_and_assign_bins(). Default is "" (no prefix).

Return type:

Table

Returns:

Table with {metric_name}_rank structs added to each array element. The table is returned with its original key restored.