gnomad_qc.v5.variant_qc.isolation_forest
Script to run the GATK isolation forest variant QC model on an AS-annotated sites VCF.
SNPs and INDELs are modeled separately (GATK --mode): for each variant type the
workflow runs ExtractVariantAnnotations -> TrainVariantAnnotationsModel (PYTHON_IFOREST)
-> ScoreVariantAnnotations (scattered over intervals). The two per-mode scored VCFs are
merged into a single result HT in the --load-iforest step.
usage: gnomad_qc.v5.variant_qc.isolation_forest.py [-h] [-o]
[--environment {rwb,batch}]
--model-id MODEL_ID
[--test]
[--test-chrom TEST_CHROM [TEST_CHROM ...]]
[--test-on-v4]
[--export-v4-test-vcf]
[--export-only]
--batch-billing-project
BATCH_BILLING_PROJECT
--gcp-billing-project
GCP_BILLING_PROJECT
[--batch-suffix BATCH_SUFFIX]
[--app-name APP_NAME]
[--driver-cores DRIVER_CORES]
[--driver-memory DRIVER_MEMORY]
[--worker-cores WORKER_CORES]
[--worker-memory WORKER_MEMORY]
[--gatk-image GATK_IMAGE]
[--scatter-count SCATTER_COUNT]
--out-vcf-name OUT_VCF_NAME
[--hyperparameters-json HYPERPARAMETERS_JSON]
[--transmitted-singletons]
[--sibling-singletons]
[--adj] [--load-iforest]
[--load-only]
[--n-partitions N_PARTITIONS]
[--header-path HEADER_PATH]
[--array-elements-required]
Named Arguments
- -o, --overwrite
Whether to overwrite existing outputs (with –load-only, only overwrite the result HT).
Default: False
- --environment
Possible choices: rwb, batch
Compute environment.
Default: “batch”
- --model-id
Model ID for the isolation forest run. Must start with ‘if’ (e.g. if_test1).
testing options
- --test
Filter to –test-chrom and use a small scatter count for testing.
Default: False
- --test-chrom
Contig(s) to score under –test/–test-on-v4. Must match the contigs in the test info VCF.
Default: [‘chr21’]
- --test-on-v4
Use gnomAD v4 sites/true-positive VCFs as inputs (v5 data unavailable).
Default: False
- --export-v4-test-vcf
Export the v4 split info HT filtered to –test-chrom as the –test-on-v4 sites VCF. Requires –test-on-v4.
Default: False
- --export-only
Stop after –export-v4-test-vcf instead of continuing to the Batch.
Default: False
batch/GCP configuration
- --batch-billing-project
Hail Batch billing project.
- --gcp-billing-project
GCP billing project for requester-pays buckets.
- --batch-suffix
String to append to the Batch name.
Default: “”
- --app-name
Job name for the batch/QoB backend.
- --driver-cores
Number of driver cores. The load step has been run with 4-8.
- --driver-memory
Driver memory. Must be ‘highmem’ for –load-iforest/–load-only.
- --worker-cores
Number of worker cores.
- --worker-memory
Worker memory.
GATK scoring
- --gatk-image
GATK docker image.
Default: “us.gcr.io/broad-gatk/gatk:4.6.1.0”
- --scatter-count
Number of intervals to scatter scoring across.
Default: 100
- --out-vcf-name
Base name for scored VCF shards.
- --hyperparameters-json
Optional GATK model hyperparameters JSON path.
training/calibration set
- --transmitted-singletons
Include transmitted singletons as a training/calibration set.
Default: False
- --sibling-singletons
Include sibling singletons as a training/calibration set.
Default: False
- --adj
Use adj genotypes for the true-positive singletons VCF.
Default: False
load results HT
- --load-iforest
Merge the scored VCFs into a variant QC result HT.
Default: False
- --load-only
Skip the GATK Batch and only run the load step.
Default: False
- --n-partitions
Number of partitions for the result HT.
Default: 5000
- --header-path
Optional header file to use when importing the scored VCFs.
- --array-elements-required
Set array_elements_required=True when importing the scored VCFs.
Default: False
Module Functions
GATK image. |
|
|
Run GATK ExtractVariantAnnotations for one variant type mode. |
|
Run GATK TrainVariantAnnotationsModel (PYTHON_IFOREST) for one variant type mode. |
|
Split the calling contigs into |
|
Run GATK ScoreVariantAnnotations for one variant type mode over one interval. |
|
Guard reuse of outputs under |
|
Add the per-mode extract/train/score jobs (SNP and INDEL) to the Batch. |
|
Export telomere/centromere intervals to a GATK-readable |
|
Export the v4 split info HT, filtered to |
Run a Batch to completion and raise if it did not succeed. |
|
|
Merge the per-mode scored VCFs into a single variant QC result HT. |
|
Raise if any input variant in the scored region is missing from the output or unscored. |
Run the isolation forest variant QC workflow. |
|
|
Get script argument parser. |
Script to run the GATK isolation forest variant QC model on an AS-annotated sites VCF.
SNPs and INDELs are modeled separately (GATK --mode): for each variant type the
workflow runs ExtractVariantAnnotations -> TrainVariantAnnotationsModel (PYTHON_IFOREST)
-> ScoreVariantAnnotations (scattered over intervals). The two per-mode scored VCFs are
merged into a single result HT in the --load-iforest step.
- gnomad_qc.v5.variant_qc.isolation_forest.DEFAULT_GATK_IMAGE = 'us.gcr.io/broad-gatk/gatk:4.6.1.0'
GATK image. 4.6.1.0 ships the conda env (scikit-learn, dill) needed for PYTHON_IFOREST.
- gnomad_qc.v5.variant_qc.isolation_forest.extract_variant_annotations_job(b, mode, sites_only_vcf, features, resource_args, calling_intervals_args, exclude_intervals, out_root, gatk_image, gcp_billing_project)[source]
Run GATK ExtractVariantAnnotations for one variant type mode.
- Parameters:
b (
Batch) – Batch to add the job to.mode (
str) – Variant type mode, ‘SNP’ or ‘INDEL’.sites_only_vcf (
str) – AS-annotated sites-only input VCF.features (
List[str]) – Features to extract for this mode.resource_args (
str) – GATK--resourceargs for the labeled sets.calling_intervals_args (
str) – GATK-Largs restricting the extracted sites.exclude_intervals (
str) – Telomere/centromere intervals excluded from the training/calibration set (-XL); these sites are not scored either.out_root (
str) – Output prefix (files written as{out_root}.*).gatk_image (
str) – GATK docker image.gcp_billing_project (
str) – GCP billing project for requester-pays buckets.
- Return type:
Job- Returns:
Job with output ResourceGroup
extract.
- gnomad_qc.v5.variant_qc.isolation_forest.train_variant_annotations_model_job(b, mode, annotations_hdf5, out_root, gatk_image, gcp_billing_project, hyperparameters_json=None)[source]
Run GATK TrainVariantAnnotationsModel (PYTHON_IFOREST) for one variant type mode.
- Parameters:
b (
Batch) – Batch to add the job to.mode (
str) – Variant type mode, ‘SNP’ or ‘INDEL’.annotations_hdf5 (
ResourceFile) – Extracted annotations HDF5 from the extract step.out_root (
str) – Output prefix (files written as{out_root}.{mode}.*).gatk_image (
str) – GATK docker image.gcp_billing_project (
str) – GCP billing project for requester-pays buckets.hyperparameters_json (
Optional[str]) – Optional model hyperparameters JSON path.
- Return type:
Job- Returns:
Job with output ResourceGroup
model.
- gnomad_qc.v5.variant_qc.isolation_forest.split_intervals_job(b, calling_intervals_args, exclude_intervals, scatter_count, gatk_image, gcp_billing_project)[source]
Split the calling contigs into
scatter_countintervals, excluding telomeres/centromeres.- Parameters:
b (
Batch) – Batch to add the job to.calling_intervals_args (
str) – GATK-Largs covering the calling contigs.exclude_intervals (
str) – Path to an intervals file to exclude (-XL).scatter_count (
int) – Number of interval shards to produce.gatk_image (
str) – GATK docker image.gcp_billing_project (
str) – GCP billing project for requester-pays buckets.
- Return type:
Job- Returns:
Job with output ResourceGroup
intervals.
- gnomad_qc.v5.variant_qc.isolation_forest.score_variant_annotations_job(b, mode, sites_only_vcf, features, extracted_vcf, model, resource_args, interval, idx, out_root, gatk_image, gcp_billing_project)[source]
Run GATK ScoreVariantAnnotations for one variant type mode over one interval.
- Parameters:
b (
Batch) – Batch to add the job to.mode (
str) – Variant type mode, ‘SNP’ or ‘INDEL’.sites_only_vcf (
str) – AS-annotated sites-only input VCF.features (
List[str]) – Features for this mode (must match the trained model).extracted_vcf (
ResourceFile) – Extracted training/calibration VCF from the extract step.model (
Union[ResourceGroup,str]) –--model-prefix: the trained model ResourceGroup from the train job, or a prefix string built from a read-in group (see isolation_forest_workflow).resource_args (
str) – GATK--resourceargs for the labeled sets.interval (
ResourceFile) – Interval file to restrict scoring to.idx (
int) – Shard index, used to identify the job in the Batch.out_root (
str) – Output prefix (files written as{out_root}.*).gatk_image (
str) – GATK docker image.gcp_billing_project (
str) – GCP billing project for requester-pays buckets.
- Return type:
Job- Returns:
Job with output ResourceGroup
output_score.
Note
GATK traversal returns every record overlapping the shard interval, and the default output mode (ANYWHERE) writes all of them, so a deletion straddling a shard boundary is emitted by both neighbors.
--variant-output-filtering STARTS_INwrites a record only if its start is in the interval, making emission exactly-once across shards. It also makes-XLstart-based, which is the definition reconcile_scored_sites uses.
- gnomad_qc.v5.variant_qc.isolation_forest.check_run_config(run_prefix, config, overwrite, write)[source]
Guard reuse of outputs under
run_prefixagainst a changed run configuration.Extract/model/score outputs are reused by path alone, so a rerun with different contigs, scatter count, or inputs would silently reuse region- or split-mismatched files. Compare against the config recorded by the first run and raise on any difference unless
overwrite.- Parameters:
run_prefix (
str) – GCS prefix for this run’s GATK outputs.config (
Dict[str,Any]) – Current run configuration.overwrite (
bool) – Whether outputs are being overwritten. Skips the comparison only in the Batch step; –load-only reuses whatever is underrun_prefix(its –overwrite only replaces the result HT), so it always compares.write (
bool) – Whether to recordconfig(the Batch step; the load step only compares).
- Return type:
None- Returns:
None.
- gnomad_qc.v5.variant_qc.isolation_forest.isolation_forest_workflow(b, sites_only_vcf, run_prefix, out_vcf_name, singletons_vcf, calling_intervals_args, exclude_intervals, scatter_count, gatk_image, gcp_billing_project, hyperparameters_json=None, overwrite=False)[source]
Add the per-mode extract/train/score jobs (SNP and INDEL) to the Batch.
Extract and train outputs are reused if they already exist (unless
overwrite).- Parameters:
b (
Batch) – Batch to add jobs to.sites_only_vcf (
str) – AS-annotated sites-only input VCF.run_prefix (
str) – GCS prefix for this run’s GATK outputs.out_vcf_name (
str) – Base name for scored VCF shards.singletons_vcf (
Optional[str]) – Optional true-positive singletons VCF.calling_intervals_args (
str) – GATK-Largs for the calling contigs.exclude_intervals (
str) – Path to the telomere/centromere exclusion intervals file.scatter_count (
int) – Number of interval shards for scoring.gatk_image (
str) – GATK docker image.gcp_billing_project (
str) – GCP billing project for requester-pays buckets.hyperparameters_json (
Optional[str]) – Optional model hyperparameters JSON path.overwrite (
bool) – Whether to rerun extract/train/score even if outputs exist.
- Return type:
None- Returns:
None.
- gnomad_qc.v5.variant_qc.isolation_forest.export_exclusion_intervals(out_path)[source]
Export telomere/centromere intervals to a GATK-readable
chr:start-endfile.Written to a stable path (not a temp file) so the same file used by GATK
-XLcan be read back during reconciliation to define the excluded region exactly.telomeres_and_centromeres is a BED import, so its intervals are half-open
[start, end), while GATK readschr:start-endas 1-based inclusive on both ends. Emitend - 1so the file covers exactly the interval.- Parameters:
out_path (
str) – Destination path for the intervals file.- Return type:
str- Returns:
out_path.
- gnomad_qc.v5.variant_qc.isolation_forest.export_v4_test_sites_vcf(contigs, out_path)[source]
Export the v4 split info HT, filtered to
contigs, as a GATK-readable sites VCF.v4’s
--export-split-info-vcfoutput is misnamed: it is exported from the unsplit info HT, so it cannot stand in for the split v5 info VCF. Read the split info HT directly instead.- Parameters:
contigs (
List[str]) – Contigs to keep.out_path (
str) – Destination VCF path.
- Return type:
None- Returns:
None.
- gnomad_qc.v5.variant_qc.isolation_forest.run_batch(b, description)[source]
Run a Batch to completion and raise if it did not succeed.
Batch.runwaits but only prints the final state, so a failed batch would otherwise fall through to the next step.- Parameters:
b (
Batch) – Batch to run.description (
str) – Batch description used in the error message.
- Return type:
None- Returns:
None.
- gnomad_qc.v5.variant_qc.isolation_forest.merge_iforest_result(run_prefix, out_vcf_name, model_id, scatter_count, n_partitions, header_path, array_elements_required, transmitted_singletons, sibling_singletons, adj)[source]
Merge the per-mode scored VCFs into a single variant QC result HT.
The SNP-pass VCF contributes SNP alleles and the INDEL-pass VCF contributes all non-SNP alleles, so the union covers every scored allele exactly once.
Only the current run’s
scatter_countshards are read (by explicit index), so a rerun with a smaller scatter count does not pick up stale shards from a prior run.- Parameters:
run_prefix (
str) – GCS prefix for this run’s GATK outputs.out_vcf_name (
str) – Base name used for scored VCF shards.model_id (
str) – Isolation forest model ID (starts withif_).scatter_count (
int) – Number of score shards per mode (must match the scoring run).n_partitions (
int) – Number of partitions for the imported HTs.header_path (
Optional[str]) – Optional VCF header file for import.array_elements_required (
bool) – Value passed to hl.import_vcf.transmitted_singletons (
bool) – Whether transmitted singletons were a training set.sibling_singletons (
bool) – Whether sibling singletons were a training set.adj (
bool) – Whether the singletons VCF used adj genotypes.
- Return type:
- Returns:
Merged variant QC result HT.
- gnomad_qc.v5.variant_qc.isolation_forest.reconcile_scored_sites(sites_only_vcf, contigs, exclude_intervals_path, scored_ht, array_elements_required=False)[source]
Raise if any input variant in the scored region is missing from the output or unscored.
GATK writes every input record within its interval to the scored VCF, so any variant absent from the merged output was dropped (e.g. a truncated shard). GATK classifies every allele as either SNP or INDEL, so every variant should also carry a SCORE.
Input and output are both split, so variants are compared on the full key.
- Parameters:
sites_only_vcf (
str) – Input sites VCF that was scored.contigs (
List[str]) – Calling contigs passed to GATK (-L).exclude_intervals_path (
str) – Telomere/centromere exclusion file passed to GATK (-XL).scored_ht (
Table) – Merged scored result HT.array_elements_required (
bool) – Value passed to hl.import_vcf for the input.
- Return type:
None- Returns:
None.