diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 925ef24..d86a437 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 0.2.18 +current_version = 0.3.1 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 cd8796a..22018a8 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -17,6 +17,22 @@ jobs: - name: Checkout repository uses: actions/checkout@v4 + - name: Free Disk Space (Ubuntu) + uses: jlumbroso/free-disk-space@main + with: + # this might remove tools that are actually needed, + # if set to "true" but frees about 6 GB + tool-cache: false + + # all of these default to true, but feel free to set to + # "false" if necessary for your workflow + android: true + dotnet: true + haskell: true + large-packages: true + docker-images: true + swap-storage: true + - name: Login to GitHub Container Registry uses: docker/login-action@v3 with: @@ -27,7 +43,8 @@ jobs: - name: Set up Docker Buildx uses: docker/setup-buildx-action@v3 - - name: Build and push image + - name: Build image + id: build run: | if [[ "${{ matrix.build_name }}" == "base" ]]; then docker_tag=bolt:${GITHUB_REF_NAME#v}; @@ -37,10 +54,26 @@ jobs: dockerfile_fn=Dockerfile.${{ matrix.build_name }}; fi; - # Build and push to GHCR umccr registry docker build \ --platform linux/amd64 \ + --load \ -f docker/${dockerfile_fn} \ -t ghcr.io/umccr/${docker_tag} \ - --push \ .; + + echo "docker_tag=ghcr.io/umccr/${docker_tag}" >> $GITHUB_OUTPUT + + - name: Smoke test + run: | + docker run --rm ${{ steps.build.outputs.docker_tag }} bolt --version + + case "${{ matrix.build_name }}" in + multiqc) docker run --rm ${{ steps.build.outputs.docker_tag }} multiqc --version ;; + gpgr) docker run --rm ${{ steps.build.outputs.docker_tag }} Rscript -e 'library(gpgr); cat("gpgr ok\n")' ;; + pcgr) docker run --rm ${{ steps.build.outputs.docker_tag }} conda run --no-capture-output -n pcgr pcgr --version ;; + snpeff) docker run --rm ${{ steps.build.outputs.docker_tag }} bash -c 'snpEff -version 2>&1 | head -1' ;; + circos) docker run --rm ${{ steps.build.outputs.docker_tag }} circos --version ;; + esac + + - name: Push image + run: docker push ${{ steps.build.outputs.docker_tag }} diff --git a/.gitignore b/.gitignore index 59a0c10..e4ef844 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,5 @@ __pycache__/ build/ venv/ working/ +data/ +workspace/ \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index a519ab0..e4a5ea6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,15 @@ # bolt changelog -## dev +## 0.3.1 + +- Fix `ModuleNotFoundError: No module named 'pkg_resources'` in `bolt:0.3.0-multiqc` — add `setuptools <81` to conda env +- Fix `merge_vcf_files` producing wrong output filename — `Path.with_suffix()` was stripping `.pass` component; use explicit path concatenation instead +- 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`) + +## 0.3.0 - [28](https://github.com/umccr/bolt/pull/28) - gpgr version bump to 2.2.12 for cancer report hypermutated flag fix @@ -10,4 +19,6 @@ - [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 uptade, with switch sv caller from GRIPSS to eSVee \ No newline at end of file +- [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/bolt/common/constants.py b/bolt/common/constants.py index e38ed3d..8c88b93 100644 --- a/bolt/common/constants.py +++ b/bolt/common/constants.py @@ -4,7 +4,12 @@ ###################################### ## Variation selection (annotation) ## ###################################### -MAX_SOMATIC_VARIANTS = 500_000 + +# Cap below PCGR's 500k limit. PCGR silently drops variants or skips HTML +# generation above 500k, and its multi-allelic decomposition can inflate +# variant count beyond what bolt outputs. 50k margin absorbs this safely. +# See: docs/adr/001-max-somatic-variants-450k.md +MAX_SOMATIC_VARIANTS = 450_000 MAX_SOMATIC_VARIANTS_GNOMAD_FILTER = 0.01 @@ -35,12 +40,51 @@ 'pathogenic', 'uncertain_significance', } -PCGR_TIERS_RESCUE = { - 'TIER_1', - 'TIER_2', +PCGR_ACTIONABILITY_TIER_RESCUE = { + '1', + '2', } +################################ +## Hypermutated report filter ## +################################ +# Values match short forms written by transfer_annotations_somatic() ('1'=TIER_1, ..., 'N'=NONCODING). +# Order is lowest clinical priority first: NONCODING dropped before TIER_1. +PCGR_TIERS_FILTERING = ( + 'N', + '4', + '3', + '2', + '1', +) + +VEP_IMPACTS_FILTER = ( + 'intergenic', + 'intronic', + 'downstream', + 'upstream', + 'impacts_other', +) + +GENOMIC_REGIONS_FILTERING = ( + 'difficult', + 'none', + 'giab_conf', +) + +HOTSPOT_FIELDS_FILTERING = ( + 'SAGE_HOTSPOT', + 'hotspot', + 'PCGR_MUTATION_HOTSPOT', +) + +RETAIN_FIELDS_FILTERING = ( + 'PANEL', + *HOTSPOT_FIELDS_FILTERING, +) + + ################################################## ## VCF FILTER tags and FORMAT, INFO annotations ## ################################################## @@ -61,6 +105,8 @@ class VcfFilter(enum.Enum): ENCODE = 'ENCODE' GNOMAD_COMMON = 'gnomAD_common' + PCGR_COUNT_LIMIT = 'PCGR_count_limit' + @property def namespace(self): return 'FILTER' @@ -77,16 +123,14 @@ class VcfInfo(enum.Enum): SAGE_NOVEL = 'SAGE_NOVEL' SAGE_RESCUE = 'SAGE_RESCUE' - PCGR_TIER = 'PCGR_TIER' + PCGR_ACTIONABILITY_TIER = 'PCGR_ACTIONABILITY_TIER' PCGR_CSQ = 'PCGR_CSQ' PCGR_MUTATION_HOTSPOT = 'PCGR_MUTATION_HOTSPOT' - PCGR_CLINVAR_CLNSIG = 'PCGR_CLINVAR_CLNSIG' + PCGR_CLINVAR_CLASSIFICATION = 'PCGR_CLINVAR_CLASSIFICATION' PCGR_COSMIC_COUNT = 'PCGR_COSMIC_COUNT' PCGR_TCGA_PANCANCER_COUNT = 'PCGR_TCGA_PANCANCER_COUNT' PCGR_ICGC_PCAWG_COUNT = 'PCGR_ICGC_PCAWG_COUNT' - CPSR_FINAL_CLASSIFICATION = 'CPSR_FINAL_CLASSIFICATION' - CPSR_PATHOGENICITY_SCORE = 'CPSR_PATHOGENICITY_SCORE' CPSR_CLINVAR_CLASSIFICATION = 'CPSR_CLINVAR_CLASSIFICATION' CPSR_CSQ = 'CPSR_CSQ' @@ -112,7 +156,7 @@ class VcfInfo(enum.Enum): GNOMAD_AF = 'gnomAD_AF' - PCGR_TIER_RESCUE = 'PCGR_TIER_RESCUE' + PCGR_ACTIONABILITY_TIER_RESCUE = 'PCGR_ACTIONABILITY_TIER_RESCUE' SAGE_HOTSPOT_RESCUE = 'SAGE_HOTSPOT_RESCUE' CLINICAL_POTENTIAL_RESCUE = 'CLINICAL_POTENTIAL_RESCUE' @@ -121,6 +165,8 @@ class VcfInfo(enum.Enum): RESCUED_FILTERS_EXISTING = 'RESCUED_FILTERS_EXISTING' RESCUED_FILTERS_PENDING = 'RESCUED_FILTERS_PENDING' + PANEL = 'PANEL' + @property def namespace(self): return 'INFO' @@ -187,6 +233,12 @@ def namespace(self): 'Description': f'gnomAD AF >= {MAX_GNOMAD_AF}', }, + VcfFilter.PCGR_COUNT_LIMIT: { + 'Description': ( + f'Manually filtered to {MAX_SOMATIC_VARIANTS} variants to stay below the ' + f'PCGR hard limit of 500,000 variants' + ), + }, # INFO VcfInfo.TUMOR_AF: { @@ -226,39 +278,42 @@ def namespace(self): 'Description': 'Variant rescued by a matching SAGE call', }, - VcfInfo.PCGR_TIER: { + VcfInfo.PCGR_ACTIONABILITY_TIER: { 'Number': '1', 'Type': 'String', 'Description': ( - 'Tier reported by PCGR with the following meaning: TIER_1: strong clinical ' - 'significance; TIER_2: potential clinical significance; TIER_3: uncertain clinical ' - 'significance; TIER_4: other coding variants; NONCODING: other non-coding variants' + 'Tier reported by PCGR: 1: strong clinical significance; ' + '2: potential clinical significance; 3: uncertain clinical significance; ' + '4: other coding variants; N: other non-coding variants' ), }, VcfInfo.PCGR_CSQ: { 'Number': '.', 'Type': 'String', 'Description': ( - 'Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|' - 'Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|' - 'CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|ALLELE_NUM|' - 'DISTANCE|STRAND|FLAGS|PICK|VARIANT_CLASS|SYMBOL_SOURCE|HGNC_ID|CANONICAL|' - 'MANE_SELECT|MANE_PLUS_CLINICAL|TSL|APPRIS|CCDS|ENSP|SWISSPROT|TREMBL|UNIPARC|' - 'UNIPROT_ISOFORM|RefSeq|DOMAINS|HGVS_OFFSET|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|' - 'gnomAD_AF|gnomAD_AFR_AF|gnomAD_AMR_AF|gnomAD_ASJ_AF|gnomAD_EAS_AF|gnomAD_FIN_AF|' - 'gnomAD_NFE_AF|gnomAD_OTH_AF|gnomAD_SAS_AF|CLIN_SIG|SOMATIC|PHENO|CHECK_REF|' - 'NearestExonJB' + 'Consequence annotations from Ensembl VEP. Format: ' + 'Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|' + 'HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|' + 'ALLELE_NUM|DISTANCE|STRAND|FLAGS|PICK|VARIANT_CLASS|SYMBOL_SOURCE|HGNC_ID|CANONICAL|' + 'MANE|MANE_SELECT|MANE_PLUS_CLINICAL|TSL|APPRIS|CCDS|ENSP|SWISSPROT|TREMBL|UNIPARC|' + 'UNIPROT_ISOFORM|RefSeq|DOMAINS|HGVS_OFFSET|gnomADe_AF|gnomADe_AFR_AF|gnomADe_AMR_AF|' + 'gnomADe_ASJ_AF|gnomADe_EAS_AF|gnomADe_FIN_AF|gnomADe_MID_AF|gnomADe_NFE_AF|' + 'gnomADe_REMAINING_AF|gnomADe_SAS_AF|gnomADg_AF|gnomADg_AFR_AF|gnomADg_AMI_AF|' + 'gnomADg_AMR_AF|gnomADg_ASJ_AF|gnomADg_EAS_AF|gnomADg_FIN_AF|gnomADg_MID_AF|' + 'gnomADg_NFE_AF|gnomADg_REMAINING_AF|gnomADg_SAS_AF|CLIN_SIG|SOMATIC|PHENO|CHECK_REF|' + 'MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE|TRANSCRIPTION_FACTORS|NearestExonJB|' + 'MaxEntScan_alt|MaxEntScan_diff|MaxEntScan_ref' ), }, VcfInfo.PCGR_MUTATION_HOTSPOT: { 'Number': '.', 'Type': 'String', - 'Description': 'Known cancer mutation hotspot, as found in cancerhotspots.org_v2, Gene|Codon|Q-value', + 'Description': 'Known cancer mutation hotspot, as found in cancerhotspots.org. Format: GeneSymbol|Entrez_ID|CodonRefAA|Alt_AA|Q-value', }, - VcfInfo.PCGR_CLINVAR_CLNSIG: { + VcfInfo.PCGR_CLINVAR_CLASSIFICATION: { 'Number': '.', 'Type': 'String', - 'Description': 'ClinVar clinical significance', + 'Description': 'ClinVar - Overall clinical significance of variant on a five-tiered scale', }, VcfInfo.PCGR_COSMIC_COUNT: { 'Number': '1', @@ -270,44 +325,27 @@ def namespace(self): 'Type': 'Integer', 'Description': 'Raw variant count across all tumor types', }, + VcfInfo.HMF_HOTSPOT: { + 'Number': '0', + 'Type': 'Flag', + 'Description': 'calculated by flag of overlapping values in field HMF from annotations/hotspots/hotspots.hmf.vcf.gz', + }, VcfInfo.PCGR_ICGC_PCAWG_COUNT: { 'Number': '1', 'Type': 'Integer', 'Description': 'Count of ICGC PCAWG hits', }, - VcfInfo.CPSR_FINAL_CLASSIFICATION: { - 'Number': '1', - 'Type': 'String', - 'Description': ( - 'Final variant classification based on the combination of CLINVAR_CLASSIFICTION (for ' - 'ClinVar-classified variants), and CPSR_CLASSIFICATION (for novel variants)' - ), - }, - VcfInfo.CPSR_PATHOGENICITY_SCORE: { - 'Number': '1', - 'Type': 'Float', - 'Description': 'Aggregated CPSR pathogenicity score', - }, VcfInfo.CPSR_CLINVAR_CLASSIFICATION: { 'Number': '1', 'Type': 'String', - 'Description': 'Clinical significance of variant on a five-tiered scale', + 'Description': 'ClinVar - Overall clinical significance of variant on a five-tiered scale', }, VcfInfo.CPSR_CSQ: { 'Number': '.', 'Type': 'String', 'Description': ( - 'Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|' - 'Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|' - 'Protein_position|Amino_acids|Codons|Existing_variation|ALLELE_NUM|DISTANCE|STRAND|' - 'FLAGS|PICK|VARIANT_CLASS|SYMBOL_SOURCE|HGNC_ID|CANONICAL|MANE_SELECT|' - 'MANE_PLUS_CLINICAL|APPRIS|CCDS|ENSP|SWISSPROT|TREMBL|UNIPARC|UNIPROT_ISOFORM|RefSeq|' - 'DOMAINS|HGVS_OFFSET|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|gnomAD_AF|gnomAD_AFR_AF|' - 'gnomAD_AMR_AF|gnomAD_ASJ_AF|gnomAD_EAS_AF|gnomAD_FIN_AF|gnomAD_NFE_AF|gnomAD_OTH_AF|' - 'gnomAD_SAS_AF|CLIN_SIG|SOMATIC|PHENO|CHECK_REF|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|' - 'MOTIF_SCORE_CHANGE|TRANSCRIPTION_FACTORS|NearestExonJB|LoF|LoF_filter|LoF_flags|' - 'LoF_info' + 'Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|ALLELE_NUM|DISTANCE|STRAND|FLAGS|PICK|VARIANT_CLASS|SYMBOL_SOURCE|HGNC_ID|CANONICAL|MANE|MANE_SELECT|MANE_PLUS_CLINICAL|TSL|APPRIS|CCDS|ENSP|SWISSPROT|TREMBL|UNIPARC|UNIPROT_ISOFORM|RefSeq|DOMAINS|HGVS_OFFSET|gnomADe_AF|gnomADe_AFR_AF|gnomADe_AMR_AF|gnomADe_ASJ_AF|gnomADe_EAS_AF|gnomADe_FIN_AF|gnomADe_MID_AF|gnomADe_NFE_AF|gnomADe_REMAINING_AF|gnomADe_SAS_AF|gnomADg_AF|gnomADg_AFR_AF|gnomADg_AMI_AF|gnomADg_AMR_AF|gnomADg_ASJ_AF|gnomADg_EAS_AF|gnomADg_FIN_AF|gnomADg_MID_AF|gnomADg_NFE_AF|gnomADg_REMAINING_AF|gnomADg_SAS_AF|CLIN_SIG|SOMATIC|PHENO|CHECK_REF|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE|TRANSCRIPTION_FACTORS|NearestExonJB|MaxEntScan_alt|MaxEntScan_diff|MaxEntScan_ref' ), }, @@ -316,7 +354,7 @@ def namespace(self): 'Type': 'Flag', 'Description': '', }, - VcfInfo.PCGR_TIER_RESCUE: { + VcfInfo.PCGR_ACTIONABILITY_TIER_RESCUE: { 'Number': '0', 'Type': 'Flag', 'Description': '', @@ -345,11 +383,17 @@ def namespace(self): }, VcfInfo.RESCUED_FILTERS_PENDING: { - 'Number': '1', + 'Number': '.', 'Type': 'String', 'Description': 'Filters pending prior to variant rescue', }, + VcfInfo.PANEL: { + 'Number': '0', + 'Type': 'Flag', + 'Description': 'UMCCR somatic panel CDS (2,000 bp padding)', + }, + # FORMAT VcfFormat.SAGE_AD: { @@ -368,13 +412,11 @@ def namespace(self): 'Description': 'Approximate read depth (reads with MQ=255 or with bad mates are filtered)', }, VcfFormat.SAGE_SB: { - 'Number': '1', + 'Number': '2', 'Type': 'Float', - 'Description': 'Strand bias - percentage of first-in-pair reads', + 'Description': 'Fragment strand bias - percentage of forward-orientation fragments (ref,alt)', }, } - - ##################### ## Other ## ##################### diff --git a/bolt/common/pcgr.py b/bolt/common/pcgr.py index fce4c56..f1ce2cb 100644 --- a/bolt/common/pcgr.py +++ b/bolt/common/pcgr.py @@ -1,9 +1,12 @@ import csv +import functools +import gzip +import itertools import pathlib import re import shutil import tempfile - +import logging import cyvcf2 @@ -11,6 +14,8 @@ from .. import util from ..common import constants +# Use the existing logger configuration +logger = logging.getLogger(__name__) def prepare_vcf_somatic(input_fp, tumor_name, normal_name, output_dir): @@ -110,13 +115,22 @@ def get_minimal_header(input_fh): return '\n'.join([filetype_line, *chrom_lines, *format_lines, column_line]) -def run_somatic(input_fp, pcgr_refdata_dir, output_dir, threads=1, pcgr_conda=None, pcgrr_conda=None, purity=None, ploidy=None, sample_id=None): +def run_somatic(input_fp, pcgr_refdata_dir, vep_dir, output_dir, chunk_nbr=None, threads=1, pcgr_threads=4, pcgr_conda=None, pcgrr_conda=None, purity=None, ploidy=None, sample_id=None, disable_estimates=False): + # threads: Nextflow process-level resource allocation (not wired to PCGR internals) + # pcgr_threads: PCGR-internal concurrency (vcfanno workers + VEP forks) + + pcgr_threads = max(1, int(pcgr_threads)) + vcfanno_threads = pcgr_threads + vep_forks = min(8, max(2, pcgr_threads)) - # NOTE(SW): Nextflow FusionFS v2.2.8 does not support PCGR output to S3; instead write to a - # temporary directory outside of the FusionFS mounted directory then manually copy across + output_dir = output_dir / f"pcgr_{chunk_nbr}" if chunk_nbr is not None else output_dir - temp_dir = tempfile.TemporaryDirectory() - pcgr_output_dir = output_dir / 'pcgr/' + if output_dir.exists(): + logger.warning(f"Output directory '{output_dir}' already exists and will be overwritten") + shutil.rmtree(output_dir) + + # Create output directory + output_dir.mkdir(parents=True, exist_ok=True) if not sample_id: sample_id = 'nosampleset' @@ -124,19 +138,18 @@ def run_somatic(input_fp, pcgr_refdata_dir, output_dir, threads=1, pcgr_conda=No command_args = [ f'--sample_id {sample_id}', f'--input_vcf {input_fp}', + f'--vep_dir {vep_dir}', + f'--refdata_dir {pcgr_refdata_dir}', f'--tumor_dp_tag TUMOR_DP', f'--tumor_af_tag TUMOR_AF', f'--control_dp_tag NORMAL_DP', f'--control_af_tag NORMAL_AF', - f'--pcgr_dir {pcgr_refdata_dir}', f'--genome_assembly grch38', f'--assay WGS', - f'--estimate_signatures', - f'--estimate_msi_status', - f'--estimate_tmb', - f'--show_noncoding', - f'--vcfanno_n_proc {threads}', - f'--vep_pick_order biotype,rank,appris,tsl,ccds,canonical,length,mane', + *([] if disable_estimates else ['--estimate_signatures', '--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', ] # NOTE(SW): VEP pick order is applied as a successive filter: @@ -161,8 +174,11 @@ def run_somatic(input_fp, pcgr_refdata_dir, output_dir, threads=1, pcgr_conda=No if ploidy: command_args.append(f'--tumor_ploidy {ploidy}') + if chunk_nbr is not None: + command_args.append(f'--no_html') + # NOTE(SW): placed here to always have output directory last - command_args.append(f'--output_dir {temp_dir.name}') + command_args.append(f'--output_dir {output_dir}') delimiter_padding = ' ' * 10 delimiter = f' \\\n{delimiter_padding}' @@ -171,7 +187,7 @@ def run_somatic(input_fp, pcgr_refdata_dir, output_dir, threads=1, pcgr_conda=No command = fr''' pcgr \ - {command_args_str} + {command_args_str} ''' if pcgr_conda: @@ -179,44 +195,62 @@ def run_somatic(input_fp, pcgr_refdata_dir, output_dir, threads=1, pcgr_conda=No command_formatting = '\n' + ' ' * 4 command = command_formatting + command_conda + command - util.execute_command(command) + # Log file path + log_file_path = output_dir / "run_somatic.log" + + # Run the command and redirect output to the log file + util.execute_command(command, log_file_path=log_file_path) - shutil.copytree(temp_dir.name, pcgr_output_dir) + pcgr_tsv_fp = pathlib.Path(output_dir) / f'{sample_id}.pcgr.grch38.snv_indel_ann.tsv.gz' + pcgr_vcf_fp = pathlib.Path(output_dir) / f'{sample_id}.pcgr.grch38.pass.vcf.gz' - return pcgr_output_dir + # Check if both files exist + if not pcgr_tsv_fp.exists(): + raise FileNotFoundError(f"Expected file {pcgr_tsv_fp} not found.") + if not pcgr_vcf_fp.exists(): + raise FileNotFoundError(f"Expected file {pcgr_vcf_fp} not found.") + return pcgr_tsv_fp, pcgr_vcf_fp -def run_germline(input_fp, panel_fp, pcgr_refdata_dir, output_dir, threads=1, pcgr_conda=None, pcgrr_conda=None, sample_id=None): + +def run_germline(input_fp, panel_fp, pcgr_refdata_dir, vep_dir, output_dir, threads=1, pcgr_threads=4, pcgr_conda=None, pcgrr_conda=None, sample_id=None): + # threads: Nextflow process-level resource allocation (not wired to CPSR internals) + # pcgr_threads: CPSR-internal concurrency (vcfanno workers) if not sample_id: sample_id = 'nosampleset' - # NOTE(SW): Nextflow FusionFS v2.2.8 does not support PCGR output to S3; instead write to a - # temporary directory outside of the FusionFS mounted directory then manually copy across - temp_dir = tempfile.TemporaryDirectory() cpsr_output_dir = output_dir / 'cpsr/' + if cpsr_output_dir.exists(): + logger.warning(f"Output directory '{cpsr_output_dir}' already exists and will be overwritten") + shutil.rmtree(cpsr_output_dir) + + # Create output directory + cpsr_output_dir.mkdir(parents=True, exist_ok=True) + command_args = [ f'--sample_id {sample_id}', f'--input_vcf {input_fp}', f'--genome_assembly grch38', f'--custom_list {panel_fp}', + f'--vep_dir {vep_dir}', + f'--refdata_dir {pcgr_refdata_dir}', # NOTE(SW): probably useful to add versioning information here; weigh against maintainence # burden f'--custom_list_name umccr_germline_panel', f'--pop_gnomad global', f'--classify_all', - f'--pcgr_dir {pcgr_refdata_dir}', - f'--vcfanno_n_proc {threads}', - f'--vep_pick_order biotype,rank,appris,tsl,ccds,canonical,length,mane', + f'--vcfanno_n_proc {pcgr_threads}', + f'--vep_pick_order biotype,rank,appris,tsl,ccds,canonical,length,mane_plus_clinical,mane_select', ] if pcgrr_conda: command_args.append(f'--pcgrr_conda {pcgrr_conda}') # NOTE(SW): placed here to always have output directory last - command_args.append(f'--output_dir {temp_dir.name}') + command_args.append(f'--output_dir {cpsr_output_dir}') delimiter_padding = ' ' * 10 delimiter = f' \\\n{delimiter_padding}' @@ -235,25 +269,24 @@ def run_germline(input_fp, panel_fp, pcgr_refdata_dir, output_dir, threads=1, pc util.execute_command(command) - shutil.copytree(temp_dir.name, cpsr_output_dir) - return cpsr_output_dir -def transfer_annotations_somatic(input_fp, tumor_name, filter_name, pcgr_dir, output_dir): +def transfer_annotations_somatic(input_fp, tumor_name, pcgr_vcf_fp, pcgr_tsv_fp, output_dir): # Set destination INFO field names and source TSV fields info_field_map = { constants.VcfInfo.PCGR_MUTATION_HOTSPOT: 'MUTATION_HOTSPOT', - constants.VcfInfo.PCGR_CLINVAR_CLNSIG: 'CLINVAR_CLNSIG', + constants.VcfInfo.PCGR_CLINVAR_CLASSIFICATION: 'CLINVAR_CLASSIFICATION', constants.VcfInfo.PCGR_TCGA_PANCANCER_COUNT: 'TCGA_PANCANCER_COUNT', constants.VcfInfo.PCGR_CSQ: 'CSQ', } - pcgr_tsv_fp = pathlib.Path(pcgr_dir) / 'nosampleset.pcgr_acmg.grch38.snvs_indels.tiers.tsv' - pcgr_vcf_fp = pathlib.Path(pcgr_dir) / 'nosampleset.pcgr_acmg.grch38.vcf.gz' + # Respect paths provided; do not override + pcgr_tsv_fp = pathlib.Path(pcgr_tsv_fp) + pcgr_vcf_fp = pathlib.Path(pcgr_vcf_fp) # Enforce matching defined and source INFO annotations - check_annotation_headers(info_field_map, pcgr_vcf_fp) + util.check_annotation_headers(info_field_map, pcgr_vcf_fp) # Gather PCGR annotation data for records pcgr_data = collect_pcgr_annotation_data(pcgr_tsv_fp, pcgr_vcf_fp, info_field_map) @@ -261,13 +294,11 @@ def transfer_annotations_somatic(input_fp, tumor_name, filter_name, pcgr_dir, ou # Open filehandles, set required header entries input_fh = cyvcf2.VCF(input_fp) - util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_TIER) + util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_ACTIONABILITY_TIER) util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_CSQ) util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_MUTATION_HOTSPOT) - util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_CLINVAR_CLNSIG) - util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_COSMIC_COUNT) + util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_CLINVAR_CLASSIFICATION) util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_TCGA_PANCANCER_COUNT) - util.add_vcf_header_entry(input_fh, constants.VcfInfo.PCGR_ICGC_PCAWG_COUNT) output_fp = output_dir / f'{tumor_name}.annotations.vcf.gz' output_fh = cyvcf2.Writer(output_fp, input_fh, 'wz') @@ -277,29 +308,25 @@ def transfer_annotations_somatic(input_fp, tumor_name, filter_name, pcgr_dir, ou # Do not process chrM since *snvs_indels.tiers.tsv does not include these annotations if record.CHROM == 'chrM': continue - # Immediately print out variants that were not annotated - if filter_name in record.FILTERS: - output_fh.write_record(record) - continue # Annotate and write - record_ann = annotate_record(record, pcgr_data) + record_ann = annotate_record(record, pcgr_data, allow_missing=True) output_fh.write_record(record_ann) + output_fh.close() def transfer_annotations_germline(input_fp, normal_name, cpsr_dir, output_dir): # Set destination INFO field names and source TSV fields + # Note: Only include fields that exist in CPSR v2.2.1 output info_field_map = { - constants.VcfInfo.CPSR_FINAL_CLASSIFICATION: 'FINAL_CLASSIFICATION', - constants.VcfInfo.CPSR_PATHOGENICITY_SCORE: 'CPSR_PATHOGENICITY_SCORE', constants.VcfInfo.CPSR_CLINVAR_CLASSIFICATION: 'CLINVAR_CLASSIFICATION', constants.VcfInfo.CPSR_CSQ: 'CSQ', } - cpsr_tsv_fp = pathlib.Path(cpsr_dir) / f'{normal_name}.cpsr.grch38.snvs_indels.tiers.tsv' + cpsr_tsv_fp = pathlib.Path(cpsr_dir) / f'{normal_name}.cpsr.grch38.classification.tsv.gz' cpsr_vcf_fp = pathlib.Path(cpsr_dir) / f'{normal_name}.cpsr.grch38.vcf.gz' # Enforce matching defined and source INFO annotations - check_annotation_headers(info_field_map, cpsr_vcf_fp) + util.check_annotation_headers(info_field_map, cpsr_vcf_fp) # Gather CPSR annotation data for records cpsr_data = collect_cpsr_annotation_data(cpsr_tsv_fp, cpsr_vcf_fp, info_field_map) @@ -307,8 +334,6 @@ def transfer_annotations_germline(input_fp, normal_name, cpsr_dir, output_dir): # Open filehandles, set required header entries input_fh = cyvcf2.VCF(input_fp) - util.add_vcf_header_entry(input_fh, constants.VcfInfo.CPSR_FINAL_CLASSIFICATION) - util.add_vcf_header_entry(input_fh, constants.VcfInfo.CPSR_PATHOGENICITY_SCORE) util.add_vcf_header_entry(input_fh, constants.VcfInfo.CPSR_CLINVAR_CLASSIFICATION) util.add_vcf_header_entry(input_fh, constants.VcfInfo.CPSR_CSQ) @@ -324,55 +349,33 @@ def transfer_annotations_germline(input_fp, normal_name, cpsr_dir, output_dir): # NOTE(SW): allow missing CPSR annotations for input variants, CPSR seems to drop some record_ann = annotate_record(record, cpsr_data, allow_missing=True) output_fh.write_record(record_ann) - - -def check_annotation_headers(info_field_map, vcf_fp): - # Ensure header descriptions from source INFO annotations match those defined here for the - # output file; force manual inspection where they do not match - vcf_fh = cyvcf2.VCF(vcf_fp) - for header_dst, header_src in info_field_map.items(): - - # Skip header lines that do not have an equivalent entry in the PCGR/CPSR VCF - try: - header_src_entry = vcf_fh.get_header_type(header_src) - except KeyError: - continue - - header_dst_entry = util.get_vcf_header_entry(header_dst) - # Remove leading and trailing quotes from source - header_src_description_unquoted = header_src_entry['Description'].strip('"') - assert header_src_description_unquoted == header_dst_entry['Description'] + output_fh.close() def collect_pcgr_annotation_data(tsv_fp, vcf_fp, info_field_map): # Gather all annotations from TSV data_tsv = dict() - with open(tsv_fp, 'r') as tsv_fh: + # Read gz or plain text based on extension + open_fn = gzip.open if str(tsv_fp).endswith('.gz') else open + 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 - # Process PCGR_TIER - # TIER_1, TIER_2, TIER_3, TIER_4, NONCODING - record_ann[constants.VcfInfo.PCGR_TIER] = record['TIER'].replace(' ', '_') - - # Count COSMIC hits - if record['COSMIC_MUTATION_ID'] == 'NA': - cosmic_count = 0 + # 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: - cosmic_count = len(record['COSMIC_MUTATION_ID'].split('&')) - record_ann[constants.VcfInfo.PCGR_COSMIC_COUNT] = cosmic_count - - # Count ICGC-PCAWG hits by taking sum of affected donors where the annotation value has - # the following format: project_code|tumor_type|affected_donors|tested_donors|frequency - icgc_pcawg_count = 0 - if record['ICGC_PCAWG_OCCURRENCE'] != 'NA': - for pcawg_hit_data in record['ICGC_PCAWG_OCCURRENCE'].split(','): - pcawrg_hit_data_fields = pcawg_hit_data.split('|') - affected_donors = int(pcawrg_hit_data_fields[2]) - icgc_pcawg_count += affected_donors - assert icgc_pcawg_count > 0 - record_ann[constants.VcfInfo.PCGR_ICGC_PCAWG_COUNT] = icgc_pcawg_count + tier_val = 'N' + record_ann[constants.VcfInfo.PCGR_ACTIONABILITY_TIER] = tier_val # Store annotation data data_tsv[key] = record_ann @@ -388,7 +391,7 @@ def collect_cpsr_annotation_data(tsv_fp, vcf_fp, info_field_map): # Gather annotations from TSV data_tsv = dict() gdot_re = re.compile('^(?P[\dXYM]+):g\.(?P\d+)(?P[A-Z]+)>(?P[A-Z]+)$') - with open(tsv_fp, 'r') as tsv_fh: + with gzip.open(tsv_fp, 'rt') as tsv_fh: for record in csv.DictReader(tsv_fh, delimiter='\t'): # Decompose CPSR 'GENOMIC_CHANGE' field into CHROM, POS, REF, and ALT re_result = gdot_re.match(record['GENOMIC_CHANGE']) @@ -410,6 +413,24 @@ def collect_cpsr_annotation_data(tsv_fp, vcf_fp, info_field_map): # Compile annotations, prefering TSV source 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" + 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')}" + pos = int(match.group('pos')) + ref = match.group('ref') + alt = match.group('alt') + return chrom, pos, ref, alt + def get_annotations_vcf(vcf_fp, info_field_map): data_vcf = dict() @@ -429,10 +450,20 @@ def get_annotations_vcf(vcf_fp, info_field_map): def get_annotation_entry_tsv(record, info_field_map): - # Set lookup key; PCGR/CPSR strips leading 'chr' from contig names - chrom = f'chr{record["CHROM"]}' - pos = int(record['POS']) - key = (chrom, pos, record['REF'], record['ALT']) + # If GENOMIC_CHANGE is present, parse it; otherwise, fallback to CHROM/POS/REF/ALT fields + chrom = pos = ref = alt = None + if 'GENOMIC_CHANGE' in record and record['GENOMIC_CHANGE']: + chrom, pos, ref, alt = parse_genomic_change(record['GENOMIC_CHANGE']) + else: + chrom = record.get('CHROM') or record.get('Chromosome') + pos = int(record.get('POS') or record.get('Start_position')) + ref = record.get('REF') + alt = record.get('ALT') + # Ensure chrom has 'chr' prefix + if chrom and not str(chrom).startswith('chr'): + chrom = f'chr{chrom}' + + key = (chrom, pos, ref, alt) record_ann = dict() for info_dst, info_src in info_field_map.items(): @@ -481,3 +512,195 @@ def annotate_record(record, annotations, *, allow_missing=False): record.INFO[info_enum.value] = v 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: + raise ValueError("max_variants must be a positive integer.") + + output_dir = pathlib.Path(output_dir / "vcf_chunks") + output_dir.mkdir(parents=True, exist_ok=True) + chunk_files = [] + chunk_number = 1 + variant_count = 0 + input_vcf = pathlib.Path(input_vcf) + base_filename = input_vcf.stem + chunk_filename = output_dir / f"{base_filename}_chunk{chunk_number}.vcf.gz" + chunk_files.append(chunk_filename) + # Open the input VCF using cyvcf2 + vcf_in = cyvcf2.VCF(input_vcf) + # Create a new VCF file for the first chunk + vcf_out = cyvcf2.Writer(str(chunk_filename), vcf_in, 'wz') + last_position = None + for record in vcf_in: + current_position = record.POS + # Check if we need to start a new chunk + if variant_count >= max_variants and (last_position is None or current_position != last_position): + # Close the current chunk file and start a new one + vcf_out.close() + chunk_number += 1 + chunk_filename = output_dir / f"{base_filename}_chunk{chunk_number}.vcf.gz" + chunk_files.append(chunk_filename) + vcf_out = cyvcf2.Writer(str(chunk_filename), vcf_in, 'wz') + variant_count = 0 + # Write the record to the current chunk + vcf_out.write_record(record) + variant_count += 1 + last_position = current_position + # Close the last chunk file + vcf_out.close() + vcf_in.close() + + for chunk_fp in chunk_files: + util.execute_command(f'bcftools index --tbi {chunk_fp}') + + logger.info(f"VCF file split into {len(chunk_files)} chunks.") + return chunk_files + +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) + 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 + # nothing to merge, so use that chunk directly as the merged output (bolt #26) + merged_vcf = merged_vcf_path.parent / f'{merged_vcf_path.name}.vcf.gz' + shutil.copy(pcgr_vcf_files[0], merged_vcf) + util.execute_command(f'bcftools index -t {merged_vcf}') + else: + merged_vcf = util.merge_vcf_files(pcgr_vcf_files, merged_vcf_path) + + return merged_vcf, merged_tsv_fp + + +def get_variant_filter_data(variant): + attribute_names = ( + 'tier', + 'difficult', + 'giab_conf', + 'intergenic', + 'intronic', + 'downstream', + 'upstream', + 'impacts_other', + ) + + data = {e: None for e in attribute_names} + + + data['tier'] = variant.INFO.get('PCGR_ACTIONABILITY_TIER') + + + info_keys = [k for k, v in variant.INFO] + + data['difficult'] = any(e.startswith('DIFFICULT') for e in info_keys) + data['giab_conf'] = 'GIAB_CONF' in info_keys + + + # NOTE(SW): GIAB_CONF always overrides DIFFICULT tags + if data['giab_conf'] and data['difficult']: + data['difficult'] = False + + + for impact in get_impacts(variant.INFO['PCGR_CSQ']): + if impact == 'intergenic_variant': + data['intergenic'] = True + elif impact == 'intron_variant': + data['intronic'] = True + elif impact == 'downstream_gene_variant': + data['downstream'] = True + elif impact == 'upstream_gene_variant': + data['upstream'] = True + elif impact: + data['impacts_other'] = True + else: + assert False + + return data + + +def get_impacts(csq_str_full): + impacts = set() + for csq_str in csq_str_full.split(','): + csq_tokens = csq_str.split('|') + impact_str = csq_tokens[1] + impacts.update(impact_str.split('&')) + return impacts + + +def determine_filter(data): + + for impact, region in get_ordering(tiers=False): + + # NOTE(SW): this is less efficient than nested loops since the outer block is reevaluated + # within what would be the inner loop each cycle; taking this route for cleaner code + + impacts_higher = get_impacts_higher(impact) + impact_filter = bool(data[impact]) and not any(bool(data[e]) for e in impacts_higher) + + region_filter = False + if region == 'none': + region_filter = not (data['difficult'] or data['giab_conf']) + else: + region_filter = data[region] + + if impact_filter and region_filter: + return (impact, region) + + return False + + +def get_variant_repr(variant): + return (variant.CHROM, variant.POS, variant.REF, tuple(variant.ALT)) + + +@functools.cache +def get_ordering(tiers=True, impacts=True, regions=True): + categories = [ + constants.PCGR_TIERS_FILTERING if tiers else None, + constants.VEP_IMPACTS_FILTER if impacts else None, + constants.GENOMIC_REGIONS_FILTERING if regions else None, + ] + + # NOTE(SW): I'm not aware of any noncoding impacts for TIER_[1-4] other than TERT but keeping + # in to be overly cautious + ordering_iter = itertools.product(*(c for c in categories if c)) + + return tuple(ordering_iter) + + +@functools.cache +def get_impacts_higher(impact): + impact_index = constants.VEP_IMPACTS_FILTER.index(impact) + if impact_index + 1 < len(constants.VEP_IMPACTS_FILTER): + impacts_higher = constants.VEP_IMPACTS_FILTER[impact_index+1:len(constants.VEP_IMPACTS_FILTER)] + else: + impacts_higher = list() + return impacts_higher diff --git a/bolt/logging_config.py b/bolt/logging_config.py new file mode 100644 index 0000000..0ddda2b --- /dev/null +++ b/bolt/logging_config.py @@ -0,0 +1,34 @@ +import logging +import sys +import pathlib +from datetime import datetime + +class IgnoreTinfoFilter(logging.Filter): + def filter(self, record): + # Exclude messages that contain the unwanted text. + if "no version information available" in record.getMessage(): + return False + return True + +def setup_logging(output_dir, script_name): + # Create a timestamp for the log file + timestamp = datetime.now().strftime('%Y%m%d_%H%M%S') + log_filename = f"{script_name}_{timestamp}.log" + log_file = pathlib.Path(output_dir) / log_filename + + # Create individual handlers. + console_handler = logging.StreamHandler(sys.stdout) + file_handler = logging.FileHandler(log_file) + + # Instantiate and attach the filter to both handlers. + tinfo_filter = IgnoreTinfoFilter() + console_handler.addFilter(tinfo_filter) + file_handler.addFilter(tinfo_filter) + + logging.basicConfig( + level=logging.DEBUG, + format='%(asctime)s - %(name)s - %(levelname)s - %(message)s', + handlers=[file_handler, console_handler] + ) + 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 b7333d9..c3b07fe 100644 --- a/bolt/util.py +++ b/bolt/util.py @@ -1,11 +1,18 @@ +import gzip import pathlib +import select import subprocess import sys import textwrap +import logging +from types import SimpleNamespace +import cyvcf2 from .common import constants +# Set up logging +logger = logging.getLogger(__name__) # TODO(SW): create note that number this assumes location of `//` def get_project_root(): @@ -15,46 +22,73 @@ def get_project_root(): return project_root -def execute_command(command): - command_prepared = command_prepare(command) +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()) - print(command_prepared) - - process = subprocess.run( - command_prepared, + process = subprocess.Popen( + prepared_command, shell=True, executable='/bin/bash', - capture_output=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True, encoding='utf-8', + bufsize=1, # line buffered ) - if process.returncode != 0: - print(process) - print(process.stderr) - sys.exit(1) - - return process - - -def command_prepare(command): - return f'set -o pipefail; {textwrap.dedent(command)}' + stdout_lines = [] + stderr_lines = [] + stream_map = { + process.stdout: (stdout_lines, logger.info), + process.stderr: (stderr_lines, logger.warning), + } + open_streams = set(stream_map) + log_file = log_file_path.open('a', encoding='utf-8') if log_file_path else None + + try: + # select multiplexes stdout and stderr in a single thread, preserving arrival order + # and preventing pipe buffer deadlock without threading races on log_file writes + while open_streams: + readable, _, _ = select.select(open_streams, [], []) + for stream in readable: + line = stream.readline() + if line: + # Filter out bash libtinfo.so.6 warnings + if 'libtinfo.so.6: no version information available' not in line: + lines, log_fn = stream_map[stream] + log_fn(line.rstrip()) + lines.append(line) + if log_file: + log_file.write(line) + log_file.flush() + else: + open_streams.discard(stream) + finally: + process.wait() + if log_file: + log_file.close() + if process.returncode != 0: + logger.error("Command failed with return code %d: %s", process.returncode, command.strip()) + raise subprocess.CalledProcessError( + process.returncode, command, + output=''.join(stdout_lines), + stderr=''.join(stderr_lines), + ) -#def count_vcf_records(fp, exclude_args=None): -# args = list() -# if exclude_args: -# args.append(f'-e \'{exclude_args}\'') -# -# args_str = ' '.join(args) -# command = f'bcftools view -H {args_str} {fp} | wc -l' -# -# result = execute_command(command) -# return int(result.stdout) - + return SimpleNamespace( + stdout=''.join(stdout_lines), + stderr=''.join(stderr_lines), + returncode=process.returncode, + pid=process.pid, + command=command, + ) def count_vcf_records(fp): result = execute_command(f'bcftools view -H {fp} | wc -l') - return int(result.stdout) + return int(result.stdout.strip()) def add_vcf_header_entry(fh, anno_enum): @@ -94,8 +128,114 @@ def get_qualified_vcf_annotation(anno_enum): assert anno_enum in constants.VcfInfo or anno_enum in constants.VcfFormat return f'{anno_enum.namespace}/{anno_enum.value}' - -#def add_vcf_filter(record, filter_enum): -# existing_filters = [e for e in record.FILTERS if e != 'PASS'] -# assert filter_enum.value not in existing_filters -# return ';'.join([*existing_filters, filter_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) + logger.info(f"Merged TSV written to: {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. + """ + 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 + 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''' + {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', + f'-o {merged_vcf}', + f'{merged_unsorted_vcf}' + ] + sort_command_args_str = delimiter.join(sort_command_args) + sort_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") + + # Optionally, remove the unsorted merged VCF file + if merged_unsorted_vcf.exists(): + merged_unsorted_vcf.unlink() + + return merged_vcf + +def check_annotation_headers(info_field_map, vcf_fp): + # Ensure header descriptions from source INFO annotations match those defined here for the + # output file; force manual inspection where they do not match + vcf_fh = cyvcf2.VCF(vcf_fp) + for header_dst, header_src in info_field_map.items(): + # Skip header lines that do not have an equivalent entry in the VCF + try: + header_src_entry = vcf_fh.get_header_type(header_src) + except KeyError: + continue + + header_dst_entry = get_vcf_header_entry(header_dst) + # Remove leading and trailing quotes from source + header_src_description_unquoted = header_src_entry['Description'].strip('"') + try: + assert header_src_description_unquoted == header_dst_entry['Description'] + except AssertionError: + print(f'Header description mismatch for {header_dst.value}') + print(f' src: {header_src_description_unquoted}') + print(f' dst: {header_dst_entry["Description"]}') + sys.exit(1) diff --git a/bolt/workflows/other/cancer_report.py b/bolt/workflows/other/cancer_report.py index 983f596..ebd2f3a 100644 --- a/bolt/workflows/other/cancer_report.py +++ b/bolt/workflows/other/cancer_report.py @@ -5,6 +5,9 @@ from ... import util +from ...logging_config import setup_logging + + @click.command(name='cancer_report') @@ -28,6 +31,10 @@ @click.option('--purple_dir', required=True, type=click.Path(exists=True)) @click.option('--virusbreakend_dir', required=True, type=click.Path(exists=True)) +@click.option('--mutpat_dir', required=True, type=click.Path(exists=True)) + +@click.option('--hrdetect_file', required=True, type=click.Path(exists=True)) +@click.option('--chord_file', required=True, type=click.Path(exists=True)) @click.option('--dragen_hrd_fp', required=False, type=click.Path(exists=True)) @@ -44,6 +51,9 @@ def entry(ctx, **kwargs): output_dir = pathlib.Path(kwargs['output_dir']) output_dir.mkdir(mode=0o755, parents=True, exist_ok=True) + script_name = pathlib.Path(__file__).stem + setup_logging(output_dir, script_name) + # Normalise SAGE variants and remove duplicates that arise for MutationalPattern compatibility decomposed_snv_vcf = normalise_and_dedup_sage_variants( kwargs['smlv_somatic_vcf_fp'], @@ -59,7 +69,7 @@ def entry(ctx, **kwargs): ) # Set other required argument values - batch_name = f'{kwargs["subject_name"]}_{kwargs["tumor_name"]}' + batch_name = f"{kwargs['subject_name']}_{kwargs['tumor_name']}" output_table_dir = output_dir / 'cancer_report_tables' # Optional dragen hrd argument @@ -99,6 +109,10 @@ def entry(ctx, **kwargs): --key_genes {kwargs['cancer_genes_fp']} \ --oncokb_genes {kwargs['oncokb_genes_fp']} \ \ + --mutpat_dir {kwargs['mutpat_dir']} \ + --hrdetect_file {kwargs['hrdetect_file']} \ + --chord_file {kwargs['chord_file']} \ + \ --img_dir {output_image_dir}/ \ --result_outdir {output_table_dir}/ \ --out_file {output_dir}/{kwargs['tumor_name']}.cancer_report.html diff --git a/bolt/workflows/smlv_germline/report.py b/bolt/workflows/smlv_germline/report.py index 0742c19..224907b 100644 --- a/bolt/workflows/smlv_germline/report.py +++ b/bolt/workflows/smlv_germline/report.py @@ -1,13 +1,17 @@ import pathlib import yaml +from ...logging_config import setup_logging import click +import logging from ... import util from ...common import pcgr +logger = logging.getLogger(__name__) + @click.command(name='report') @click.pass_context @@ -22,6 +26,7 @@ @click.option('--germline_panel_list_fp', required=True, type=click.Path(exists=True)) @click.option('--pcgr_data_dir', required=True, type=click.Path(exists=True)) +@click.option('--vep_dir', required=True, type=click.Path(exists=True)) @click.option('--threads', required=True, type=int, default=1) @@ -35,6 +40,9 @@ def entry(ctx, **kwargs): output_dir = pathlib.Path(kwargs['output_dir']) output_dir.mkdir(mode=0o755, parents=True, exist_ok=True) + # Set up logging + script_name = pathlib.Path(__file__).stem + setup_logging(output_dir, script_name) # BCFtools stats run_bcftool_stats(kwargs['vcf_unfiltered_fp'], kwargs['normal_name'], output_dir) @@ -69,6 +77,7 @@ def entry(ctx, **kwargs): cpsr_prep_fp, kwargs['germline_panel_list_fp'], kwargs['pcgr_data_dir'], + kwargs['vep_dir'], output_dir, threads=kwargs['threads'], pcgr_conda=kwargs['pcgr_conda'], diff --git a/bolt/workflows/smlv_somatic/annotate.py b/bolt/workflows/smlv_somatic/annotate.py index f163d5f..01d28b2 100644 --- a/bolt/workflows/smlv_somatic/annotate.py +++ b/bolt/workflows/smlv_somatic/annotate.py @@ -1,14 +1,13 @@ +import logging import pathlib - - import click import cyvcf2 - from ... import util from ...common import constants from ...common import pcgr - +from ...logging_config import setup_logging +logger = logging.getLogger(__name__) @click.command(name='annotate') @click.pass_context @@ -23,6 +22,7 @@ @click.option('--pon_dir', required=True, type=click.Path(exists=True)) @click.option('--pcgr_data_dir', required=True, type=click.Path(exists=True)) +@click.option('--vep_dir', required=True, type=click.Path(exists=True)) @click.option('--pcgr_conda', required=False, type=str) @click.option('--pcgrr_conda', required=False, type=str) @@ -30,6 +30,7 @@ @click.option('--threads', required=False, default=4, type=int) @click.option('--output_dir', required=True, type=click.Path()) +@click.option('--pcgr_variant_chunk_size', required=False, type=int, help='Override maximum variants per PCGR chunk.') def entry(ctx, **kwargs): '''Annotate variants with information from several sources\f @@ -44,6 +45,9 @@ def entry(ctx, **kwargs): output_dir = pathlib.Path(kwargs['output_dir']) output_dir.mkdir(mode=0o755, parents=True, exist_ok=True) + script_name = pathlib.Path(__file__).stem + setup_logging(output_dir, script_name) + # Set all FILTER="." to FILTER="PASS" as required by PURPLE filter_pass_fp = set_filter_pass(kwargs['vcf_fp'], kwargs['tumor_name'], output_dir) @@ -73,55 +77,67 @@ def entry(ctx, **kwargs): ) # Annotate with cancer-related and functional information from a range of sources using PCGR - # - Select variants to process - there is an upper limit for PCGR of around 500k # - 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] # - Known mutation hotspot [INFO/PCGR_MUTATION_HOTSPOT] - # - ClinVar clinical significant [INFO/PCGR_CLINVAR_CLNSIG] - # - Hits in COSMIC [INFO/PCGR_COSMIC_COUNT] + # - ClinVar clinical significant [INFO/PCGR_CLNSIG] # - Hits in TCGA [INFO/PCGR_TCGA_PANCANCER_COUNT] # - Hits in PCAWG [INFO/PCGR_ICGC_PCAWG_COUNT] - # Set selected data or full input - selection_data = select_variants( - pon_fp, - kwargs['tumor_name'], - kwargs['cancer_genes_fp'], - output_dir, - ) - - if not (pcgr_prep_input_fp := selection_data.get('filtered')): - pcgr_prep_input_fp = selection_data['selected'] # Prepare VCF for PCGR annotation pcgr_prep_fp = pcgr.prepare_vcf_somatic( - pcgr_prep_input_fp, + pon_fp, kwargs['tumor_name'], kwargs['normal_name'], output_dir, ) - # Run PCGR - pcgr_dir = pcgr.run_somatic( - pcgr_prep_fp, - kwargs['pcgr_data_dir'], - output_dir, - threads=kwargs['threads'], - pcgr_conda=kwargs['pcgr_conda'], - pcgrr_conda=kwargs['pcgrr_conda'], - ) + 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}") + + # Run PCGR in chunks if exceeding the maximum allowed for somatic variants + chunk_size = kwargs.get('pcgr_variant_chunk_size') + if chunk_size is not None and chunk_size <= 0: + raise click.BadParameter('must be a positive integer', param_hint='--pcgr_variant_chunk_size') + chunk_size = chunk_size or constants.MAX_SOMATIC_VARIANTS + + if total_variants > chunk_size: + vcf_chunks = pcgr.split_vcf(pcgr_prep_fp, output_dir, max_variants=chunk_size) + pcgr_tsv_fp, pcgr_vcf_fp = pcgr.run_somatic_chunk( + vcf_chunks, + kwargs['pcgr_data_dir'], + kwargs['vep_dir'], + output_dir, + pcgr_output_dir, + kwargs['threads'], + kwargs['pcgr_conda'], + kwargs['pcgrr_conda'], + ) + else: + pcgr_tsv_fp, pcgr_vcf_fp = pcgr.run_somatic( + pcgr_prep_fp, + kwargs['pcgr_data_dir'], + kwargs['vep_dir'], + pcgr_output_dir, + chunk_nbr=None, + threads=kwargs['threads'], + pcgr_conda=kwargs['pcgr_conda'], + pcgrr_conda=kwargs['pcgrr_conda'], + ) # Transfer PCGR annotations to full set of variants pcgr.transfer_annotations_somatic( - selection_data['selected'], + pon_fp, kwargs['tumor_name'], - selection_data.get('filter_name'), - pcgr_dir, + pcgr_vcf_fp, + pcgr_tsv_fp, output_dir, ) - + logger.info("Annotation process completed") def set_filter_pass(input_fp, tumor_name, output_dir): output_fp = output_dir / f'{tumor_name}.set_filter_pass.vcf.gz' @@ -136,7 +152,6 @@ def set_filter_pass(input_fp, tumor_name, output_dir): return output_fp - def general_annotations(input_fp, tumor_name, threads, annotations_dir, output_dir): toml_fp = pathlib.Path(annotations_dir) / 'vcfanno_annotations.toml' @@ -169,7 +184,6 @@ def panel_of_normal_annotations(input_fp, tumor_name, threads, pon_dir, output_d util.execute_command(command) return output_fp - def select_variants(input_fp, tumor_name, cancer_genes_fp, output_dir): # Exclude variants until we hopefully move the needle below the threshold diff --git a/bolt/workflows/smlv_somatic/filter.py b/bolt/workflows/smlv_somatic/filter.py index e46f285..aa4f7c4 100644 --- a/bolt/workflows/smlv_somatic/filter.py +++ b/bolt/workflows/smlv_somatic/filter.py @@ -1,12 +1,14 @@ import click import pathlib - +import logging import cyvcf2 from ... import util from ...common import constants +from ...logging_config import setup_logging +logger = logging.getLogger(__name__) @click.command(name='filter') @@ -26,6 +28,9 @@ def entry(ctx, **kwargs): output_dir = pathlib.Path(kwargs['output_dir']) output_dir.mkdir(mode=0o755, parents=True, exist_ok=True) + script_name = pathlib.Path(__file__).stem + setup_logging(output_dir, script_name) + # Open input VCF and set required header entries for output in_fh = cyvcf2.VCF(kwargs['vcf_fp']) header_filters = ( @@ -37,7 +42,7 @@ def entry(ctx, **kwargs): constants.VcfFilter.ENCODE, constants.VcfFilter.GNOMAD_COMMON, constants.VcfInfo.SAGE_HOTSPOT_RESCUE, - constants.VcfInfo.PCGR_TIER_RESCUE, + constants.VcfInfo.PCGR_ACTIONABILITY_TIER_RESCUE, constants.VcfInfo.CLINICAL_POTENTIAL_RESCUE, constants.VcfInfo.RESCUED_FILTERS_EXISTING, constants.VcfInfo.RESCUED_FILTERS_PENDING, @@ -126,7 +131,11 @@ def set_filter_data(record, tumor_index): # PON filter ## # NOTE(SW): 'max' is inclusive - keeps variants with 0 to n-1 PON hits; preserved from Umccrise - pon_count = record.INFO.get(constants.VcfInfo.PON_COUNT.value, 0) + pon_count = get_record_value( + record, + constants.VcfInfo.PON_COUNT.value, + default=0 + ) if pon_count >= constants.PON_HIT_THRESHOLD: filters.append(constants.VcfFilter.PON) @@ -141,7 +150,14 @@ def set_filter_data(record, tumor_index): ## # NOTE(SW): rounding is essential here for accurate comparison; cyvcf2 floating-point error # means INFO/gnomAD_AF=0.01 can be represented as 0.009999999776482582 - gnomad_af = round(record.INFO.get(constants.VcfInfo.GNOMAD_AF.value, 0), 3) + gnomad_af = round( + get_record_value( + record, + constants.VcfInfo.GNOMAD_AF.value, + default=0.0 + ), + 3, + ) if gnomad_af >= constants.MAX_GNOMAD_AF: filters.append(constants.VcfFilter.GNOMAD_COMMON) @@ -162,15 +178,15 @@ def set_filter_data(record, tumor_index): ## # PCGR tier rescue ## - pcgr_tier = record.INFO.get(constants.VcfInfo.PCGR_TIER.value) - if pcgr_tier in constants.PCGR_TIERS_RESCUE: - info_rescue.append(constants.VcfInfo.PCGR_TIER_RESCUE) + pcgr_tier = record.INFO.get(constants.VcfInfo.PCGR_ACTIONABILITY_TIER.value) + if pcgr_tier in constants.PCGR_ACTIONABILITY_TIER_RESCUE: + info_rescue.append(constants.VcfInfo.PCGR_ACTIONABILITY_TIER_RESCUE) ## # SAGE hotspot rescue ## # NOTE(SW): effectively reverts any FILTERs that may have been applied above - if record.INFO.get(constants.VcfInfo.SAGE_HOTSPOT.value) is not None: + if get_record_value(record, constants.VcfInfo.SAGE_HOTSPOT.value) is not None: info_rescue.append(constants.VcfInfo.SAGE_HOTSPOT_RESCUE) ## @@ -181,19 +197,26 @@ def set_filter_data(record, tumor_index): # single CLINICAL_POTENTIAL_RESCUE flag # Get ClinVar clinical significance entries - clinvar_clinsig = record.INFO.get(constants.VcfInfo.PCGR_CLINVAR_CLNSIG.value, '') + clinvar_clinsig = get_record_value( + record, + constants.VcfInfo.PCGR_CLINVAR_CLASSIFICATION.value, + default='', + ) clinvar_clinsigs = clinvar_clinsig.split(',') # Hit counts in relevant reference somatic mutation databases - cosmic_count = record.INFO.get(constants.VcfInfo.PCGR_COSMIC_COUNT.value, 0) - tcga_pancancer_count = record.INFO.get(constants.VcfInfo.PCGR_TCGA_PANCANCER_COUNT.value, 0) - icgc_pcawg_count = record.INFO.get(constants.VcfInfo.PCGR_ICGC_PCAWG_COUNT.value, 0) + tcga_pancancer_count = get_record_value( + record, + constants.VcfInfo.PCGR_TCGA_PANCANCER_COUNT.value, + default=0 + ) + hmf_present = get_record_value(record, constants.VcfInfo.HMF_HOTSPOT.value) + pcgr_hotspot_present = get_record_value(record, constants.VcfInfo.PCGR_MUTATION_HOTSPOT.value) + if ( - record.INFO.get(constants.VcfInfo.HMF_HOTSPOT.value) is not None or - record.INFO.get(constants.VcfInfo.PCGR_MUTATION_HOTSPOT.value) is not None or + hmf_present is not None or + pcgr_hotspot_present is not None or any(e in clinvar_clinsigs for e in constants.CLINVAR_CLINSIGS_RESCUE) or - cosmic_count >= constants.MIN_COSMIC_COUNT_RESCUE or - tcga_pancancer_count >= constants.MIN_TCGA_PANCANCER_COUNT_RESCUE or - icgc_pcawg_count >= constants.MIN_ICGC_PCAWG_COUNT_RESCUE + tcga_pancancer_count >= constants.MIN_TCGA_PANCANCER_COUNT_RESCUE ): info_rescue.append(constants.VcfInfo.CLINICAL_POTENTIAL_RESCUE) @@ -228,3 +251,22 @@ def set_filter_data(record, tumor_index): filters_existing = [e for e in record.FILTERS if e != 'PASS'] assert all(e not in filters_existing for e in filters_value) record.FILTER = ';'.join([*filters_existing, *filters_value]) + + +def get_record_value(record, key, *, default=None): + ''' + Return a INFO value, replacing '.', '' and missing entries by default value. + + This prevents placeholder values emitted by pcgr from triggering rescues + or filters. For example, an INFO placeholder '.' for + PCGR_TCGA_PANCANCER_COUNT will resolve to the integer 0 when called with + `get_record_value(record, 'PCGR_TCGA_PANCANCER_COUNT', default=0)`. + ''' + value = record.INFO.get(key) + if value is None: + return default + if value == '' or value == '.': + return default + if not isinstance(value, (int, float, str)): + logger.error(f'record ID: {record.ID} INFO value for {key}: {value} ({type(value)}, not int, float, str)') + return value \ No newline at end of file diff --git a/bolt/workflows/smlv_somatic/report.py b/bolt/workflows/smlv_somatic/report.py index bdf908f..c37f8c1 100644 --- a/bolt/workflows/smlv_somatic/report.py +++ b/bolt/workflows/smlv_somatic/report.py @@ -1,3 +1,4 @@ +import collections import csv import json import pathlib @@ -6,11 +7,14 @@ import click import cyvcf2 import yaml +import logging from ... import util from ...common import constants from ...common import pcgr +from ...logging_config import setup_logging +logger = logging.getLogger(__name__) @click.command(name='report') @@ -27,6 +31,7 @@ @click.option('--pcgrr_conda', required=False, type=str) @click.option('--pcgr_data_dir', required=False, type=str) +@click.option('--vep_dir', required=True, type=click.Path(exists=True)) @click.option('--purple_purity_fp', required=True, type=click.Path(exists=True)) @click.option('--cancer_genes_fp', required=True, type=click.Path(exists=True)) @@ -45,6 +50,9 @@ def entry(ctx, **kwargs): output_dir = pathlib.Path(kwargs['output_dir']) output_dir.mkdir(mode=0o755, parents=True, exist_ok=True) + script_name = pathlib.Path(__file__).stem + setup_logging(output_dir, script_name) + # BCFtools stats bcftools_vcf_fp = bcftools_stats_prepare(kwargs['vcf_fp'], kwargs['tumor_name'], output_dir) run_bcftools_stats(bcftools_vcf_fp, kwargs['tumor_name'], output_dir) @@ -62,7 +70,7 @@ def entry(ctx, **kwargs): # Variant type counts # NOTE(SW): this is intended to preserve counts in the MultiQC report variant_counts_types_dragen = count_variant_types(kwargs['vcf_dragen_fp']) - variant_counts_types_bolt= count_variant_types(kwargs['vcf_fp']) + variant_counts_types_bolt = count_variant_types(kwargs['vcf_fp']) # NOTE(SW): using pass variants only for now @@ -107,17 +115,29 @@ 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_input_vcf_fp = kwargs['vcf_fp'] + else: + pcgr_input_vcf_fp = select_pcgr_variants( + kwargs['vcf_fp'], + kwargs['cancer_genes_fp'], + kwargs['tumor_name'], + output_dir, + ) + pcgr_prep_fp = pcgr.prepare_vcf_somatic( - kwargs['vcf_fp'], + 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'], - output_dir, + kwargs['vep_dir'], + pcgr_output_dir, threads=kwargs['threads'], pcgr_conda=kwargs['pcgr_conda'], pcgrr_conda=kwargs['pcgrr_conda'], @@ -143,7 +163,7 @@ def bcftools_stats_prepare(input_fp, tumor_name, output_dir): elif record.INFO.get('SAGE_NOVEL') is not None: record.QUAL = None else: - assert False + raise AssertionError(f'Record at {record.CHROM}:{record.POS} has neither SQ nor SAGE_NOVEL — cannot determine QUAL') output_fh.write_record(record) @@ -153,8 +173,8 @@ def bcftools_stats_prepare(input_fp, tumor_name, output_dir): def run_bcftools_stats(input_fp, tumor_name, output_dir): output_fp = output_dir / f'{tumor_name}.somatic.bcftools_stats.txt' command = fr''' - bcftools stats {input_fp} | \ - sed '6 s#{input_fp}$#{tumor_name}#' > {output_fp} + bcftools stats '{input_fp}' | \ + sed '6 s#{input_fp}$#{tumor_name}#' > '{output_fp}' ''' util.execute_command(command) @@ -281,9 +301,106 @@ def count_variant_process(vcf_fp): if not record.FILTER or rescued_filters: counts['filter_pass'] += 1 + counts['is_hypermutated'] = counts['dragen'] > constants.MAX_SOMATIC_VARIANTS return counts +def select_pcgr_variants(vcf_fp, cancer_genes_fp, tumor_name, output_dir): + """Filter variants for hypermutated samples to stay below PCGR's 500k limit. + + Retained variants (hotspot / panel) are skipped for tiered filtering. The remaining + variants are classified by (tier, impact, region) and categories are dropped in + priority order — NONCODING first, TIER_1 last — until the count falls within + MAX_SOMATIC_VARIANTS. Raises RuntimeError if the limit cannot be reached (i.e. + retained variants alone exceed it). + + Returns the path to the pass VCF (variants that survived filtering). A second + traceability VCF with FILTER=PCGR_count_limit set on dropped variants is written + alongside it. + """ + # Annotate variants in UMCCR somatic gene panel + fp_annotated_out = output_dir / f'{tumor_name}.umccr_panel_variants_annotated.vcf.gz' + util.execute_command(fr''' + bcftools annotate \ + --annotations <(awk 'BEGIN {{ OFS="\t" }} {{ print $1, $2-2000, $3+2000, "1" }}' {cancer_genes_fp}) \ + --header-line '{util.get_vcf_header_line(constants.VcfInfo.PANEL)}' \ + --columns CHROM,FROM,TO,{constants.VcfInfo.PANEL.value} \ + --output {fp_annotated_out} \ + {vcf_fp} + ''') + + # Set filter category for each variant + variants_sorted = collections.defaultdict(list) + variant_count = 0 + + 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): + continue + + data = pcgr.get_variant_filter_data(variant) + variant_filter = pcgr.determine_filter(data) + if not variant_filter: + raise AssertionError( + f'determine_filter returned no category for variant {variant_repr} (data={data})' + ) + + filter_category = (data['tier'], *variant_filter) + variants_sorted[filter_category].append(variant_repr) + + # Determine the filter categories needed to bring the count under MAX_SOMATIC_VARIANTS + filter_sum = 0 + filter_categories = list() + for key in pcgr.get_ordering(): + if (variant_count - filter_sum) <= constants.MAX_SOMATIC_VARIANTS: + break + filter_sum += len(variants_sorted.get(key, [])) + filter_categories.append(key) + + filter_variants = set() + for key in filter_categories: + filter_variants.update(variants_sorted[key]) + + expected_output = variant_count - len(filter_variants) + if expected_output > constants.MAX_SOMATIC_VARIANTS: + raise RuntimeError( + f'select_pcgr_variants failed to cap variants for {tumor_name}: ' + f'{expected_output} > {constants.MAX_SOMATIC_VARIANTS}' + ) + + logger.info('%s: select_pcgr_variants total=%d filtered=%d output=%d', + tumor_name, variant_count, len(filter_variants), expected_output) + + # Write passing variants; write all variants (with FILTER set) for traceability. + # Re-opening the original vcf_fp (not fp_annotated_out) so that the temporary + # PANEL annotation used for tiered filtering is not propagated to PCGR input. + fh_in = cyvcf2.VCF(vcf_fp) + util.add_vcf_header_entry(fh_in, constants.VcfFilter.PCGR_COUNT_LIMIT) + + # NOTE(SW): creating an additional VCF with all records for traceability + fp_out = output_dir / f'{tumor_name}.pcgr_hypermutated.pass.vcf.gz' + fp_set_out = output_dir / f'{tumor_name}.pcgr_hypermutated.filters_set.vcf.gz' + + fh_out = cyvcf2.Writer(fp_out, fh_in, 'wz') + fh_set_out = cyvcf2.Writer(fp_set_out, fh_in, 'wz') + + for variant in fh_in: + variant_repr = pcgr.get_variant_repr(variant) + if variant_repr not in filter_variants: + # Write only passing + fh_out.write_record(variant) + else: + variant.FILTER = constants.VcfFilter.PCGR_COUNT_LIMIT.value + # Write all variants including those with FILTER set + fh_set_out.write_record(variant) + + fh_out.close() + fh_set_out.close() + + return fp_out + + def parse_purple_purity_file(fp): with open(fp, 'r') as fh: entries = list(csv.DictReader(fh, delimiter='\t')) diff --git a/bolt/workflows/smlv_somatic/rescue.py b/bolt/workflows/smlv_somatic/rescue.py index c92df94..ff69e6e 100644 --- a/bolt/workflows/smlv_somatic/rescue.py +++ b/bolt/workflows/smlv_somatic/rescue.py @@ -7,10 +7,14 @@ import click import cyvcf2 +import logging from ... import util from ...common import constants +from ...logging_config import setup_logging + +logger = logging.getLogger(__name__) @click.command(name='rescue') @@ -37,6 +41,9 @@ def entry(ctx, **kwargs): output_dir = pathlib.Path(kwargs['output_dir']) output_dir.mkdir(mode=0o755, parents=True, exist_ok=True) + script_name = pathlib.Path(__file__).stem + setup_logging(output_dir, script_name) + # Select PASS SAGE variants in hotspots and then split into existing and novel calls sage_pass_vcf_fp = select_sage_pass_hotspot( kwargs['sage_vcf_fp'], @@ -106,20 +113,28 @@ def annotate_existing_sage_calls(input_fp, tumor_name, sage_vcf_fp, output_dir): # Get input file handle input_fh = cyvcf2.VCF(input_fp) + # Define expected SAGE annotations used for VCF header consistency check + info_field_map_sage = { + constants.VcfInfo.SAGE_HOTSPOT: 'SAGE_HOTSPOT', + constants.VcfInfo.SAGE_NOVEL: 'SAGE_NOVEL', + constants.VcfInfo.SAGE_RESCUE: 'SAGE_RESCUE', + constants.VcfFormat.SAGE_AD: 'SAGE_AD', + constants.VcfFormat.SAGE_AF: 'SAGE_AF', + constants.VcfFormat.SAGE_DP: 'SAGE_DP', + constants.VcfFormat.SAGE_SB: 'SAGE_SB', + constants.VcfFilter.SAGE_LOWCONF: 'SAGE_LOWCONF', + } + + util.check_annotation_headers(info_field_map_sage, input_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) - util.add_vcf_header_entry(input_fh, constants.VcfInfo.SAGE_HOTSPOT) util.add_vcf_header_entry(input_fh, constants.VcfInfo.SAGE_RESCUE) - - # TODO(SW): check that defined header descriptions match those in the SAGE fp; collect as list - # here and iterate to check and then add to input_fp header also in another loop - util.add_vcf_header_entry(input_fh, constants.VcfFormat.SAGE_AD) util.add_vcf_header_entry(input_fh, constants.VcfFormat.SAGE_AF) util.add_vcf_header_entry(input_fh, constants.VcfFormat.SAGE_DP) util.add_vcf_header_entry(input_fh, constants.VcfFormat.SAGE_SB) - # Open output file and use header from input file output_fp = output_dir / f'{tumor_name}.anno.vcf.gz' output_fh = cyvcf2.Writer(output_fp, input_fh, 'wz') diff --git a/bolt/workflows/sv_somatic/annotate.py b/bolt/workflows/sv_somatic/annotate.py index 2988e53..8577922 100644 --- a/bolt/workflows/sv_somatic/annotate.py +++ b/bolt/workflows/sv_somatic/annotate.py @@ -5,10 +5,13 @@ import click import cyvcf2 import pysam +import logging from ... import util +logger = logging.getLogger(__name__) + @click.command(name='annotate') @click.pass_context diff --git a/conda/env/bolt_env.yml b/conda/env/bolt_env.yml index e27bd24..fe7200b 100644 --- a/conda/env/bolt_env.yml +++ b/conda/env/bolt_env.yml @@ -15,4 +15,6 @@ dependencies: - pybedtools - python >=3.10 - pyyaml + - setuptools <81 - vcfanno ==0.3.5 + - ncurses>=6.3 diff --git a/docker/Dockerfile.gpgr b/docker/Dockerfile.gpgr index 841c9d0..d22d5ec 100644 --- a/docker/Dockerfile.gpgr +++ b/docker/Dockerfile.gpgr @@ -18,7 +18,7 @@ RUN \ RUN \ conda install --prefix /env/ \ - 'r-gpgr ==2.2.12' \ + 'r-gpgr ==2.3.1' \ 'r-sigrap ==0.1.1' \ 'bioconductor-bsgenome.hsapiens.ucsc.hg38 ==1.4.5' \ 'bioconductor-txdb.hsapiens.ucsc.hg38.knowngene ==3.16.0' \ diff --git a/docker/Dockerfile.pcgr b/docker/Dockerfile.pcgr index 4ea4932..c63e694 100644 --- a/docker/Dockerfile.pcgr +++ b/docker/Dockerfile.pcgr @@ -15,13 +15,13 @@ RUN \ conda create \ --solver libmamba \ --name pcgr \ - --file https://raw.githubusercontent.com/sigven/pcgr/v1.4.1/conda/env/lock/pcgr-linux-64.lock + --file https://raw.githubusercontent.com/sigven/pcgr/refs/tags/v2.2.5/conda/env/lock/pcgr-linux-64.lock RUN \ conda create \ --solver libmamba \ --name pcgrr \ - --file https://raw.githubusercontent.com/sigven/pcgr/v1.4.1/conda/env/lock/pcgrr-linux-64.lock + --file https://raw.githubusercontent.com/sigven/pcgr/refs/tags/v2.2.5/conda/env/lock/pcgrr-linux-64.lock COPY ./conda/env/bolt_env.yml /tmp/ RUN \ @@ -44,5 +44,5 @@ FROM quay.io/bioconda/base-glibc-busybox-bash:2.1.0 # Copy Conda install and all environments COPY --from=build /opt/conda/ /opt/conda/ -ENV PATH="/opt/conda/envs/bolt/bin:/opt/conda/bin:${PATH}" +ENV PATH="/opt/conda/envs/bolt/bin:/opt/conda/envs/pcgr/bin:/opt/conda/bin:${PATH}" ENV LD_LIBRARY_PATH="/opt/conda/lib/" diff --git a/docs/adr/001-max-somatic-variants-450k.md b/docs/adr/001-max-somatic-variants-450k.md new file mode 100644 index 0000000..70695a2 --- /dev/null +++ b/docs/adr/001-max-somatic-variants-450k.md @@ -0,0 +1,58 @@ +# ADR-001: Cap somatic variants at 450k before PCGR (not 500k) + +**Date:** 2026-07-22 +**Status:** Accepted +**Context:** bolt #35, sash #52 +**Deciders:** Team (oral decision) + +## Context + +PCGR has a 500,000 variant threshold that triggers two problematic behaviours: + +1. **Python side** (`pcgr/main.py`, `pcgr_vars.MAX_VARIANTS_FOR_REPORT = 500_000`): + When input variants exceed 500k, PCGR silently drops intergenic, intronic, upstream_gene, and downstream_gene variants. This makes the TMB calculation an underestimate and removes our control over which variants appear in the report. + +2. **R side** (`pcgrr/R/main.R` ~line 954): + If variants are still ≥ 500k after the Python-side filtering, the HTML report is **silently not generated** — no error, just no output file. + +Neither triggers a hard failure. Both are silent. This makes the 500k boundary dangerous: a sample could pass through bolt, enter PCGR, and produce either a misleading report (missing variants, wrong TMB) or no report at all — with no error in the logs. + +## Decision + +Set `MAX_SOMATIC_VARIANTS = 450_000` in bolt. This is the threshold used by `select_pcgr_variants()` for tiered filtering and by `count_variant_process()` for the `is_hypermutated` flag. + +### Why 450k and not 500k? + +1. **Bolt's filter is coarse-grained.** It drops entire variant categories at once (e.g. all NONCODING_INTERGENIC variants). Output is always ≤ `MAX_SOMATIC_VARIANTS`, but it can undershoot significantly. + +2. **pcgrr's boundary check is `< 500000` (strict less-than).** If bolt outputs exactly 500,000 variants, pcgrr would skip HTML report generation. The 50k margin ensures we never land on this boundary. + +3. **PCGR's Python filter uses `> 500000`.** So even 500,001 variants trigger silent intergenic/intronic removal. The margin gives a buffer against any minor variant-count inflation during PCGR's annotation pipeline. + +### Why not lower (e.g. 400k)? + +Lower thresholds would cause more samples to be flagged as hypermutated and have their variants filtered, potentially losing clinically relevant variants in TIER_1/2 categories for samples that PCGR could actually handle fine. + +## Consequences + +- Samples with > 450k PASS somatic variants get tiered filtering (lowest-priority categories dropped first until under 450k) +- If tiered filtering cannot bring the count below 450k (because retained variants — hotspots, panel genes — alone exceed it), PCGR is skipped entirely for that sample +- Real case: L2100242 had 595k PASS variants; retained variants alone exceeded 450k → PCGR skipped (sash #52) +- The PCGR HTML report, MAF output, and VCF2MAF are lost for skipped samples; cancer report (gpgr) is unaffected + +## Alternatives Considered + +| Option | Pros | Cons | +|--------|------|------| +| 500k (match PCGR exactly) | Fewer samples filtered | Risk of silent report loss at boundary; TMB underestimates | +| 450k (chosen) | Safe margin; bolt controls filtering priority | Slightly more aggressive filtering on borderline samples | +| 400k | Extra safety | Unnecessary — 50k margin already covers all known edge cases | +| No limit (let PCGR handle it) | Simplest | PCGR's filtering is indiscriminate (drops all intergenic regardless of tier); TMB affected; report may vanish | + +## References + +- PCGR source: `sigven/pcgr`, `pcgr/main.py` ~line 559, `pcgrr/R/main.R` ~line 954 +- `pcgr_vars.MAX_VARIANTS_FOR_REPORT = 500_000` +- sash #52: hypermutated sample handling +- bolt #26: single-chunk PCGR merge guard +- Real failure: L2100242 (595k PASS variants, PANEL+hotspot retained variants alone > 450k) diff --git a/pyproject.toml b/pyproject.toml index 99b7714..6dd2fba 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,7 +8,7 @@ include = ["bolt*"] [project] name = "bolt" -version = "0.2.18" +version = "0.3.1" authors = [ {name = "Stephen Watts", email = "stephen.watts@umccr.org"}, ] diff --git a/tests/test_pcgr_hypermutated.py b/tests/test_pcgr_hypermutated.py new file mode 100644 index 0000000..b195668 --- /dev/null +++ b/tests/test_pcgr_hypermutated.py @@ -0,0 +1,636 @@ +"""Tests for hypermutated sample handling — tier ordering fix and variant trimming.""" +import gzip +import pathlib +import shutil +import tempfile +import unittest +from unittest.mock import patch + +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 + + +# 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] + + +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 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 _count_vcf(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_exactly_at_limit_nothing_dropped(self): + """When variant count == MAX_SOMATIC_VARIANTS, no filtering occurs. + + This documents the boundary: the check is `<=` so exactly-at-limit passes through. + Ensures bolt doesn't accidentally trigger PCGR's own 500k internal filter when + MAX_SOMATIC_VARIANTS < 500k. + """ + with tempfile.TemporaryDirectory() as tmp: + # 5 variants, limit=5 → no filtering needed + v = [] + for i in range(1, 4): # 3 TIER_1 + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(4, 6): # 2 NONCODING + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + count = self._run(v, limit=5, tmp=tmp) + self.assertEqual(count, 5) # All survive — exactly at limit + + def test_one_over_limit_drops_lowest_priority_category(self): + """When variant count is limit+1, the lowest-priority category is dropped entirely. + + This tests the coarse-grained nature of the filter: we drop whole categories, + so output may undershoot the limit significantly. This is by design — it keeps + the logic simple and deterministic, and the 450k→500k margin ensures PCGR's + own internal filter (which drops intergenic/intronic indiscriminately) never fires. + """ + with tempfile.TemporaryDirectory() as tmp: + # 6 variants (limit=5): 3 TIER_1 + 3 NONCODING + # One over limit → all 3 NONCODING dropped → output = 3 (undershoots limit) + v = [] + for i in range(1, 4): # 3 TIER_1 + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=1;PCGR_CSQ={_csq("intron_variant")}')) + for i in range(4, 7): # 3 NONCODING + v.append((i*10, f'PCGR_ACTIONABILITY_TIER=N;PCGR_CSQ={_csq("intergenic_variant")}')) + count = self._run(v, limit=5, tmp=tmp) + # All 3 NONCODING dropped (whole category), only 3 TIER_1 remain + self.assertEqual(count, 3) + self.assertLessEqual(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_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 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') + + +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 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) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_smlv_somatic_filter.py b/tests/test_smlv_somatic_filter.py index 41a0e96..f872227 100644 --- a/tests/test_smlv_somatic_filter.py +++ b/tests/test_smlv_somatic_filter.py @@ -46,10 +46,10 @@ def get_record_from_str(variant_str): def get_record( chrom='chr1', - pos='.', + pos='1', vid='.', - ref='.', - alt='.', + ref='A', + alt='T', qual='.', vfilter='.', info_data=None, @@ -160,15 +160,15 @@ def test_common_population_filter(self): def test_pcgr_tier_rescue(self): pcgr_tiers = [ - 'TIER_1', - 'TIER_2', + '1', + '2', ] - rescue_tag_str = bolt_constants.VcfInfo.PCGR_TIER_RESCUE.value + rescue_tag_str = bolt_constants.VcfInfo.PCGR_ACTIONABILITY_TIER_RESCUE.value for pcgr_tier in pcgr_tiers: record = get_record( **self.records['filter_min_af9.9'], - info_data={'PCGR_TIER': pcgr_tier}, + info_data={'PCGR_ACTIONABILITY_TIER': pcgr_tier}, ) smlv_somatic_filter.set_filter_data(record, 0) assert not record.FILTER @@ -189,9 +189,7 @@ def test_clinical_potential_rescue_general(self): info_data_sets = [ {'HMF_HOTSPOT': ''}, {'PCGR_MUTATION_HOTSPOT': ''}, - {'PCGR_COSMIC_COUNT': 11}, {'PCGR_TCGA_PANCANCER_COUNT': 6}, - {'PCGR_ICGC_PCAWG_COUNT': 6}, ] rescue_tag_str = bolt_constants.VcfInfo.CLINICAL_POTENTIAL_RESCUE.value @@ -217,7 +215,7 @@ def test_clinical_potential_rescue_clinvar_clinsig(self): for clinsig in clinsigs: record = get_record( **self.records['filter_min_af9.9'], - info_data={'PCGR_CLINVAR_CLNSIG': clinsig}, + info_data={'PCGR_CLINVAR_CLASSIFICATION': clinsig}, ) smlv_somatic_filter.set_filter_data(record, 0) assert not record.FILTER