diff --git a/.bumpversion.cfg b/.bumpversion.cfg index d86a437..8e74266 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 0.3.1 +current_version = 0.3.2 commit = True tag = False parse = (?P\d+)\.(?P\d+)\.(?P[a-z0-9+]+) diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index 22018a8..9fd84a5 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -18,7 +18,7 @@ jobs: uses: actions/checkout@v4 - name: Free Disk Space (Ubuntu) - uses: jlumbroso/free-disk-space@main + uses: jlumbroso/free-disk-space@v1.3.0 with: # this might remove tools that are actually needed, # if set to "true" but frees about 6 GB diff --git a/CHANGELOG.md b/CHANGELOG.md index e4a5ea6..0b43c32 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,23 @@ # bolt changelog +## 0.3.2 + +- Fix: `rescue.py`'s SAGE VCF header-consistency check compared the DRAGEN input VCF's headers against the expected SAGE header descriptions instead of the actual SAGE VCF's — the check ran but validated nothing (introduced [24](https://github.com/umccr/bolt/pull/24)). Now checks `sage_vcf_fp`. First-ever test coverage added for `rescue.py` and for `util.check_annotation_headers` +- Fix: hypermutated `report` samples (tiered-selection path) never disabled `--estimate_msi`/`--estimate_tmb` — TMB/MSI were computed on a trimmed variant subset with no indication in the report. Now passes `disable_estimates=True` when `select_pcgr_variants` trims input, matching `run_somatic_chunk`'s existing chunked-input behavior +- Fix: `get_annotations_vcf` crashed with `AssertionError` on a duplicate PCGR VCF key (variant mapping to multiple transcripts) — now logs a warning and keeps the first entry ([34](https://github.com/umccr/bolt/pull/34)) +- Fix: `collect_pcgr_annotation_data` / `collect_cpsr_annotation_data` had the same duplicate-key `AssertionError` on their TSV paths — somatic path now keeps the most actionable tier on a duplicate, germline/CPSR path keeps the first entry; both warn instead of crashing +- Test: `TestGetAnnotationsVcf`, `TestCollectPcgrAnnotationData`, `TestCollectCpsrAnnotationData` — duplicate-key regression coverage for all three annotation-collection paths +- Test: `TestMergeVcfFiles` — bcftools-guarded integration test proving `merge_vcf_files` (`bcftools merge -m all`) losslessly reassembles PCGR sites-only chunk VCFs (no loss/dup, position-sorted, indexed output). Locks the sites-only invariant that keeps `bcftools merge` safe for the hypermutated chunk-merge path; documented with a `NOTE` in `merge_vcf_files` +- [32](https://github.com/umccr/bolt/pull/32) - Fix `PCGR_MUTATION_HOTSPOT=.` (dot placeholder) treated as truthy in retention check — was preventing tiered filtering from running for any sample with >450k PASS variants ([sash#52](https://github.com/umccr/sash/issues/52)) +- [32](https://github.com/umccr/bolt/pull/32) - Graceful PCGR skip when `select_pcgr_variants` cannot cap variants to `MAX_SOMATIC_VARIANTS` — logs warning and continues without cancer report; non-PCGR outputs still publish ([sash#52](https://github.com/umccr/sash/issues/52)) +- Fix: disable `--estimate_msi`/`--estimate_tmb` in chunked PCGR annotation runs — estimates on partial VCFs are not meaningful +- Fix: pin `jlumbroso/free-disk-space` CI action to `v1.3.0` (was `@main`) +- Test: `TestSelectPcgrVariants` — 8 integration tests covering tiered trimming, hotspot retention, PANEL retention, NONCODING-first drop order, and the `PCGR_MUTATION_HOTSPOT=.` regression +- Test: `TestEntrySkipsPcgrOnOverflow` — 2 tests: entry() skips PCGR on `RuntimeError` from unresolvable overflow; entry() calls PCGR normally when within limit +- Test: `TestSelectPcgrVariantsRaisesOnUnresolvableOverflow` — asserts `RuntimeError` when retained variants alone exceed `MAX_SOMATIC_VARIANTS` +- Test: `TestRunSomaticCommandArgs` — 2 tests: `--estimate_signatures` absent from all `run_somatic` commands; `disable_estimates=True` suppresses `--estimate_msi`/`--estimate_tmb` +- Test: `TestRunSomaticChunkArgMapping.test_disable_estimates_passed_to_run_somatic` — `run_somatic_chunk` passes `disable_estimates=True` to every `run_somatic` call + ## 0.3.1 - Fix `ModuleNotFoundError: No module named 'pkg_resources'` in `bolt:0.3.0-multiqc` — add `setuptools <81` to conda env @@ -7,18 +25,27 @@ - Fix VCF writers not closed in `transfer_annotations_somatic` and `transfer_annotations_germline` — BGZip output could be truncated - Fix `split_vcf` writing uncompressed plain `.vcf` chunks — now uses `.vcf.gz` with `wz` mode - Fix `PCGR_ACTIONABILITY_TIER` VCF header description — updated to match stored short-form values (`1`,`2`,`3`,`4`,`N`) -- Add regression test for chunk file compression (`test_chunks_are_gzipped`) +- Fix `split_vcf` chunks not tabix-indexed, causing PCGR to fail reading them +- Fix `build.yaml` and `Dockerfile.pcgr` build issues +- Remove unused `logging` import and fix `PCGR_MAX_SOMATIC_VARIANTS` header description in `constants.py` +- Bump `r-gpgr` to 2.3.1 in `Dockerfile.gpgr` +- Add CI smoke tests to catch Docker image startup failures before push +- [31](https://github.com/umccr/bolt/pull/31) - Drop `--estimate_signatures` from PCGR somatic invocation — signature analysis comes from gpgr/sigrap downstream; keep `--estimate_msi`/`--estimate_tmb` ([sash#57](https://github.com/umccr/sash/issues/57)) +- Test: regression test for chunk file compression (`test_chunks_are_gzipped`) +- Test: regression test for chunk tabix indexing (`test_chunks_are_tabix_indexed`) +- Test: `TestTierOrdering` — 3 tests verifying `PCGR_TIERS_FILTERING` uses short forms (`N`,`4`,`3`,`2`,`1`) and NONCODING precedes TIER_1 in `get_ordering()` +- Test: `TestSplitVcf` — 4 tests: chunking above/below limit, `.vcf.gz` compression, `.tbi` indexing +- Test: `TestSelectPcgrVariants` — initial 6 integration tests for tiered trimming logic +- Test: `TestGetVariantFilterData`, `TestDetermineFilter`, `TestGetImpacts` — 14 unit tests covering variant attribute extraction and filter-category determination +- Test: `TestCountVariantProcess` — 4 tests: `is_hypermutated` flag, DRAGEN count, SAGE_NOVEL exclusion, annotation-filter exclusion +- Test: `TestRunSomaticChunkArgMapping` — asserts `pcgr_conda` is not shifted into `pcgr_threads` position on positional arg mapping ## 0.3.0 - [28](https://github.com/umccr/bolt/pull/28) - gpgr version bump to 2.2.12 for cancer report hypermutated flag fix - [17](https://github.com/umccr/bolt/pull/17) - change dragen HRD file optional - - [14](https://github.com/umccr/bolt/pull/14) - gpgr version bump to 2.2.0 - - [3](https://github.com/scwatts/bolt/pull/3) - Improve PCGR / CPSR argument handling - - [6](https://github.com/umccr/bolt/pull/6) - Change oncoanalyser v2.0.0 update, with switch sv caller from GRIPSS to eSVee - - [9](https://github.com/umccr/bolt/pull/9) Add hypermutation sample handling diff --git a/TESTING.md b/TESTING.md new file mode 100644 index 0000000..1676ea2 --- /dev/null +++ b/TESTING.md @@ -0,0 +1,74 @@ +# Testing + +## Running tests + +Canonical (matches CI): + +```bash +python -m unittest discover tests/ --buffer +``` + +Also supported: + +```bash +python -m pytest tests/ -v +``` + +Single test: + +```bash +python -m pytest tests/test_smlv_somatic_filter.py::TestSmlvSomaticFilter::test_min_af_filter -v +``` + +Tests in `tests/` are pure Python + in-memory `cyvcf2` logic and must run +without any bioinformatics binary installed (no `bcftools`, `pcgr`, `cpsr`, +`vcfanno`, `snpEff`, `gpgr`, VEP). Functions that shell out to those tools are +either mocked/patched in tests or are not unit-tested (see below). The one +exception is `TestMergeVcfFiles`, an integration test that exercises the real +`merge_vcf_files` → `bcftools merge` path; it is guarded with +`@unittest.skipUnless(shutil.which('bcftools'), ...)`, so it runs in the conda +CI env and skips cleanly (never fails) where `bcftools` is absent. + +## Test Coverage + +| Module | Status | Functions covered | Test file | +|---|---|---|---| +| `bolt/util.py` | Partial | `get_vcf_header_entry`, `get_vcf_header_line`, `get_qualified_vcf_annotation`, `add_vcf_header_entry`, `merge_tsv_files`, `merge_vcf_files` (bcftools-guarded lossless/sorted integration test), `check_annotation_headers` | `tests/test_util.py` | +| `bolt/common/pcgr.py` | Partial | `get_ordering`, `get_impacts`, `determine_filter`, `get_variant_filter_data`, `split_vcf`, `run_somatic_chunk` (arg-mapping regression) | `tests/test_pcgr.py` | +| `bolt/common/pcgr.py` | Partial | `parse_genomic_change`, `get_impacts_higher`, `get_annotation_entry_tsv`, `compile_annotation_data`, `annotate_record`, `get_annotations_vcf` (duplicate-key regression), `collect_pcgr_annotation_data` (duplicate-key tier resolution), `collect_cpsr_annotation_data` (duplicate-key regression) | `tests/test_pcgr_annotation.py` | +| `bolt/workflows/smlv_somatic/filter.py` | Partial | `set_filter_data` | `tests/test_smlv_somatic_filter.py` | +| `bolt/workflows/smlv_somatic/report.py` | Partial | `select_pcgr_variants`, `count_variant_process`, `entry` overflow handling, `entry` `disable_estimates` branching | `tests/test_smlv_somatic_report.py` | +| `bolt/workflows/smlv_somatic/rescue.py` | Partial | `annotate_existing_sage_calls` (SAGE VCF header-consistency check only) | `tests/test_smlv_somatic_rescue.py` | + +"Partial" means the module has meaningful test coverage for its pure/testable logic, +but not every function in the file is tested (see below for what is excluded and why). + +## Untested / not unit-testable in CI + +These require live bioinformatics binaries (bcftools, PCGR/CPSR, vcfanno, snpEff, +gpgr) or full end-to-end subprocess orchestration, so they are excluded from the +unit test suite: + +| Function/module | Reason | +|---|---| +| `bolt/util.py: count_vcf_records` | Shells out to `bcftools view` | +| `bolt/util.py: execute_command` | Spawns real subprocesses via `/bin/bash` | +| `bolt/common/pcgr.py: prepare_vcf_somatic` / `prepare_vcf_germline` | Shells out to `bcftools index`/`bcftools view`/`bcftools annotate` | +| `bolt/common/pcgr.py: run_somatic` / `run_somatic_chunk` (execution path) / `run_germline` | Invoke `pcgr`/`cpsr` CLI directly | +| `bolt/common/pcgr.py: transfer_annotations_somatic` / `transfer_annotations_germline` | Depend on real PCGR/CPSR TSV+VCF output files | +| `bolt/common/pcgr.py: merging_pcgr_files` | Wraps `merge_vcf_files`/`merge_tsv_files` (bcftools-dependent) | +| `bolt/workflows/smlv_somatic/annotate.py` | Orchestrates vcfanno + PON + PCGR subprocess pipeline; no unit tests | +| `bolt/workflows/smlv_somatic/rescue.py` (all functions except the header check) | SAGE hotspot recall orchestrates `bcftools isec`/`concat`/`annotate` subprocesses end-to-end; no unit tests | +| `bolt/workflows/smlv_somatic/prepare.py` | bcftools-based VCF prep; no unit tests | +| `bolt/workflows/smlv_germline/prepare.py` | bcftools-based panel region selection; no unit tests | +| `bolt/workflows/smlv_germline/report.py` | bcftools stats + CPSR report generation; no unit tests | +| `bolt/workflows/sv_somatic/annotate.py` | snpEff subprocess annotation; no unit tests | +| `bolt/workflows/sv_somatic/prioritise.py` | Wraps `external/prioritize_sv.py`; no unit tests | +| `bolt/workflows/other/cancer_report.py` | Invokes `gpgr.R canrep` subprocess; no unit tests | +| `bolt/workflows/other/multiqc_report.py` | Invokes `multiqc` subprocess; no unit tests | +| `bolt/workflows/other/purple_baf_plot.py` | Invokes `circos` subprocess; no unit tests | +| `bolt/external/prioritize_sv.py` | Vendored third-party AstraZeneca SV annotation logic; no unit tests | + +Contributions raising coverage for these are welcome, but will require mocking +subprocess calls (as done in `tests/test_smlv_somatic_report.py` for +`select_pcgr_variants`) or dedicated fixture VCFs/TSVs under `tests/fixtures/`. diff --git a/bolt/common/pcgr.py b/bolt/common/pcgr.py index f1ce2cb..0726a27 100644 --- a/bolt/common/pcgr.py +++ b/bolt/common/pcgr.py @@ -146,7 +146,8 @@ def run_somatic(input_fp, pcgr_refdata_dir, vep_dir, output_dir, chunk_nbr=None, f'--control_af_tag NORMAL_AF', f'--genome_assembly grch38', f'--assay WGS', - *([] if disable_estimates else ['--estimate_signatures', '--estimate_msi', '--estimate_tmb']), + # NOTE: --estimate_signatures dropped (umccr/sash#57); keep msi/tmb. + *([] if disable_estimates else ['--estimate_msi', '--estimate_tmb']), f'--vcfanno_n_proc {vcfanno_threads}', f'--vep_n_forks {vep_forks}', f'--vep_pick_order biotype,rank,appris,tsl,ccds,canonical,length,mane_plus_clinical,mane_select', @@ -352,6 +353,25 @@ def transfer_annotations_germline(input_fp, normal_name, cpsr_dir, output_dir): output_fh.close() +# Lower value = more clinically actionable. Used to resolve duplicate PCGR TSV +# rows for the same variant (PCGR can emit one row per transcript mapping). +_TIER_ORDER = {'1': 0, '2': 1, '3': 2, '4': 3, 'N': 4} + + +def _normalise_tier(raw_tier): + # Normalize PCGR actionability tier to simple values: '1','2','3','4','N' + tier_norm = (raw_tier or '').strip().replace('_', ' ').upper() + if tier_norm in ('TIER 1', 'TIER1', '1'): + return '1' + elif tier_norm in ('TIER 2', 'TIER2', '2'): + return '2' + elif tier_norm in ('TIER 3', 'TIER3', '3'): + return '3' + elif tier_norm in ('TIER 4', 'TIER4', '4'): + return '4' + return 'N' + + def collect_pcgr_annotation_data(tsv_fp, vcf_fp, info_field_map): # Gather all annotations from TSV data_tsv = dict() @@ -360,23 +380,19 @@ def collect_pcgr_annotation_data(tsv_fp, vcf_fp, info_field_map): with open_fn(tsv_fp, 'rt') as tsv_fh: for record in csv.DictReader(tsv_fh, delimiter='\t'): key, record_ann = get_annotation_entry_tsv(record, info_field_map) - assert key not in data_tsv - - # Normalize PCGR actionability tier to simple values: '1','2','3','4','N' - raw_tier = (record.get('ACTIONABILITY_TIER') or '').strip() - tier_norm = raw_tier.replace('_', ' ').upper() - if tier_norm in ('TIER 1','TIER1','1'): - tier_val = '1' - elif tier_norm in ('TIER 2','TIER2','2'): - tier_val = '2' - elif tier_norm in ('TIER 3','TIER3','3'): - tier_val = '3' - elif tier_norm in ('TIER 4','TIER4','4'): - tier_val = '4' - else: - tier_val = 'N' + tier_val = _normalise_tier(record.get('ACTIONABILITY_TIER')) record_ann[constants.VcfInfo.PCGR_ACTIONABILITY_TIER] = tier_val + # NOTE(QC): PCGR can emit multiple TSV rows for the same variant when it + # maps to multiple transcripts. Keep the most actionable entry (lowest + # tier number) rather than asserting. See https://github.com/umccr/bolt/pull/34 + if key in data_tsv: + existing_tier = data_tsv[key][constants.VcfInfo.PCGR_ACTIONABILITY_TIER] + if _TIER_ORDER[tier_val] >= _TIER_ORDER[existing_tier]: + logger.warning(f'Duplicate PCGR TSV key {key}: keeping tier {existing_tier}, skipping tier {tier_val}') + continue + logger.warning(f'Duplicate PCGR TSV key {key}: replacing tier {existing_tier} with more actionable tier {tier_val}') + # Store annotation data data_tsv[key] = record_ann @@ -404,7 +420,15 @@ def collect_cpsr_annotation_data(tsv_fp, vcf_fp, info_field_map): record['ALT'] = re_result.group('alt') key, record_ann = get_annotation_entry_tsv(record, info_field_map) - assert key not in data_tsv + + # NOTE(QC): CPSR can emit multiple TSV rows for the same variant when it + # maps to multiple transcripts (same class of issue as the PCGR somatic + # TSV path above). CPSR rows carry no actionability tier to break ties + # on, so keep the first entry and warn. See https://github.com/umccr/bolt/pull/34 + if key in data_tsv: + logger.warning(f'Duplicate CPSR TSV key {key}: keeping first entry') + continue + data_tsv[key] = record_ann # Gather annotations from VCF @@ -414,18 +438,14 @@ def collect_cpsr_annotation_data(tsv_fp, vcf_fp, info_field_map): return compile_annotation_data(data_tsv, data_vcf) def parse_genomic_change(genomic_change): - """ - Parse a genomic change string, e.g., "3:g.41224645T>C" - Returns a tuple: (chrom, pos, ref, alt) - """ - # Regular expression for the format "chrom:g.posRef>Alt" + # Format: "chrom:g.posRef>Alt" e.g. "3:g.41224645T>C" pattern = r'^(?P\w+):g\.(?P\d+)(?P\w+)>(?P\w+)$' match = re.match(pattern, genomic_change) if not match: raise ValueError(f"Format not recognized: {genomic_change}") - - # Get values and format as needed - chrom = f"chr{match.group('chrom')}" + chrom = match.group('chrom') + if not chrom.startswith('chr'): + chrom = f'chr{chrom}' pos = int(match.group('pos')) ref = match.group('ref') alt = match.group('alt') @@ -439,7 +459,14 @@ def get_annotations_vcf(vcf_fp, info_field_map): assert len(record.ALT) == 1 [alt] = record.ALT key = (f'chr{record.CHROM}', record.POS, record.REF, alt) - assert key not in data_vcf + + # NOTE(QC): PCGR can emit duplicate variant entries in its output VCF when a + # variant maps to multiple transcripts. Keep the first entry and warn rather + # than asserting, consistent with the TSV-side dedup. + # https://github.com/umccr/bolt/pull/34 + if key in data_vcf: + logger.warning(f'Duplicate PCGR VCF key {key}: keeping first entry') + continue data_vcf[key] = dict() for info_dst, info_src in info_field_map.items(): @@ -505,7 +532,7 @@ def annotate_record(record, annotations, *, allow_missing=False): if allow_missing: return record else: - assert key not in annotations + assert False, f'Missing annotation key: {key}' # Transfer annotations for info_enum, v in annotations[key].items(): @@ -514,11 +541,6 @@ def annotate_record(record, annotations, *, allow_missing=False): return record def split_vcf(input_vcf, output_dir, *, max_variants=None): - """ - Splits a VCF file into multiple chunks, each containing up to max_variants variants. - Each chunk includes the VCF header. - Ensures no overlapping positions between chunks. - """ if max_variants is None: max_variants = constants.MAX_SOMATIC_VARIANTS elif max_variants <= 0: @@ -566,27 +588,25 @@ def split_vcf(input_vcf, output_dir, *, max_variants=None): def run_somatic_chunk(vcf_chunks, pcgr_data_dir, vep_dir, output_dir, pcgr_output_dir, max_threads, pcgr_conda, pcgrr_conda): pcgr_tsv_files = [] pcgr_vcf_files = [] - - # Process each chunk sequentially + for chunk_number, vcf_file in enumerate(vcf_chunks, start=1): - pcgr_tsv_fp, pcgr_vcf_fp = run_somatic(vcf_file, pcgr_data_dir, vep_dir, pcgr_output_dir, chunk_nbr=chunk_number, threads=max_threads, pcgr_conda=pcgr_conda, pcgrr_conda=pcgrr_conda) + pcgr_tsv_fp, pcgr_vcf_fp = run_somatic(vcf_file, pcgr_data_dir, vep_dir, pcgr_output_dir, chunk_nbr=chunk_number, threads=max_threads, pcgr_conda=pcgr_conda, pcgrr_conda=pcgrr_conda, disable_estimates=True) if pcgr_tsv_fp: pcgr_tsv_files.append(pcgr_tsv_fp) if pcgr_vcf_fp: pcgr_vcf_files.append(pcgr_vcf_fp) - + merged_vcf_fp, merged_tsv_fp = merging_pcgr_files(output_dir, pcgr_vcf_files, pcgr_tsv_files) return merged_tsv_fp, merged_vcf_fp + def merging_pcgr_files(output_dir, pcgr_vcf_files, pcgr_tsv_files): pcgr_dir = pathlib.Path(output_dir) / 'pcgr' pcgr_dir.mkdir(exist_ok=True) - # Merge all TSV files into a single file in the pcgr directory merged_tsv_fp = pcgr_dir / "nosampleset.pcgr_acmg.grch38.snvs_indels.tiers.tsv.gz" util.merge_tsv_files(pcgr_tsv_files, merged_tsv_fp) - # Step 5: Merge all VCF files into a single file in the pcgr directory merged_vcf_path = pcgr_dir / "nosampleset.pcgr.grch38.pass" if len(pcgr_vcf_files) == 1: # NOTE(QC): bcftools merge requires 2+ inputs; with a single chunk there is diff --git a/bolt/logging_config.py b/bolt/logging_config.py index 0ddda2b..66b80ee 100644 --- a/bolt/logging_config.py +++ b/bolt/logging_config.py @@ -28,7 +28,8 @@ def setup_logging(output_dir, script_name): logging.basicConfig( level=logging.DEBUG, format='%(asctime)s - %(name)s - %(levelname)s - %(message)s', - handlers=[file_handler, console_handler] + handlers=[file_handler, console_handler], + force=True, ) logger = logging.getLogger(__name__) logger.info("Logging setup complete") \ No newline at end of file diff --git a/bolt/util.py b/bolt/util.py index c3b07fe..6baf634 100644 --- a/bolt/util.py +++ b/bolt/util.py @@ -23,7 +23,6 @@ def get_project_root(): def execute_command(command, log_file_path=None): - # set -e: exit on error, -u: exit on unset variable, -o pipefail: pipeline fails if any command fails prepared_command = f'set -euo pipefail; {textwrap.dedent(command)}' logger.info("Executing command: %s", command.strip()) @@ -129,15 +128,10 @@ def get_qualified_vcf_annotation(anno_enum): return f'{anno_enum.namespace}/{anno_enum.value}' def merge_tsv_files(tsv_files, merged_tsv_fp): - """ - Merge gzipped TSV files into a single gzipped TSV. - """ - with gzip.open(merged_tsv_fp, 'wt', encoding='utf-8') as merged_tsv: for i, tsv_file in enumerate(tsv_files): with gzip.open(tsv_file, 'rt', encoding='utf-8') as infile: for line_number, line in enumerate(infile): - # Skip header except for the first file if i > 0 and line_number == 0: continue merged_tsv.write(line) @@ -145,43 +139,32 @@ def merge_tsv_files(tsv_files, merged_tsv_fp): def merge_vcf_files(vcf_files, merged_vcf_fp): - """ - Merges multiple VCF files into a single sorted VCF file using bcftools. - - Parameters: - - vcf_files: List of paths to VCF files to be merged. - - merged_vcf_fp: Path to the output merged VCF file (without extension). - - Returns: - - Path to the sorted merged VCF file. - """ + # NOTE(QC): reassembles PCGR hypermutated chunk outputs. `bcftools merge` is + # a multi-sample joiner and would fail with "Duplicate sample names" on + # same-named inputs that carry a genotype column — but this is safe here + # because PCGR inputs/outputs are sites-only (pcgr.prepare_vcf_somatic / + # get_minimal_header strip all FORMAT and sample columns), so merge just + # produces the union of disjoint chunks. Do not switch to a genotype-bearing + # VCF upstream without revisiting this. Locked by tests/test_util.py + # TestMergeVcfFiles (bcftools-guarded lossless/sorted regression test). merged_vcf_fp = pathlib.Path(merged_vcf_fp) merged_unsorted_vcf = merged_vcf_fp.parent / f'{merged_vcf_fp.name}.unsorted.vcf.gz' merged_vcf = merged_vcf_fp.parent / f'{merged_vcf_fp.name}.vcf.gz' - # Prepare the bcftools merge command arguments + delimiter_padding = ' ' * 10 + delimiter = f' \\\n{delimiter_padding}' + command_args = [ 'bcftools merge', '-m all', '-Oz', f'-o {merged_unsorted_vcf}', ] + [str(vcf_file) for vcf_file in vcf_files] - - # Format the command for readability - delimiter_padding = ' ' * 10 - delimiter = f' \\\n{delimiter_padding}' command_args_str = delimiter.join(command_args) - - command = f''' + execute_command(f''' {command_args_str} - ''' - - # Run the bcftools merge command - logger.info("Running bcftools merge...") - execute_command(command) - logger.info(f"Merged VCF written to: {merged_unsorted_vcf}") + ''') - # Sort the merged VCF file sort_command_args = [ 'bcftools sort', '-Oz', @@ -189,30 +172,12 @@ def merge_vcf_files(vcf_files, merged_vcf_fp): f'{merged_unsorted_vcf}' ] sort_command_args_str = delimiter.join(sort_command_args) - sort_command = f''' + execute_command(f''' {sort_command_args_str} - ''' - - logger.info("Sorting merged VCF file...") - execute_command(sort_command) - logger.info(f"Sorted merged VCF written to: {merged_vcf}") - - # Index the sorted merged VCF file - index_command_args = [ - 'bcftools index', - '-t', - f'{merged_vcf}' - ] - index_command_args_str = delimiter.join(index_command_args) - index_command = f''' - {index_command_args_str} - ''' + ''') - logger.info("Indexing sorted merged VCF file...") - execute_command(index_command) - logger.info(f"Indexed merged VCF file: {merged_vcf}.tbi") + execute_command(f'bcftools index -t {merged_vcf}') - # Optionally, remove the unsorted merged VCF file if merged_unsorted_vcf.exists(): merged_unsorted_vcf.unlink() diff --git a/bolt/workflows/smlv_somatic/annotate.py b/bolt/workflows/smlv_somatic/annotate.py index 01d28b2..78c07df 100644 --- a/bolt/workflows/smlv_somatic/annotate.py +++ b/bolt/workflows/smlv_somatic/annotate.py @@ -80,12 +80,11 @@ def entry(ctx, **kwargs): # - Set tumor and normal AF and DP in INFO for PCGR and remove all other annotations # - Run PCGR on minimal VCF (pcgr_prep_fp) # - Transfer selected PCGR annotations to unfiltered VCF (selected_fp) - # - PCGR ACMG TIER [INFO/PCGR_TIER] - # - VEP consequence [INFO/PCR_CSQ] + # - PCGR actionability tier [INFO/PCGR_ACTIONABILITY_TIER] + # - VEP consequence [INFO/PCGR_CSQ] # - Known mutation hotspot [INFO/PCGR_MUTATION_HOTSPOT] - # - ClinVar clinical significant [INFO/PCGR_CLNSIG] + # - ClinVar classification [INFO/PCGR_CLINVAR_CLASSIFICATION] # - Hits in TCGA [INFO/PCGR_TCGA_PANCANCER_COUNT] - # - Hits in PCAWG [INFO/PCGR_ICGC_PCAWG_COUNT] # Prepare VCF for PCGR annotation pcgr_prep_fp = pcgr.prepare_vcf_somatic( @@ -97,7 +96,7 @@ def entry(ctx, **kwargs): pcgr_output_dir = output_dir / 'pcgr' total_variants = util.count_vcf_records(pcgr_prep_fp) - print(f"Total number of variants in the input VCF: {total_variants}") + logger.info(f"Total number of variants in the input VCF: {total_variants}") # Run PCGR in chunks if exceeding the maximum allowed for somatic variants chunk_size = kwargs.get('pcgr_variant_chunk_size') diff --git a/bolt/workflows/smlv_somatic/report.py b/bolt/workflows/smlv_somatic/report.py index c37f8c1..ccda347 100644 --- a/bolt/workflows/smlv_somatic/report.py +++ b/bolt/workflows/smlv_somatic/report.py @@ -115,36 +115,50 @@ def entry(ctx, **kwargs): # PCGR report purple_data = parse_purple_purity_file(kwargs['purple_purity_fp']) - if variant_counts_process['filter_pass'] <= constants.MAX_SOMATIC_VARIANTS: + pcgr_skipped = False + pcgr_input_trimmed = variant_counts_process['filter_pass'] > constants.MAX_SOMATIC_VARIANTS + if not pcgr_input_trimmed: pcgr_input_vcf_fp = kwargs['vcf_fp'] else: - pcgr_input_vcf_fp = select_pcgr_variants( - kwargs['vcf_fp'], - kwargs['cancer_genes_fp'], + try: + pcgr_input_vcf_fp = select_pcgr_variants( + kwargs['vcf_fp'], + kwargs['cancer_genes_fp'], + kwargs['tumor_name'], + output_dir, + ) + except RuntimeError as e: + # NOTE(QC): tiered filtering could not bring PASS count below + # MAX_SOMATIC_VARIANTS (sash #52). Skip PCGR; sash marks the + # PCGR emits as optional so downstream reports still publish. + logger.warning(f'Skipping PCGR for {kwargs["tumor_name"]}: {e}') + pcgr_skipped = True + + if not pcgr_skipped: + pcgr_prep_fp = pcgr.prepare_vcf_somatic( + pcgr_input_vcf_fp, kwargs['tumor_name'], + kwargs['normal_name'], output_dir, ) - pcgr_prep_fp = pcgr.prepare_vcf_somatic( - pcgr_input_vcf_fp, - kwargs['tumor_name'], - kwargs['normal_name'], - output_dir, - ) - - pcgr_output_dir = output_dir / 'pcgr' - pcgr.run_somatic( - pcgr_prep_fp, - kwargs['pcgr_data_dir'], - kwargs['vep_dir'], - pcgr_output_dir, - threads=kwargs['threads'], - pcgr_conda=kwargs['pcgr_conda'], - pcgrr_conda=kwargs['pcgrr_conda'], - purity=purple_data['purity'], - ploidy=purple_data['ploidy'], - sample_id=kwargs['tumor_name'], - ) + pcgr_output_dir = output_dir / 'pcgr' + pcgr.run_somatic( + pcgr_prep_fp, + kwargs['pcgr_data_dir'], + kwargs['vep_dir'], + pcgr_output_dir, + threads=kwargs['threads'], + pcgr_conda=kwargs['pcgr_conda'], + pcgrr_conda=kwargs['pcgrr_conda'], + purity=purple_data['purity'], + ploidy=purple_data['ploidy'], + sample_id=kwargs['tumor_name'], + # NOTE(QC): pcgr_input_vcf_fp is a tiered-selection subset when + # pcgr_input_trimmed, so TMB/MSI estimates on it are not meaningful — + # same reasoning as run_somatic_chunk's chunked inputs. + disable_estimates=pcgr_input_trimmed, + ) def bcftools_stats_prepare(input_fp, tumor_name, output_dir): @@ -336,7 +350,8 @@ def select_pcgr_variants(vcf_fp, cancer_genes_fp, tumor_name, output_dir): for variant_count, variant in enumerate(cyvcf2.VCF(fp_annotated_out), 1): variant_repr = pcgr.get_variant_repr(variant) - if any(variant.INFO.get(e) for e in constants.RETAIN_FIELDS_FILTERING): + # NOTE(QC): exclude '.' — PCGR writes it as a missing-value placeholder for String fields; cyvcf2 returns it truthy (sash #52). + if any(variant.INFO.get(e) not in (None, '.') for e in constants.RETAIN_FIELDS_FILTERING): continue data = pcgr.get_variant_filter_data(variant) diff --git a/bolt/workflows/smlv_somatic/rescue.py b/bolt/workflows/smlv_somatic/rescue.py index ff69e6e..e544f72 100644 --- a/bolt/workflows/smlv_somatic/rescue.py +++ b/bolt/workflows/smlv_somatic/rescue.py @@ -125,7 +125,7 @@ def annotate_existing_sage_calls(input_fp, tumor_name, sage_vcf_fp, output_dir): constants.VcfFilter.SAGE_LOWCONF: 'SAGE_LOWCONF', } - util.check_annotation_headers(info_field_map_sage, input_fp) + util.check_annotation_headers(info_field_map_sage, sage_vcf_fp) # Add header entries so that they are included in the output file via templating done below util.add_vcf_header_entry(input_fh, constants.VcfFilter.SAGE_LOWCONF) diff --git a/pyproject.toml b/pyproject.toml index 6dd2fba..65c0a2f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,7 +8,7 @@ include = ["bolt*"] [project] name = "bolt" -version = "0.3.1" +version = "0.3.2" authors = [ {name = "Stephen Watts", email = "stephen.watts@umccr.org"}, ] diff --git a/tests/helpers.py b/tests/helpers.py new file mode 100644 index 0000000..1ba472d --- /dev/null +++ b/tests/helpers.py @@ -0,0 +1,46 @@ +"""Shared VCF-construction fixtures used across test modules.""" +import pathlib +import tempfile + +import cyvcf2 + + +# Minimal CSQ: only tokens[1] (consequence) is read by get_impacts() +def _csq(consequence): + return f'A|{consequence}|.|.|.|.|.|.|.|.|.|.|.|.|.|.|.|.|.|.|.' + + +# Minimal VCF header with all INFO fields used by select_pcgr_variants +HEADER = ( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' +) + + +def _write_vcf(path, variants): + with open(path, 'w') as fh: + fh.write(HEADER) + for pos, info in variants: + fh.write(f'chr1\t{pos}\t.\tA\tT\t.\tPASS\t{info}\n') + + +def _count_vcf(fp): + return sum(1 for _ in cyvcf2.VCF(str(fp))) + + +def _make_variant(info_str): + """Return a cyvcf2 Variant built from info_str using the test VCF header.""" + with tempfile.TemporaryDirectory() as tmp: + vcf_path = pathlib.Path(tmp) / 'test.vcf' + _write_vcf(vcf_path, [(100, info_str)]) + return list(cyvcf2.VCF(str(vcf_path)))[0] diff --git a/tests/test_pcgr.py b/tests/test_pcgr.py new file mode 100644 index 0000000..418e1b3 --- /dev/null +++ b/tests/test_pcgr.py @@ -0,0 +1,428 @@ +"""Tests for bolt/common/pcgr.py — tier ordering, filter categorisation, chunking.""" +import gzip +import pathlib +import tempfile +import unittest +from unittest.mock import patch + +import bolt.common.constants as constants +import bolt.common.pcgr as pcgr +import bolt.util as util + +from tests.helpers import _csq, _count_vcf, _make_variant, _write_vcf + + +class TestTierOrdering(unittest.TestCase): + """Verify the PCGR_TIERS_FILTERING fix: values and priority order.""" + + def test_noncoding_filtered_before_tier1(self): + """N (NONCODING) entries must all precede '1' (TIER_1) entries in get_ordering().""" + ordering = pcgr.get_ordering() + tiers = [key[0] for key in ordering] + n_idx = [i for i, t in enumerate(tiers) if t == 'N'] + t1_idx = [i for i, t in enumerate(tiers) if t == '1'] + self.assertTrue(n_idx, 'No NONCODING (N) entries in get_ordering()') + self.assertTrue(t1_idx, 'No TIER_1 (1) entries in get_ordering()') + self.assertLess(max(n_idx), min(t1_idx), + 'All NONCODING entries must precede all TIER_1 entries') + + def test_no_long_form_tier_values(self): + """PCGR_TIERS_FILTERING must use short forms ('1'-'4', 'N'), not 'TIER_1' etc.""" + for v in constants.PCGR_TIERS_FILTERING: + self.assertNotIn('TIER_', v, + f"Found long-form tier value '{v}' — must be short form") + + def test_priority_order(self): + """Full ordering: N before 4 before 3 before 2 before 1.""" + expected = ('N', '4', '3', '2', '1') + self.assertEqual(constants.PCGR_TIERS_FILTERING, expected) + + +class TestGetImpacts(unittest.TestCase): + """Unit tests for pcgr.get_impacts() — CSQ string parsing.""" + + def test_single_consequence(self): + csq = _csq('intron_variant') + self.assertEqual(pcgr.get_impacts(csq), {'intron_variant'}) + + def test_multi_consequences_ampersand(self): + """A single CSQ entry with two consequences joined by & returns both.""" + csq = _csq('intron_variant&upstream_gene_variant') + self.assertEqual(pcgr.get_impacts(csq), {'intron_variant', 'upstream_gene_variant'}) + + def test_multiple_csq_entries_union(self): + """Comma-separated CSQ entries — returns the union of all consequences.""" + csq = f'{_csq("intron_variant")},{_csq("intergenic_variant")}' + self.assertEqual(pcgr.get_impacts(csq), {'intron_variant', 'intergenic_variant'}) + + +class TestDetermineFilter(unittest.TestCase): + """Unit tests for pcgr.determine_filter() — filter category determination.""" + + def _data(self, **overrides): + base = { + 'tier': None, + 'difficult': False, + 'giab_conf': False, + 'intergenic': None, + 'intronic': None, + 'downstream': None, + 'upstream': None, + 'impacts_other': None, + } + base.update(overrides) + return base + + def test_intergenic_difficult(self): + data = self._data(intergenic=True, difficult=True) + self.assertEqual(pcgr.determine_filter(data), ('intergenic', 'difficult')) + + def test_intergenic_no_region(self): + data = self._data(intergenic=True, difficult=False, giab_conf=False) + self.assertEqual(pcgr.determine_filter(data), ('intergenic', 'none')) + + def test_intergenic_giab_conf(self): + data = self._data(intergenic=True, giab_conf=True) + self.assertEqual(pcgr.determine_filter(data), ('intergenic', 'giab_conf')) + + def test_intronic_supersedes_intergenic(self): + """When both intergenic and intronic are present, intronic wins (higher priority).""" + data = self._data(intergenic=True, intronic=True, difficult=True) + self.assertEqual(pcgr.determine_filter(data), ('intronic', 'difficult')) + + def test_impacts_other_highest_priority(self): + """impacts_other is the last to be filtered — it wins over all other impacts.""" + data = self._data( + intergenic=True, intronic=True, downstream=True, + upstream=True, impacts_other=True, difficult=True, + ) + self.assertEqual(pcgr.determine_filter(data), ('impacts_other', 'difficult')) + + def test_no_impact_returns_false(self): + """A variant with no recognisable impact cannot be categorised.""" + data = self._data() # all impacts None + self.assertFalse(pcgr.determine_filter(data)) + + def test_giab_conf_region(self): + data = self._data(impacts_other=True, giab_conf=True) + self.assertEqual(pcgr.determine_filter(data), ('impacts_other', 'giab_conf')) + + +class TestGetVariantFilterData(unittest.TestCase): + """Unit tests for pcgr.get_variant_filter_data() — data extraction from VCF records.""" + + def test_tier_extracted(self): + info = f'PCGR_ACTIONABILITY_TIER=2;PCGR_CSQ={_csq("intron_variant")}' + data = pcgr.get_variant_filter_data(_make_variant(info)) + self.assertEqual(data['tier'], '2') + + def test_intergenic_impact(self): + info = f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}' + data = pcgr.get_variant_filter_data(_make_variant(info)) + self.assertTrue(data['intergenic']) + self.assertFalse(data['intronic']) + self.assertFalse(data['downstream']) + self.assertFalse(data['upstream']) + self.assertFalse(data['impacts_other']) + + def test_intronic_impact(self): + info = f'PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}' + data = pcgr.get_variant_filter_data(_make_variant(info)) + self.assertTrue(data['intronic']) + self.assertFalse(data['intergenic']) + + def test_giab_conf_overrides_difficult(self): + """GIAB_CONF flag must clear the difficult flag even when DIFFICULT_* is also present.""" + info = f'GIAB_CONF;DIFFICULT_segdup;PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}' + data = pcgr.get_variant_filter_data(_make_variant(info)) + self.assertTrue(data['giab_conf']) + self.assertFalse(data['difficult']) + + def test_difficult_without_giab(self): + info = f'DIFFICULT_segdup;PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}' + data = pcgr.get_variant_filter_data(_make_variant(info)) + self.assertTrue(data['difficult']) + self.assertFalse(data['giab_conf']) + + +class TestSplitVcf(unittest.TestCase): + """Tests for pcgr.split_vcf() — chunking the annotation path for large VCFs. + + split_vcf() is the annotate-path strategy for hypermutated samples: it divides + a VCF into ≤MAX_SOMATIC_VARIANTS chunks so each chunk can be run through PCGR + independently. Tested 2026-05-13 with a synthetic 550k VCF: 550k → 450k + 100k. + """ + + def test_chunks_above_limit(self): + """VCF exceeding the limit is split into correctly-sized chunks.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'input.vcf' + # 25 variants, limit=10 → expect 3 chunks (10, 10, 5) + v = [(i * 10, f'PCGR_CSQ={_csq("intron_variant")}') for i in range(1, 26)] + _write_vcf(vcf_fp, v) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 10): + chunks = pcgr.split_vcf(vcf_fp, tmp_path) + + self.assertEqual(len(chunks), 3) + counts = [_count_vcf(c) for c in chunks] + self.assertLessEqual(max(counts), 10) + self.assertEqual(sum(counts), 25) + + def test_no_chunking_within_limit(self): + """VCF within the limit produces a single chunk containing all variants.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'input.vcf' + v = [(i * 10, f'PCGR_CSQ={_csq("intron_variant")}') for i in range(1, 6)] + _write_vcf(vcf_fp, v) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 10): + chunks = pcgr.split_vcf(vcf_fp, tmp_path) + + self.assertEqual(len(chunks), 1) + self.assertEqual(_count_vcf(chunks[0]), 5) + + def test_chunks_are_gzipped(self): + """Chunk files must be .vcf.gz — plain .vcf chunks violate CLAUDE.md and waste disk.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'input.vcf' + v = [(i * 10, f'PCGR_CSQ={_csq("intron_variant")}') for i in range(1, 26)] + _write_vcf(vcf_fp, v) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 10): + chunks = pcgr.split_vcf(vcf_fp, tmp_path) + + for chunk in chunks: + self.assertTrue(str(chunk).endswith('.vcf.gz'), + f'Expected .vcf.gz chunk, got: {chunk.name}') + + def test_chunks_are_tabix_indexed(self): + """Each .vcf.gz chunk must have a .tbi index — PCGR v2.2.5 requires it.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'input.vcf' + v = [(i * 10, f'PCGR_CSQ={_csq("intron_variant")}') for i in range(1, 26)] + _write_vcf(vcf_fp, v) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 10): + chunks = pcgr.split_vcf(vcf_fp, tmp_path) + + for chunk in chunks: + tbi = pathlib.Path(str(chunk) + '.tbi') + self.assertTrue(tbi.exists(), f'Missing tabix index for {chunk.name}') + + +class TestRunSomaticChunkArgMapping(unittest.TestCase): + """Regression test: run_somatic_chunk must forward args as keywords to run_somatic. + + Before the fix, run_somatic_chunk called run_somatic positionally (6 args), + skipping pcgr_threads. This caused pcgr_conda ('pcgr') to land in the + pcgr_threads slot → ValueError: invalid literal for int() with base 10: 'pcgr'. + """ + + def test_pcgr_conda_not_shifted_into_pcgr_threads(self): + """pcgr_conda must reach run_somatic as pcgr_conda, not as pcgr_threads.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'chunk.vcf' + _write_vcf(vcf_fp, [(10, f'PCGR_CSQ={_csq("intron_variant")}')]) + + captured = {} + + def fake_run_somatic(*args, **kwargs): + captured['args'] = args + captured['kwargs'] = kwargs + return (None, None) + + with patch('bolt.common.pcgr.run_somatic', side_effect=fake_run_somatic), \ + patch('bolt.common.pcgr.merging_pcgr_files', + return_value=(tmp_path / 'out.vcf', tmp_path / 'out.tsv')): + pcgr.run_somatic_chunk( + [vcf_fp], + pcgr_data_dir=tmp_path / 'pcgr_data', + vep_dir=tmp_path / 'vep', + output_dir=tmp_path, + pcgr_output_dir=tmp_path / 'pcgr_output', + max_threads=4, + pcgr_conda='pcgr_env', + pcgrr_conda='pcgrr_env', + ) + + kw = captured['kwargs'] + self.assertEqual(kw.get('pcgr_conda'), 'pcgr_env', + 'pcgr_conda was not forwarded — likely shifted into pcgr_threads') + self.assertEqual(kw.get('pcgrr_conda'), 'pcgrr_env', + 'pcgrr_conda was not forwarded correctly') + self.assertEqual(kw.get('threads'), 4, + 'threads (max_threads) was not forwarded correctly') + self.assertEqual(kw.get('chunk_nbr'), 1, + 'chunk_nbr was not forwarded correctly') + + def test_disable_estimates_passed_to_run_somatic(self): + """run_somatic_chunk must pass disable_estimates=True to every run_somatic call. + + Chunked PCGR runs must not include --estimate_msi/--estimate_tmb per chunk — + those flags produce per-chunk partial estimates that are meaningless after merging. + Fixed in bolt 0.3.2 (umccr/sash#57 + disable_estimates wiring). + """ + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'chunk.vcf' + _write_vcf(vcf_fp, [(10, f'PCGR_CSQ={_csq("intron_variant")}')]) + + captured = {} + + def fake_run_somatic(*args, **kwargs): + captured['kwargs'] = kwargs + return (None, None) + + with patch('bolt.common.pcgr.run_somatic', side_effect=fake_run_somatic), \ + patch('bolt.common.pcgr.merging_pcgr_files', + return_value=(tmp_path / 'out.tsv', tmp_path / 'out.vcf')): + pcgr.run_somatic_chunk( + [vcf_fp], + pcgr_data_dir=tmp_path / 'pcgr_data', + vep_dir=tmp_path / 'vep', + output_dir=tmp_path, + pcgr_output_dir=tmp_path / 'pcgr_output', + max_threads=4, + pcgr_conda='pcgr_env', + pcgrr_conda='pcgrr_env', + ) + + self.assertTrue( + captured['kwargs'].get('disable_estimates'), + 'run_somatic_chunk must forward disable_estimates=True — ' + 'per-chunk MSI/TMB estimates are meaningless after merge', + ) + + +class TestMergingPcgrFiles(unittest.TestCase): + """Regression test for bolt #26: bcftools merge requires 2+ inputs. + + When a sample's variants fit in a single PCGR chunk, run_somatic_chunk still + called merging_pcgr_files() -> util.merge_vcf_files() unconditionally, which + invoked `bcftools merge` on a single VCF and errored (Usage: bcftools merge + [options] [...]). + """ + + def _write_gz_vcf(self, path, variants): + vcf_path = path.with_suffix('') + _write_vcf(vcf_path, variants) + util.execute_command(f'bcftools view -Oz -o {path} {vcf_path}') + util.execute_command(f'bcftools index -t {path}') + + def test_single_chunk_skips_bcftools_merge(self): + """A single VCF chunk must bypass bcftools merge and pass through directly.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'chunk1.vcf.gz' + self._write_gz_vcf(vcf_fp, [(10, f'PCGR_CSQ={_csq("intron_variant")}')]) + + tsv_fp = tmp_path / 'chunk1.tsv.gz' + with gzip.open(tsv_fp, 'wt') as fh: + fh.write('col1\tcol2\nval1\tval2\n') + + merged_vcf, merged_tsv = pcgr.merging_pcgr_files(tmp_path, [vcf_fp], [tsv_fp]) + + self.assertTrue(pathlib.Path(merged_vcf).exists()) + self.assertEqual(_count_vcf(merged_vcf), 1) + self.assertTrue(pathlib.Path(f'{merged_vcf}.tbi').exists(), + 'Single-chunk pass-through VCF must still be tabix indexed') + self.assertTrue(pathlib.Path(merged_tsv).exists()) + + def test_multiple_chunks_still_merge(self): + """Two or more chunks must still go through bcftools merge as before.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf1_fp = tmp_path / 'chunk1.vcf.gz' + vcf2_fp = tmp_path / 'chunk2.vcf.gz' + self._write_gz_vcf(vcf1_fp, [(10, f'PCGR_CSQ={_csq("intron_variant")}')]) + self._write_gz_vcf(vcf2_fp, [(20, f'PCGR_CSQ={_csq("intron_variant")}')]) + + tsv1_fp = tmp_path / 'chunk1.tsv.gz' + tsv2_fp = tmp_path / 'chunk2.tsv.gz' + with gzip.open(tsv1_fp, 'wt') as fh: + fh.write('col1\tcol2\nval1\tval2\n') + with gzip.open(tsv2_fp, 'wt') as fh: + fh.write('col1\tcol2\nval3\tval4\n') + + merged_vcf, merged_tsv = pcgr.merging_pcgr_files( + tmp_path, [vcf1_fp, vcf2_fp], [tsv1_fp, tsv2_fp] + ) + + self.assertTrue(pathlib.Path(merged_vcf).exists()) + self.assertTrue(pathlib.Path(merged_tsv).exists()) + + +class TestRunSomaticCommandArgs(unittest.TestCase): + """Verify run_somatic builds the correct PCGR command-line arguments.""" + + def test_estimate_signatures_absent_from_command(self): + """--estimate_signatures must not appear in the PCGR command (dropped in sash#57). + + --estimate_msi and --estimate_tmb must still be present for non-chunked runs. + """ + captured = {} + + def fake_execute(cmd, **kwargs): + captured['cmd'] = cmd + # create the output files run_somatic expects to find after pcgr runs + output_dir.mkdir(parents=True, exist_ok=True) + (output_dir / 'nosampleset.pcgr.grch38.snv_indel_ann.tsv.gz').touch() + (output_dir / 'nosampleset.pcgr.grch38.pass.vcf.gz').touch() + + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + output_dir = tmp_path / 'output' + + with patch('bolt.common.pcgr.util.execute_command', side_effect=fake_execute): + pcgr.run_somatic( + input_fp=tmp_path / 'input.vcf.gz', + pcgr_refdata_dir=tmp_path / 'refdata', + vep_dir=tmp_path / 'vep', + output_dir=output_dir, + ) + + self.assertIn('cmd', captured, 'execute_command was not called') + self.assertNotIn('--estimate_signatures', captured['cmd'], + '--estimate_signatures must be absent (dropped in sash#57)') + self.assertIn('--estimate_msi', captured['cmd'], + '--estimate_msi must still be present for non-chunked runs') + self.assertIn('--estimate_tmb', captured['cmd'], + '--estimate_tmb must still be present for non-chunked runs') + + def test_disable_estimates_suppresses_msi_tmb(self): + """disable_estimates=True must suppress --estimate_msi and --estimate_tmb.""" + captured = {} + + def fake_execute(cmd, **kwargs): + captured['cmd'] = cmd + output_dir.mkdir(parents=True, exist_ok=True) + (output_dir / 'nosampleset.pcgr.grch38.snv_indel_ann.tsv.gz').touch() + (output_dir / 'nosampleset.pcgr.grch38.pass.vcf.gz').touch() + + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + output_dir = tmp_path / 'output' + + with patch('bolt.common.pcgr.util.execute_command', side_effect=fake_execute): + pcgr.run_somatic( + input_fp=tmp_path / 'input.vcf.gz', + pcgr_refdata_dir=tmp_path / 'refdata', + vep_dir=tmp_path / 'vep', + output_dir=output_dir, + disable_estimates=True, + ) + + self.assertNotIn('--estimate_msi', captured['cmd'], + '--estimate_msi must be absent when disable_estimates=True') + self.assertNotIn('--estimate_tmb', captured['cmd'], + '--estimate_tmb must be absent when disable_estimates=True') + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_pcgr_annotation.py b/tests/test_pcgr_annotation.py new file mode 100644 index 0000000..17b2ec8 --- /dev/null +++ b/tests/test_pcgr_annotation.py @@ -0,0 +1,428 @@ +"""Unit tests for bolt/common/pcgr.py pure annotation-transfer helpers. + +Covers only binary-free / cyvcf2-in-memory logic. Functions requiring +bcftools/PCGR/vcfanno subprocesses are intentionally NOT tested here. +Does NOT duplicate tests already in test_pcgr.py / test_smlv_somatic_report.py +(get_ordering, get_impacts, determine_filter, get_variant_filter_data, +split_vcf, run_somatic_chunk, count_variant_process, select_pcgr_variants). +""" +import pathlib +import tempfile +import unittest + +import cyvcf2 + +import bolt.common.constants as constants +import bolt.common.pcgr as pcgr +import bolt.util as util + + +# Minimal VCF header with INFO tags exercised by annotate_record tests +HEADER = ( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + f'##INFO=\n' + f'##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' +) + + +def _write_vcf(path, variants): + with open(path, 'w') as fh: + fh.write(HEADER) + for chrom, pos, ref, alt, info in variants: + fh.write(f'{chrom}\t{pos}\t.\t{ref}\t{alt}\t.\tPASS\t{info}\n') + + +def _make_variant_with_tags(tags, chrom='chr1', pos=100, ref='A', alt='T', info_str='.'): + """Return a cyvcf2 Variant with a custom set of INFO enum tags registered.""" + with tempfile.TemporaryDirectory() as tmp: + vcf_path = pathlib.Path(tmp) / 'test.vcf' + _write_vcf(vcf_path, [(chrom, pos, ref, alt, info_str)]) + fh = cyvcf2.VCF(str(vcf_path)) + for tag_enum in tags: + util.add_vcf_header_entry(fh, tag_enum) + return list(fh)[0] + + +def _make_variant(chrom='chr1', pos=100, ref='A', alt='T', info_str='.'): + """Return a cyvcf2 Variant with default INFO tags registered.""" + return _make_variant_with_tags( + [constants.VcfInfo.PCGR_ACTIONABILITY_TIER, constants.VcfInfo.PCGR_CSQ], + chrom=chrom, pos=pos, ref=ref, alt=alt, info_str=info_str, + ) + + +class TestParseGenomicChange(unittest.TestCase): + """Unit tests for pcgr.parse_genomic_change().""" + + def test_standard_substitution(self): + chrom, pos, ref, alt = pcgr.parse_genomic_change('3:g.41224645T>C') + self.assertEqual(chrom, 'chr3') + self.assertEqual(pos, 41224645) + self.assertEqual(ref, 'T') + self.assertEqual(alt, 'C') + + def test_chr_prefix_added(self): + chrom, _, _, _ = pcgr.parse_genomic_change('17:g.7674220C>T') + self.assertTrue(chrom.startswith('chr'), f"Expected 'chr' prefix, got: {chrom}") + + def test_multichar_alleles(self): + chrom, pos, ref, alt = pcgr.parse_genomic_change('1:g.100ACGT>TTTT') + self.assertEqual(chrom, 'chr1') + self.assertEqual(ref, 'ACGT') + self.assertEqual(alt, 'TTTT') + + def test_returns_int_pos(self): + _, pos, _, _ = pcgr.parse_genomic_change('3:g.41224645T>C') + self.assertIsInstance(pos, int) + + def test_raises_on_garbage_input(self): + with self.assertRaises(ValueError): + pcgr.parse_genomic_change('garbage') + + def test_raises_on_partial_format(self): + with self.assertRaises(ValueError): + pcgr.parse_genomic_change('3:g.41224645') + + +class TestGetImpactsHigher(unittest.TestCase): + """Unit tests for pcgr.get_impacts_higher(). + + Expectations are derived from constants.VEP_IMPACTS_FILTER rather than + hard-coded literals, so the test tracks the source of truth. + """ + + def test_returns_higher_impacts_for_first_element(self): + first = constants.VEP_IMPACTS_FILTER[0] + higher = pcgr.get_impacts_higher(first) + expected = list(constants.VEP_IMPACTS_FILTER[1:]) + self.assertEqual(list(higher), expected) + + def test_returns_empty_for_last_element(self): + last = constants.VEP_IMPACTS_FILTER[-1] + higher = pcgr.get_impacts_higher(last) + self.assertEqual(list(higher), []) + + def test_middle_element(self): + if len(constants.VEP_IMPACTS_FILTER) < 3: + self.skipTest('VEP_IMPACTS_FILTER has fewer than 3 elements') + mid_idx = len(constants.VEP_IMPACTS_FILTER) // 2 + mid = constants.VEP_IMPACTS_FILTER[mid_idx] + higher = pcgr.get_impacts_higher(mid) + expected = list(constants.VEP_IMPACTS_FILTER[mid_idx + 1:]) + self.assertEqual(list(higher), expected) + + def test_does_not_include_self(self): + impact = constants.VEP_IMPACTS_FILTER[0] + higher = pcgr.get_impacts_higher(impact) + self.assertNotIn(impact, higher) + + +class TestGetAnnotationEntryTsv(unittest.TestCase): + """Unit tests for pcgr.get_annotation_entry_tsv().""" + + # Map used in tests: one field from INFO enum present in VCF_HEADER_ENTRIES + INFO_FIELD_MAP = { + constants.VcfInfo.PCGR_CSQ: 'CSQ', + } + + def test_genomic_change_key_extraction(self): + record = { + 'GENOMIC_CHANGE': '3:g.41224645T>C', + 'CSQ': 'some_consequence', + } + key, record_ann = pcgr.get_annotation_entry_tsv(record, self.INFO_FIELD_MAP) + self.assertEqual(key, ('chr3', 41224645, 'T', 'C')) + + def test_chrom_pos_ref_alt_fallback_with_chr_prefix(self): + record = { + 'CHROM': '1', + 'POS': '100', + 'REF': 'A', + 'ALT': 'T', + 'CSQ': 'some_consequence', + } + key, record_ann = pcgr.get_annotation_entry_tsv(record, self.INFO_FIELD_MAP) + chrom, pos, ref, alt = key + self.assertTrue(str(chrom).startswith('chr'), f"Expected 'chr' prefix, got: {chrom}") + + def test_na_values_omitted_from_record_ann(self): + record = { + 'GENOMIC_CHANGE': '3:g.41224645T>C', + 'CSQ': 'NA', + } + _, record_ann = pcgr.get_annotation_entry_tsv(record, self.INFO_FIELD_MAP) + self.assertNotIn(constants.VcfInfo.PCGR_CSQ, record_ann) + + def test_falsy_values_omitted_from_record_ann(self): + record = { + 'GENOMIC_CHANGE': '3:g.41224645T>C', + 'CSQ': '', + } + _, record_ann = pcgr.get_annotation_entry_tsv(record, self.INFO_FIELD_MAP) + self.assertNotIn(constants.VcfInfo.PCGR_CSQ, record_ann) + + def test_valid_value_present_in_record_ann(self): + record = { + 'GENOMIC_CHANGE': '3:g.41224645T>C', + 'CSQ': 'A|intron_variant|...', + } + _, record_ann = pcgr.get_annotation_entry_tsv(record, self.INFO_FIELD_MAP) + self.assertIn(constants.VcfInfo.PCGR_CSQ, record_ann) + self.assertEqual(record_ann[constants.VcfInfo.PCGR_CSQ], 'A|intron_variant|...') + + +class TestCompileAnnotationData(unittest.TestCase): + """Unit tests for pcgr.compile_annotation_data().""" + + def test_tsv_wins_on_key_field_collision(self): + key = ('chr1', 100, 'A', 'T') + data_tsv = {key: {constants.VcfInfo.PCGR_CSQ: 'tsv_value'}} + data_vcf = {key: {constants.VcfInfo.PCGR_CSQ: 'vcf_value'}} + result = pcgr.compile_annotation_data(data_tsv, data_vcf) + self.assertEqual(result[key][constants.VcfInfo.PCGR_CSQ], 'tsv_value') + + def test_vcf_only_key_added(self): + key_tsv = ('chr1', 100, 'A', 'T') + key_vcf = ('chr2', 200, 'G', 'C') + data_tsv = {key_tsv: {constants.VcfInfo.PCGR_CSQ: 'tsv_value'}} + data_vcf = {key_vcf: {constants.VcfInfo.PCGR_CSQ: 'vcf_value'}} + result = pcgr.compile_annotation_data(data_tsv, data_vcf) + self.assertIn(key_vcf, result) + self.assertEqual(result[key_vcf][constants.VcfInfo.PCGR_CSQ], 'vcf_value') + + def test_vcf_only_field_added_to_existing_key(self): + key = ('chr1', 100, 'A', 'T') + data_tsv = {key: {constants.VcfInfo.PCGR_ACTIONABILITY_TIER: '1'}} + data_vcf = {key: {constants.VcfInfo.PCGR_CSQ: 'vcf_csq'}} + result = pcgr.compile_annotation_data(data_tsv, data_vcf) + self.assertEqual(result[key][constants.VcfInfo.PCGR_ACTIONABILITY_TIER], '1') + self.assertEqual(result[key][constants.VcfInfo.PCGR_CSQ], 'vcf_csq') + + def test_empty_inputs(self): + result = pcgr.compile_annotation_data({}, {}) + self.assertEqual(result, {}) + + def test_tsv_only_key_preserved(self): + key = ('chr1', 100, 'A', 'T') + data_tsv = {key: {constants.VcfInfo.PCGR_CSQ: 'tsv_value'}} + result = pcgr.compile_annotation_data(data_tsv, {}) + self.assertIn(key, result) + self.assertEqual(result[key][constants.VcfInfo.PCGR_CSQ], 'tsv_value') + + +class TestAnnotateRecord(unittest.TestCase): + """Unit tests for pcgr.annotate_record().""" + + def test_annotations_written_to_matching_record(self): + variant = _make_variant_with_tags( + [constants.VcfInfo.PCGR_ACTIONABILITY_TIER], + chrom='chr1', pos=100, ref='A', alt='T', + ) + annotations = { + ('chr1', 100, 'A', 'T'): { + constants.VcfInfo.PCGR_ACTIONABILITY_TIER: '2', + } + } + result = pcgr.annotate_record(variant, annotations, allow_missing=True) + self.assertEqual( + result.INFO.get(constants.VcfInfo.PCGR_ACTIONABILITY_TIER.value), + '2', + ) + + def test_missing_key_allow_missing_true_returns_record_unchanged(self): + variant = _make_variant(chrom='chr1', pos=100, ref='A', alt='T') + annotations = {} + result = pcgr.annotate_record(variant, annotations, allow_missing=True) + self.assertIsNotNone(result) + + def test_missing_key_allow_missing_false_raises_assertion(self): + variant = _make_variant(chrom='chr1', pos=100, ref='A', alt='T') + annotations = {} + with self.assertRaises(AssertionError): + pcgr.annotate_record(variant, annotations, allow_missing=False) + + def test_key_uses_chrom_pos_ref_alt(self): + """annotate_record must match on the exact (CHROM, POS, REF, ALT) tuple.""" + variant = _make_variant_with_tags( + [constants.VcfInfo.PCGR_ACTIONABILITY_TIER], + chrom='chr5', pos=999, ref='G', alt='C', + ) + annotations = { + ('chr5', 999, 'G', 'C'): { + constants.VcfInfo.PCGR_ACTIONABILITY_TIER: '3', + } + } + result = pcgr.annotate_record(variant, annotations, allow_missing=False) + self.assertEqual( + result.INFO.get(constants.VcfInfo.PCGR_ACTIONABILITY_TIER.value), + '3', + ) + + +class TestGetAnnotationsVcf(unittest.TestCase): + """Unit tests for pcgr.get_annotations_vcf() duplicate-key handling. + + PCGR strips the 'chr' prefix from its own output VCF, so fixtures here use + bare contig names ('1', not 'chr1') to match real PCGR output. + """ + + PCGR_VCF_HEADER = ( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + f'##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + ) + + def _write_pcgr_vcf(self, path, rows): + with open(path, 'w') as fh: + fh.write(self.PCGR_VCF_HEADER) + for chrom, pos, ref, alt, info in rows: + fh.write(f'{chrom}\t{pos}\t.\t{ref}\t{alt}\t.\tPASS\t{info}\n') + + def test_duplicate_vcf_key_keeps_first_no_crash(self): + """A duplicate variant in PCGR's VCF output must be skipped, not raise AssertionError.""" + info_field_map = {constants.VcfInfo.PCGR_CSQ: 'PCGR_CSQ'} + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'pcgr.vcf' + self._write_pcgr_vcf(vcf_fp, [ + ('1', 100, 'A', 'T', 'PCGR_CSQ=first'), + ('1', 100, 'A', 'T', 'PCGR_CSQ=second'), # duplicate transcript mapping + ]) + result = pcgr.get_annotations_vcf(vcf_fp, info_field_map) + + key = ('chr1', 100, 'A', 'T') + self.assertEqual(len(result), 1) + self.assertEqual(result[key][constants.VcfInfo.PCGR_CSQ], 'first') + + def test_non_duplicate_vcf_keys_all_present(self): + """Distinct variants are all retained unaffected by dedup handling.""" + info_field_map = {constants.VcfInfo.PCGR_CSQ: 'PCGR_CSQ'} + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'pcgr.vcf' + self._write_pcgr_vcf(vcf_fp, [ + ('1', 100, 'A', 'T', 'PCGR_CSQ=v1'), + ('1', 200, 'C', 'G', 'PCGR_CSQ=v2'), + ]) + result = pcgr.get_annotations_vcf(vcf_fp, info_field_map) + + self.assertEqual(len(result), 2) + self.assertEqual(result[('chr1', 100, 'A', 'T')][constants.VcfInfo.PCGR_CSQ], 'v1') + self.assertEqual(result[('chr1', 200, 'C', 'G')][constants.VcfInfo.PCGR_CSQ], 'v2') + + +def _write_empty_pcgr_vcf(path): + """Header-only sites-only VCF, so get_annotations_vcf() returns {}.""" + with open(path, 'w') as fh: + fh.write( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + f'##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + ) + + +def _write_tsv(path, header, rows): + import csv as _csv + with open(path, 'w', newline='') as fh: + writer = _csv.DictWriter(fh, fieldnames=header, delimiter='\t') + writer.writeheader() + for row in rows: + writer.writerow(row) + + +class TestCollectPcgrAnnotationData(unittest.TestCase): + """Unit tests for pcgr.collect_pcgr_annotation_data() duplicate-key tier resolution.""" + + HEADER = ['GENOMIC_CHANGE', 'CSQ', 'ACTIONABILITY_TIER'] + + def _collect(self, tmp, rows): + tsv_fp = pathlib.Path(tmp) / 'pcgr.tsv' + vcf_fp = pathlib.Path(tmp) / 'pcgr.vcf' + _write_tsv(tsv_fp, self.HEADER, rows) + _write_empty_pcgr_vcf(vcf_fp) + info_field_map = {constants.VcfInfo.PCGR_CSQ: 'CSQ'} + return pcgr.collect_pcgr_annotation_data(tsv_fp, vcf_fp, info_field_map) + + def test_duplicate_key_keeps_more_actionable_tier_no_crash(self): + """A second, more-actionable-tier row for the same variant replaces the first.""" + with tempfile.TemporaryDirectory() as tmp: + result = self._collect(tmp, [ + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'transcript_a', 'ACTIONABILITY_TIER': 'TIER 3'}, + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'transcript_b', 'ACTIONABILITY_TIER': 'TIER 1'}, + ]) + key = ('chr1', 100, 'A', 'T') + self.assertEqual(len(result), 1) + self.assertEqual(result[key][constants.VcfInfo.PCGR_CSQ], 'transcript_b') + self.assertEqual(result[key][constants.VcfInfo.PCGR_ACTIONABILITY_TIER], '1') + + def test_duplicate_key_skips_less_actionable_tier(self): + """A second, less-actionable-tier row for the same variant is dropped.""" + with tempfile.TemporaryDirectory() as tmp: + result = self._collect(tmp, [ + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'transcript_a', 'ACTIONABILITY_TIER': 'TIER 1'}, + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'transcript_b', 'ACTIONABILITY_TIER': 'TIER 3'}, + ]) + key = ('chr1', 100, 'A', 'T') + self.assertEqual(len(result), 1) + self.assertEqual(result[key][constants.VcfInfo.PCGR_CSQ], 'transcript_a') + self.assertEqual(result[key][constants.VcfInfo.PCGR_ACTIONABILITY_TIER], '1') + + def test_non_duplicate_keys_all_present(self): + with tempfile.TemporaryDirectory() as tmp: + result = self._collect(tmp, [ + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'v1', 'ACTIONABILITY_TIER': 'TIER 1'}, + {'GENOMIC_CHANGE': '1:g.200C>G', 'CSQ': 'v2', 'ACTIONABILITY_TIER': 'TIER 2'}, + ]) + self.assertEqual(len(result), 2) + self.assertIn(('chr1', 100, 'A', 'T'), result) + self.assertIn(('chr1', 200, 'C', 'G'), result) + + +class TestCollectCpsrAnnotationData(unittest.TestCase): + """Unit tests for pcgr.collect_cpsr_annotation_data() duplicate-key handling.""" + + HEADER = ['GENOMIC_CHANGE', 'CSQ'] + + def _collect(self, tmp, rows): + tsv_fp = pathlib.Path(tmp) / 'cpsr.tsv.gz' + vcf_fp = pathlib.Path(tmp) / 'cpsr.vcf' + import gzip as _gzip + import csv as _csv + with _gzip.open(tsv_fp, 'wt', newline='') as fh: + writer = _csv.DictWriter(fh, fieldnames=self.HEADER, delimiter='\t') + writer.writeheader() + for row in rows: + writer.writerow(row) + _write_empty_pcgr_vcf(vcf_fp) + info_field_map = {constants.VcfInfo.PCGR_CSQ: 'CSQ'} + return pcgr.collect_cpsr_annotation_data(tsv_fp, vcf_fp, info_field_map) + + def test_duplicate_key_keeps_first_no_crash(self): + """A duplicate CPSR TSV row for the same variant must be skipped, not raise.""" + with tempfile.TemporaryDirectory() as tmp: + result = self._collect(tmp, [ + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'first'}, + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'second'}, + ]) + key = ('chr1', 100, 'A', 'T') + self.assertEqual(len(result), 1) + self.assertEqual(result[key][constants.VcfInfo.PCGR_CSQ], 'first') + + def test_non_duplicate_keys_all_present(self): + with tempfile.TemporaryDirectory() as tmp: + result = self._collect(tmp, [ + {'GENOMIC_CHANGE': '1:g.100A>T', 'CSQ': 'v1'}, + {'GENOMIC_CHANGE': '1:g.200C>G', 'CSQ': 'v2'}, + ]) + self.assertEqual(len(result), 2) + self.assertIn(('chr1', 100, 'A', 'T'), result) + self.assertIn(('chr1', 200, 'C', 'G'), result) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_smlv_somatic_report.py b/tests/test_smlv_somatic_report.py new file mode 100644 index 0000000..78dd9c2 --- /dev/null +++ b/tests/test_smlv_somatic_report.py @@ -0,0 +1,417 @@ +"""Tests for bolt/workflows/smlv_somatic/report.py — hypermutated variant trimming.""" +import pathlib +import shutil +import tempfile +import unittest +from unittest.mock import patch + +from click.testing import CliRunner + +import cyvcf2 + +import bolt.common.constants as constants +import bolt.common.pcgr as pcgr +import bolt.util as util +import bolt.workflows.smlv_somatic.report as report_mod + +from tests.helpers import _csq, _write_vcf + + +class TestSelectPcgrVariants(unittest.TestCase): + """Integration tests for select_pcgr_variants() trimming logic.""" + + def _run(self, variants, limit, tmp): + """Run select_pcgr_variants with a small MAX_SOMATIC_VARIANTS limit.""" + vcf_fp = pathlib.Path(tmp) / 'input.vcf' + _write_vcf(vcf_fp, variants) + cancer_genes = pathlib.Path(tmp) / 'genes.bed' + cancer_genes.write_text('chr1\t1\t9999999\n') + + # Mock bcftools annotate: copy input to the expected output path + orig_execute = util.execute_command + def fake_execute(cmd, **_): + import re + m = re.search(r'--output\s+(\S+)', cmd) + if m and 'bcftools annotate' in cmd: + shutil.copy(str(vcf_fp), m.group(1)) + else: + orig_execute(cmd) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', limit), \ + patch('bolt.util.execute_command', side_effect=fake_execute): + out_fp = report_mod.select_pcgr_variants( + vcf_fp, cancer_genes, 'TUMOR', pathlib.Path(tmp) + ) + return sum(1 for _ in cyvcf2.VCF(str(out_fp))) + + def test_output_within_limit(self): + """Output must never exceed MAX_SOMATIC_VARIANTS.""" + with tempfile.TemporaryDirectory() as tmp: + # 15 variants: 2 hotspot + 5 TIER_1 + 5 TIER_3 + 3 NONCODING + v = [] + for i in range(1, 3): # hotspot + v.append((i*10, f'HMF_HOTSPOT;PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(3, 8): # TIER_1 intronic + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(8, 13): # TIER_3 intronic + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=3;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(13, 16): # NONCODING intergenic + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + count = self._run(v, limit=10, tmp=tmp) + self.assertLessEqual(count, 10) + + def test_noncoding_dropped_before_tier1(self): + """With limit = total - 3, the 3 NONCODING variants should be dropped (not TIER_1).""" + with tempfile.TemporaryDirectory() as tmp: + v = [] + for i in range(1, 6): # 5 TIER_1 intronic + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(6, 9): # 3 NONCODING intergenic + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + # limit=5: should drop the 3 NONCODING to get to 5 + count = self._run(v, limit=5, tmp=tmp) + self.assertEqual(count, 5) + + def test_hotspots_never_dropped_by_tiered_filter(self): + """Hotspot variants must survive tiered filtering. + + Note: HMF_HOTSPOT is not in RETAIN_FIELDS_FILTERING — these variants survive because + they are TIER_1 (highest priority), not via the hotspot retention path. + """ + with tempfile.TemporaryDirectory() as tmp: + v = [] + for i in range(1, 3): # 2 HMF_HOTSPOT TIER_1 variants (survive via tier priority) + v.append((i*10, f'HMF_HOTSPOT;PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(3, 13): # 10 NONCODING (should all be filtered) + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + # limit=2: only the 2 TIER_1 variants should remain + count = self._run(v, limit=2, tmp=tmp) + self.assertEqual(count, 2) + + def test_all_within_limit_nothing_filtered(self): + """When total variants are below the limit, nothing is dropped.""" + with tempfile.TemporaryDirectory() as tmp: + v = [(i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}') + for i in range(1, 6)] # 5 NONCODING + count = self._run(v, limit=10, tmp=tmp) + self.assertEqual(count, 5) + + def test_retained_variants_bypass_tiered_filter(self): + """Variants with PANEL or SAGE_HOTSPOT bypass tiered filtering and always survive.""" + with tempfile.TemporaryDirectory() as tmp: + v = [] + for i in range(1, 3): # 2 SAGE_HOTSPOT + v.append((i*10, 'SAGE_HOTSPOT')) + for i in range(3, 5): # 2 PANEL-only + v.append((i*10, 'PANEL')) + for i in range(5, 15): # 10 NONCODING that get dropped + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + count = self._run(v, limit=4, tmp=tmp) + self.assertEqual(count, 4) + + def test_pcgr_mutation_hotspot_real_value_is_retained(self): + """A real PCGR_MUTATION_HOTSPOT value (non-dot) must retain the variant.""" + with tempfile.TemporaryDirectory() as tmp: + v = [] + for i in range(1, 3): # 2 real hotspot variants — must survive + v.append((i*10, f'PCGR_MUTATION_HOTSPOT=GRCH38_1_{i}_A_T;PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(3, 8): # 5 NONCODING — dropped + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + count = self._run(v, limit=2, tmp=tmp) + self.assertEqual(count, 2) + + def test_pcgr_mutation_hotspot_dot_not_treated_as_retained(self): + """PCGR_MUTATION_HOTSPOT=. must not retain variants — '.' is a missing-value placeholder. + + cyvcf2 returns the string '.' (truthy) for String INFO fields written as '=.' by PCGR on + every non-hotspot variant. Without the fix, any(variant.INFO.get(e) ...) always returns + True and ALL variants are treated as retained, so tiered filtering never drops anything and + RuntimeError fires for any sample with >450k variants (sash #52 root cause). + """ + with tempfile.TemporaryDirectory() as tmp: + v = [] + for i in range(1, 3): # 2 real SAGE_HOTSPOT — must be retained + v.append((i*10, f'SAGE_HOTSPOT;PCGR_MUTATION_HOTSPOT=.;PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(3, 8): # 5 NONCODING with PCGR_MUTATION_HOTSPOT=. — must be droppable + v.append((i*10, f'PCGR_MUTATION_HOTSPOT=.;PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + # limit=2: only the 2 real SAGE_HOTSPOT variants survive; the 5 dot-placeholder ones are dropped + count = self._run(v, limit=2, tmp=tmp) + self.assertEqual(count, 2) + + def test_filters_set_vcf_marks_dropped_variants(self): + """The traceability VCF marks filtered-out variants with PCGR_count_limit. + + The function drops entire categories, so we need two distinct categories: + - 3 TIER_1 intronic (high priority — kept) + - 2 NONCODING intergenic (lowest priority — dropped as a whole category) + """ + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'input.vcf' + v = [] + for i in range(1, 4): # 3 TIER_1 intronic — kept + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(4, 6): # 2 NONCODING intergenic — dropped whole category + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + _write_vcf(vcf_fp, v) + cancer_genes = tmp_path / 'genes.bed' + cancer_genes.write_text('chr1\t1\t9999999\n') + + orig_execute = util.execute_command + def fake_execute(cmd, **_): + import re + m = re.search(r'--output\s+(\S+)', cmd) + if m and 'bcftools annotate' in cmd: + shutil.copy(str(vcf_fp), m.group(1)) + else: + orig_execute(cmd) + + # limit=3: the 2 NONCODING category is dropped, 3 TIER_1 survive + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 3), \ + patch('bolt.util.execute_command', side_effect=fake_execute): + report_mod.select_pcgr_variants(vcf_fp, cancer_genes, 'TUMOR', tmp_path) + + filters_set_fp = tmp_path / 'TUMOR.pcgr_hypermutated.filters_set.vcf.gz' + self.assertTrue(filters_set_fp.exists(), 'filters_set VCF not created') + + all_records = list(cyvcf2.VCF(str(filters_set_fp))) + self.assertEqual(len(all_records), 5, 'filters_set VCF should contain all input variants') + + filter_tag = constants.VcfFilter.PCGR_COUNT_LIMIT.value + dropped = [r for r in all_records if filter_tag in (r.FILTERS or [])] + self.assertEqual(len(dropped), 2, 'Expected 2 NONCODING variants marked with PCGR_count_limit') + + +class TestSelectPcgrVariantsRaisesOnUnresolvableOverflow(unittest.TestCase): + """select_pcgr_variants raises RuntimeError when all variants are retained (hotspots).""" + + def _fake_execute(self, vcf_fp): + def _run(cmd, **_): + import re + m = re.search(r'--output\s+(\S+)', cmd) + if m and 'bcftools annotate' in cmd: + shutil.copy(str(vcf_fp), m.group(1)) + else: + util.execute_command(cmd) + return _run + + def test_raises_when_all_variants_are_hotspots(self): + """All SAGE_HOTSPOT variants are RETAIN_FIELDS — tiered filtering cannot drop any; RuntimeError expected.""" + # RETAIN_FIELDS_FILTERING includes SAGE_HOTSPOT — use that flag, not HMF_HOTSPOT + HOTSPOT_INFO = f'SAGE_HOTSPOT;PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}' + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + vcf_fp = tmp_path / 'input.vcf' + variants = [(i * 10, HOTSPOT_INFO) for i in range(1, 6)] + _write_vcf(vcf_fp, variants) + cancer_genes = tmp_path / 'genes.bed' + cancer_genes.write_text('chr1\t1\t9999999\n') + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 3), \ + patch('bolt.util.execute_command', side_effect=self._fake_execute(vcf_fp)): + with self.assertRaises(RuntimeError): + report_mod.select_pcgr_variants(vcf_fp, cancer_genes, 'TUMOR', tmp_path) + + +class TestCountVariantProcess(unittest.TestCase): + """Verify count_variant_process counts and is_hypermutated flag (bolt #27). + + is_hypermutated must use the 'dragen' count (raw, pre-bolt-filter), not + 'filter_pass'. A sample with many DRAGEN variants that are mostly filtered + away must still be flagged as hypermutated. + """ + + # Minimal header for count_variant_process: needs FILTER tags + SAGE_NOVEL INFO + COUNT_HEADER = ( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + f'##FILTER=\n' + f'##FILTER=\n' + f'##FILTER=\n' + f'##INFO=\n' + f'##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + ) + + def _write_count_vcf(self, path, rows): + """rows: list of (pos, filter_str, info_str) tuples.""" + with open(path, 'w') as fh: + fh.write(self.COUNT_HEADER) + for pos, filt, info in rows: + fh.write(f'chr1\t{pos}\t.\tA\tT\t.\t{filt}\t{info}\n') + + def test_is_hypermutated_uses_dragen_count(self): + """is_hypermutated=True when dragen count > MAX_SOMATIC_VARIANTS even if filter_pass is below.""" + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'test.vcf' + min_af = constants.VcfFilter.MIN_AF.value + # 3 DRAGEN PASS variants + 2 filtered by bolt (MIN_AF) — filter_pass=3, dragen=5 + rows = [(i * 10, 'PASS', '.') for i in range(1, 4)] + rows += [(i * 10 + 5, min_af, '.') for i in range(1, 3)] + self._write_count_vcf(vcf_fp, rows) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 4): + counts = report_mod.count_variant_process(vcf_fp) + + self.assertEqual(counts['dragen'], 5) + self.assertEqual(counts['filter_pass'], 3) + # dragen(5) > MAX(4) → hypermutated, even though filter_pass(3) ≤ MAX(4) + self.assertTrue(counts['is_hypermutated']) + + def test_is_hypermutated_false_when_dragen_within_limit(self): + """is_hypermutated=False when dragen count ≤ MAX_SOMATIC_VARIANTS.""" + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'test.vcf' + rows = [(i * 10, 'PASS', '.') for i in range(1, 4)] + self._write_count_vcf(vcf_fp, rows) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 10): + counts = report_mod.count_variant_process(vcf_fp) + + self.assertEqual(counts['dragen'], 3) + self.assertFalse(counts['is_hypermutated']) + + def test_sage_novel_excluded_from_dragen_count(self): + """SAGE_NOVEL variants are not counted as DRAGEN variants.""" + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'test.vcf' + sage_novel_info = constants.VcfInfo.SAGE_NOVEL.value + rows = [ + (10, 'PASS', '.'), # dragen + (20, 'PASS', sage_novel_info), # sage novel — not dragen + (30, 'PASS', '.'), # dragen + ] + self._write_count_vcf(vcf_fp, rows) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 100): + counts = report_mod.count_variant_process(vcf_fp) + + self.assertEqual(counts['dragen'], 2) + self.assertEqual(counts['sage'], 3) + + def test_annotation_filter_excluded_from_annotated_count(self): + """Variants with bolt annotation filters are excluded from annotated count.""" + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'test.vcf' + annot_filter = constants.VcfFilter.MAX_VARIANTS_NON_PASS.value + rows = [ + (10, 'PASS', '.'), # annotated + (20, annot_filter, '.'), # not annotated (bolt annotation filter) + ] + self._write_count_vcf(vcf_fp, rows) + + with patch('bolt.common.constants.MAX_SOMATIC_VARIANTS', 100): + counts = report_mod.count_variant_process(vcf_fp) + + self.assertEqual(counts['annotated'], 1) + self.assertEqual(counts['dragen'], 2) + + +_PASS_COUNTS = {'pass': {'snps': 0, 'indels': 0, 'others': 0, 'total': 0}} + + +def _cli_args(dummy, output_dir): + """Return CliRunner args list for report entry(); all file paths point to dummy.""" + d = str(dummy) + return [ + '--tumor_name', 'TUMOR', + '--normal_name', 'NORMAL', + '--vcf_fp', d, + '--vcf_filters_fp', d, + '--vcf_dragen_fp', d, + '--vep_dir', str(dummy.parent), + '--purple_purity_fp', d, + '--cancer_genes_fp', d, + '--giab_regions_fp', d, + '--genome_fp', d, + '--threads', '1', + '--output_dir', str(output_dir), + ] + + +class TestEntrySkipsPcgrOnOverflow(unittest.TestCase): + """entry() catches RuntimeError from select_pcgr_variants and skips PCGR entirely.""" + + def test_run_somatic_not_called_when_cap_exceeded(self): + """When select_pcgr_variants raises RuntimeError, run_somatic must not be called.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + dummy = tmp_path / 'dummy.vcf.gz' + dummy.touch() + + with patch.object(report_mod, 'bcftools_stats_prepare', return_value=dummy), \ + patch.object(report_mod, 'run_bcftools_stats'), \ + patch.object(report_mod, 'allele_frequencies'), \ + patch.object(report_mod, 'count_variant_types', return_value=_PASS_COUNTS), \ + patch.object(report_mod, 'count_variant_process', + return_value={'filter_pass': constants.MAX_SOMATIC_VARIANTS + 1}), \ + patch.object(report_mod, 'parse_purple_purity_file', + return_value={'purity': 0.8, 'ploidy': 2.0}), \ + patch.object(report_mod, 'select_pcgr_variants', + side_effect=RuntimeError('595416 > 450000')), \ + patch.object(pcgr, 'prepare_vcf_somatic') as mock_prep, \ + patch.object(pcgr, 'run_somatic') as mock_run: + result = CliRunner().invoke(report_mod.entry, _cli_args(dummy, tmp_path / 'out')) + + self.assertEqual(result.exit_code, 0, result.output) + mock_prep.assert_not_called() + mock_run.assert_not_called() + + def test_run_somatic_called_when_within_limit(self): + """When PASS count is within limit, run_somatic must be called normally.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + dummy = tmp_path / 'dummy.vcf.gz' + dummy.touch() + fake_prep_output = tmp_path / 'prep.vcf.gz' + fake_prep_output.touch() + + with patch.object(report_mod, 'bcftools_stats_prepare', return_value=dummy), \ + patch.object(report_mod, 'run_bcftools_stats'), \ + patch.object(report_mod, 'allele_frequencies'), \ + patch.object(report_mod, 'count_variant_types', return_value=_PASS_COUNTS), \ + patch.object(report_mod, 'count_variant_process', + return_value={'filter_pass': constants.MAX_SOMATIC_VARIANTS - 1}), \ + patch.object(report_mod, 'parse_purple_purity_file', + return_value={'purity': 0.8, 'ploidy': 2.0}), \ + patch.object(pcgr, 'prepare_vcf_somatic', return_value=fake_prep_output), \ + patch.object(pcgr, 'run_somatic') as mock_run: + result = CliRunner().invoke(report_mod.entry, _cli_args(dummy, tmp_path / 'out')) + + self.assertEqual(result.exit_code, 0, result.output) + mock_run.assert_called_once() + self.assertIs(mock_run.call_args.kwargs['disable_estimates'], False) + + def test_run_somatic_called_with_disable_estimates_when_trimmed(self): + """When select_pcgr_variants trims a hypermutated sample, estimates must be disabled.""" + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + dummy = tmp_path / 'dummy.vcf.gz' + dummy.touch() + fake_selected_output = tmp_path / 'selected.vcf.gz' + fake_selected_output.touch() + fake_prep_output = tmp_path / 'prep.vcf.gz' + fake_prep_output.touch() + + with patch.object(report_mod, 'bcftools_stats_prepare', return_value=dummy), \ + patch.object(report_mod, 'run_bcftools_stats'), \ + patch.object(report_mod, 'allele_frequencies'), \ + patch.object(report_mod, 'count_variant_types', return_value=_PASS_COUNTS), \ + patch.object(report_mod, 'count_variant_process', + return_value={'filter_pass': constants.MAX_SOMATIC_VARIANTS + 1}), \ + patch.object(report_mod, 'parse_purple_purity_file', + return_value={'purity': 0.8, 'ploidy': 2.0}), \ + patch.object(report_mod, 'select_pcgr_variants', + return_value=fake_selected_output), \ + patch.object(pcgr, 'prepare_vcf_somatic', return_value=fake_prep_output), \ + patch.object(pcgr, 'run_somatic') as mock_run: + result = CliRunner().invoke(report_mod.entry, _cli_args(dummy, tmp_path / 'out')) + + self.assertEqual(result.exit_code, 0, result.output) + mock_run.assert_called_once() + self.assertIs(mock_run.call_args.kwargs['disable_estimates'], True) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_smlv_somatic_rescue.py b/tests/test_smlv_somatic_rescue.py new file mode 100644 index 0000000..cd66902 --- /dev/null +++ b/tests/test_smlv_somatic_rescue.py @@ -0,0 +1,75 @@ +"""Unit tests for bolt/workflows/smlv_somatic/rescue.py. + +Covers only the SAGE VCF header-consistency check in +annotate_existing_sage_calls(), which is reachable without a bcftools binary +(the check runs before any subprocess call — a SystemExit from a mismatch +short-circuits execution before bcftools is ever invoked). The rest of +rescue.py orchestrates bcftools subprocesses end-to-end and is not +unit-tested here. +""" +import pathlib +import tempfile +import unittest + +import bolt.common.constants as constants +import bolt.workflows.smlv_somatic.rescue as rescue + + +def _write_input_vcf_no_sage_fields(path): + """A plain DRAGEN-style VCF with no SAGE_* header lines at all.""" + with open(path, 'w') as fh: + fh.write( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + '##FORMAT=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTUMOR\n' + 'chr1\t100\t.\tA\tT\t.\tPASS\t.\tAD\t10,5\n' + ) + + +def _write_sage_vcf_mismatched_hotspot_description(path): + """A SAGE-style VCF whose SAGE_HOTSPOT description does not match constants.py.""" + with open(path, 'w') as fh: + fh.write( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + '##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + 'chr1\t100\t.\tA\tT\t.\tPASS\tSAGE_HOTSPOT\n' + ) + + +class TestAnnotateExistingSageCallsHeaderCheck(unittest.TestCase): + """Regression test for rescue.py's SAGE VCF header-consistency check. + + annotate_existing_sage_calls(input_fp, tumor_name, sage_vcf_fp, output_dir) + must validate header descriptions against sage_vcf_fp (the SAGE VCF), not + input_fp (the DRAGEN VCF being annotated). Prior to the fix, the check was + wired to input_fp, which has no SAGE_* header lines at all — so a mismatch + in the real SAGE VCF's headers went undetected. + """ + + def test_header_mismatch_in_sage_vcf_triggers_exit(self): + """A mismatched SAGE_HOTSPOT description in sage_vcf_fp must raise SystemExit. + + input_fp has no SAGE_* headers at all, so if the check were (incorrectly) + run against input_fp instead, every field would be silently skipped and + no SystemExit would be raised — proving the check now reads sage_vcf_fp. + """ + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + input_fp = tmp_path / 'input.vcf' + sage_vcf_fp = tmp_path / 'sage.vcf' + _write_input_vcf_no_sage_fields(input_fp) + _write_sage_vcf_mismatched_hotspot_description(sage_vcf_fp) + + with self.assertRaises(SystemExit): + rescue.annotate_existing_sage_calls( + input_fp, 'TUMOR', sage_vcf_fp, tmp_path, + ) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_util.py b/tests/test_util.py new file mode 100644 index 0000000..bdd4ab0 --- /dev/null +++ b/tests/test_util.py @@ -0,0 +1,337 @@ +"""Unit tests for bolt/util.py — VCF header helpers and merge helpers. + +Covers binary-free logic plus a bcftools-guarded integration test for +merge_vcf_files (see TestMergeVcfFiles). The remaining bcftools-dependent +functions (count_vcf_records, execute_command) are intentionally NOT tested +here. +""" +import gzip +import pathlib +import shutil +import subprocess +import tempfile +import unittest + +import cyvcf2 + +import bolt.common.constants as constants +import bolt.util as util + + +# Minimal VCF header covering enum members exercised below +HEADER = ( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + f'##FILTER=\n' + f'##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' +) + + +def _write_vcf(path, variants): + with open(path, 'w') as fh: + fh.write(HEADER) + for pos, info in variants: + fh.write(f'chr1\t{pos}\t.\tA\tT\t.\tPASS\t{info}\n') + + +def _make_vcf_handle(tmp_dir): + """Return a cyvcf2.VCF handle opened on a minimal VCF written to tmp_dir.""" + vcf_path = pathlib.Path(tmp_dir) / 'test.vcf' + _write_vcf(vcf_path, [(100, '.')]) + return cyvcf2.VCF(str(vcf_path)) + + +class TestGetVcfHeaderEntry(unittest.TestCase): + """Unit tests for util.get_vcf_header_entry().""" + + def test_id_matches_enum_value(self): + entry = util.get_vcf_header_entry(constants.VcfInfo.SAGE_NOVEL) + self.assertEqual(entry['ID'], constants.VcfInfo.SAGE_NOVEL.value) + + def test_merges_header_entries_fields(self): + entry = util.get_vcf_header_entry(constants.VcfInfo.SAGE_NOVEL) + expected = constants.VCF_HEADER_ENTRIES[constants.VcfInfo.SAGE_NOVEL] + for key, value in expected.items(): + self.assertEqual(entry[key], value) + + def test_filter_enum_entry(self): + entry = util.get_vcf_header_entry(constants.VcfFilter.MIN_AF) + self.assertEqual(entry['ID'], constants.VcfFilter.MIN_AF.value) + self.assertIn('Description', entry) + + +class TestGetVcfHeaderLine(unittest.TestCase): + """Unit tests for util.get_vcf_header_line().""" + + def test_filter_line_format(self): + line = util.get_vcf_header_line(constants.VcfFilter.MIN_AF) + self.assertTrue(line.startswith('##FILTER=<')) + self.assertIn(f'ID={constants.VcfFilter.MIN_AF.value}', line) + + def test_info_line_format(self): + line = util.get_vcf_header_line(constants.VcfInfo.SAGE_NOVEL) + self.assertTrue(line.startswith('##INFO=<')) + self.assertIn(f'ID={constants.VcfInfo.SAGE_NOVEL.value}', line) + self.assertIn('Number=', line) + self.assertIn('Type=', line) + + def test_format_line_format(self): + line = util.get_vcf_header_line(constants.VcfFormat.SAGE_AD) + self.assertTrue(line.startswith('##FORMAT=<')) + self.assertIn(f'ID={constants.VcfFormat.SAGE_AD.value}', line) + + +class TestGetQualifiedVcfAnnotation(unittest.TestCase): + """Unit tests for util.get_qualified_vcf_annotation().""" + + def test_info_enum_qualified(self): + result = util.get_qualified_vcf_annotation(constants.VcfInfo.SAGE_NOVEL) + self.assertEqual(result, f'INFO/{constants.VcfInfo.SAGE_NOVEL.value}') + + def test_format_enum_qualified(self): + result = util.get_qualified_vcf_annotation(constants.VcfFormat.SAGE_AD) + self.assertEqual(result, f'FORMAT/{constants.VcfFormat.SAGE_AD.value}') + + def test_filter_enum_raises(self): + with self.assertRaises(AssertionError): + util.get_qualified_vcf_annotation(constants.VcfFilter.MIN_AF) + + +class TestAddVcfHeaderEntry(unittest.TestCase): + """Unit tests for util.add_vcf_header_entry().""" + + def test_info_tag_added_to_header(self): + with tempfile.TemporaryDirectory() as tmp: + fh = _make_vcf_handle(tmp) + util.add_vcf_header_entry(fh, constants.VcfInfo.PCGR_ACTIONABILITY_TIER) + header_type = fh.get_header_type(constants.VcfInfo.PCGR_ACTIONABILITY_TIER.value) + self.assertEqual(header_type['ID'], constants.VcfInfo.PCGR_ACTIONABILITY_TIER.value) + + def test_filter_tag_added_to_header(self): + with tempfile.TemporaryDirectory() as tmp: + fh = _make_vcf_handle(tmp) + util.add_vcf_header_entry(fh, constants.VcfFilter.PON) + # FILTER header lines live under BCF_HL_FLT (order=0); get_header_type + # defaults to [INFO, FORMAT] (order=[1, 2]) so FILTER lookups need order=[0]. + header_type = fh.get_header_type(constants.VcfFilter.PON.value, order=[0]) + self.assertEqual(header_type['ID'], constants.VcfFilter.PON.value) + + def test_format_tag_added_to_header(self): + with tempfile.TemporaryDirectory() as tmp: + fh = _make_vcf_handle(tmp) + util.add_vcf_header_entry(fh, constants.VcfFormat.SAGE_DP) + header_type = fh.get_header_type(constants.VcfFormat.SAGE_DP.value) + self.assertEqual(header_type['ID'], constants.VcfFormat.SAGE_DP.value) + + +class TestMergeTsvFiles(unittest.TestCase): + """Unit tests for util.merge_tsv_files().""" + + def _write_gz_tsv(self, path, rows): + with gzip.open(path, 'wt', encoding='utf-8') as fh: + for row in rows: + fh.write(row + '\n') + + def test_single_header_row_retained(self): + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + tsv1 = tmp_path / 'a.tsv.gz' + tsv2 = tmp_path / 'b.tsv.gz' + self._write_gz_tsv(tsv1, ['CHROM\tPOS', 'chr1\t100']) + self._write_gz_tsv(tsv2, ['CHROM\tPOS', 'chr1\t200']) + + merged_fp = tmp_path / 'merged.tsv.gz' + util.merge_tsv_files([tsv1, tsv2], merged_fp) + + with gzip.open(merged_fp, 'rt', encoding='utf-8') as fh: + lines = [line.rstrip('\n') for line in fh] + + self.assertEqual(lines.count('CHROM\tPOS'), 1) + + def test_all_data_rows_present_in_order(self): + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + tsv1 = tmp_path / 'a.tsv.gz' + tsv2 = tmp_path / 'b.tsv.gz' + self._write_gz_tsv(tsv1, ['CHROM\tPOS', 'chr1\t100', 'chr1\t101']) + self._write_gz_tsv(tsv2, ['CHROM\tPOS', 'chr1\t200']) + + merged_fp = tmp_path / 'merged.tsv.gz' + util.merge_tsv_files([tsv1, tsv2], merged_fp) + + with gzip.open(merged_fp, 'rt', encoding='utf-8') as fh: + lines = [line.rstrip('\n') for line in fh] + + self.assertEqual( + lines, + ['CHROM\tPOS', 'chr1\t100', 'chr1\t101', 'chr1\t200'], + ) + + def test_output_is_gzipped(self): + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + tsv1 = tmp_path / 'a.tsv.gz' + self._write_gz_tsv(tsv1, ['CHROM\tPOS', 'chr1\t100']) + + merged_fp = tmp_path / 'merged.tsv.gz' + util.merge_tsv_files([tsv1], merged_fp) + + # gzip files start with the magic number 0x1f 0x8b + with open(merged_fp, 'rb') as fh: + magic = fh.read(2) + self.assertEqual(magic, b'\x1f\x8b') + + +class TestCheckAnnotationHeaders(unittest.TestCase): + """Unit tests for util.check_annotation_headers().""" + + def _write_vcf_with_sage_hotspot(self, path, description): + with open(path, 'w') as fh: + fh.write( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + f'##INFO=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + 'chr1\t100\t.\tA\tT\t.\tPASS\t.\n' + ) + + def _write_vcf_without_sage_fields(self, path): + with open(path, 'w') as fh: + fh.write( + '##fileformat=VCFv4.2\n' + '##FILTER=\n' + '##contig=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + 'chr1\t100\t.\tA\tT\t.\tPASS\t.\n' + ) + + def test_matching_description_returns_normally(self): + expected = constants.VCF_HEADER_ENTRIES[constants.VcfInfo.SAGE_HOTSPOT]['Description'] + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'matching.vcf' + self._write_vcf_with_sage_hotspot(vcf_fp, expected) + # No exception/SystemExit raised + util.check_annotation_headers( + {constants.VcfInfo.SAGE_HOTSPOT: 'SAGE_HOTSPOT'}, vcf_fp, + ) + + def test_mismatched_description_exits(self): + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'mismatched.vcf' + self._write_vcf_with_sage_hotspot(vcf_fp, 'a totally different description') + with self.assertRaises(SystemExit): + util.check_annotation_headers( + {constants.VcfInfo.SAGE_HOTSPOT: 'SAGE_HOTSPOT'}, vcf_fp, + ) + + def test_field_absent_from_target_vcf_is_skipped(self): + with tempfile.TemporaryDirectory() as tmp: + vcf_fp = pathlib.Path(tmp) / 'no_sage.vcf' + self._write_vcf_without_sage_fields(vcf_fp) + # SAGE_HOTSPOT has no header entry in this VCF at all; must be + # skipped rather than raising, so no exception/SystemExit here. + util.check_annotation_headers( + {constants.VcfInfo.SAGE_HOTSPOT: 'SAGE_HOTSPOT'}, vcf_fp, + ) + + +@unittest.skipUnless(shutil.which('bcftools'), 'bcftools not available') +class TestMergeVcfFiles(unittest.TestCase): + """Integration tests for util.merge_vcf_files(). + + merge_vcf_files reassembles PCGR hypermutated chunk outputs with + `bcftools merge -m all`. These chunks are sites-only VCFs (no FORMAT or + sample columns) because pcgr.prepare_vcf_somatic / get_minimal_header strip + them. bcftools merge only fails with "Duplicate sample names" when inputs + carry a same-named genotype column; on sites-only inputs it correctly + produces the union. These tests lock that invariant in: any regression that + reintroduces a sample column (which would break the merge) is caught here. + """ + + # Sites-only header (no FORMAT, no sample column) — mirrors get_minimal_header + _HEADER = ( + '##fileformat=VCFv4.2\n' + '##contig=\n' + '##contig=\n' + '##INFO=\n' + '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + ) + + def _write_chunk(self, path, records): + """Write a sites-only VCF, then bgzip + tabix-index it via bcftools. + + `records` is an iterable of (chrom, pos, ref, alt) tuples. Positions + within a chunk are sorted before writing so indexing succeeds; chunks + may be mutually out of order to exercise the cross-chunk sort. + """ + plain = pathlib.Path(f'{path}.plain.vcf') + with open(plain, 'w') as fh: + fh.write(self._HEADER) + for chrom, pos, ref, alt in sorted(records, key=lambda r: (r[0], r[1])): + fh.write(f'{chrom}\t{pos}\t.\t{ref}\t{alt}\t.\tPASS\tPCGR_TIER=1\n') + subprocess.run(['bcftools', 'view', '-Oz', '-o', str(path), str(plain)], check=True) + subprocess.run(['bcftools', 'index', '-t', str(path)], check=True) + return records + + def _read_keys(self, vcf_fp): + return [ + (record.CHROM, record.POS, record.REF, record.ALT[0]) + for record in cyvcf2.VCF(str(vcf_fp)) + ] + + def test_merge_is_lossless_and_sorted(self): + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + # Disjoint positions, and chunks deliberately out of order relative + # to each other so the merge must interleave/sort across chunks. + chunk_a = self._write_chunk( + tmp_path / 'chunk_a.vcf.gz', + [('chr1', 300, 'G', 'A'), ('chr1', 100, 'A', 'T'), ('chr2', 50, 'C', 'G')], + ) + chunk_b = self._write_chunk( + tmp_path / 'chunk_b.vcf.gz', + [('chr1', 200, 'C', 'G'), ('chr1', 400, 'T', 'C')], + ) + + merged_vcf = util.merge_vcf_files( + [tmp_path / 'chunk_a.vcf.gz', tmp_path / 'chunk_b.vcf.gz'], + tmp_path / 'merged.pass', + ) + + merged_keys = self._read_keys(merged_vcf) + expected_keys = list(chunk_a) + list(chunk_b) + + # No loss, no duplication: exact multiset match + self.assertEqual(len(merged_keys), len(expected_keys)) + self.assertCountEqual(merged_keys, expected_keys) + # Position-sorted output (cross-chunk interleave) + self.assertEqual( + merged_keys, + sorted(merged_keys, key=lambda k: (k[0], k[1])), + ) + + def test_merge_output_is_bgzipped_and_indexed(self): + with tempfile.TemporaryDirectory() as tmp: + tmp_path = pathlib.Path(tmp) + self._write_chunk(tmp_path / 'chunk_a.vcf.gz', [('chr1', 100, 'A', 'T')]) + self._write_chunk(tmp_path / 'chunk_b.vcf.gz', [('chr1', 200, 'C', 'G')]) + + merged_vcf = util.merge_vcf_files( + [tmp_path / 'chunk_a.vcf.gz', tmp_path / 'chunk_b.vcf.gz'], + tmp_path / 'merged.pass', + ) + + # BGZF/gzip magic bytes + with open(merged_vcf, 'rb') as fh: + self.assertEqual(fh.read(2), b'\x1f\x8b') + # merge_vcf_files tabix-indexes its output + self.assertTrue(pathlib.Path(f'{merged_vcf}.tbi').exists()) + # Intermediate unsorted file is cleaned up + self.assertFalse((tmp_path / 'merged.pass.unsorted.vcf.gz').exists()) + + +if __name__ == '__main__': + unittest.main()