diff --git a/CITATION.cff b/CITATION.cff index b303859..55d2fd2 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,8 +1,8 @@ cff-version: 1.2.0 message: "If you use this software, please cite it as below." title: "vartracker" -version: "2.2.0" -date-released: "2026-05-05" +version: "2.2.1" +date-released: "2026-05-06" license: MIT repository-code: "https://github.com/charlesfoster/vartracker" url: "https://github.com/charlesfoster/vartracker" @@ -22,7 +22,7 @@ preferred-citation: - family-names: Foster given-names: Charles title: "vartracker" - version: "2.2.0" + version: "2.2.1" doi: "10.5281/zenodo.18452274" url: "https://github.com/charlesfoster/vartracker" - date-released: "2026-05-05" + date-released: "2026-05-06" diff --git a/README.md b/README.md index b92c28b..30a3e45 100755 --- a/README.md +++ b/README.md @@ -167,7 +167,7 @@ Docker is a self-contained reproducible option. If you publish the image, record set it when running to include it in the run manifest: ```bash -export VARTRACKER_CONTAINER_IMAGE=ghcr.io/your-org/vartracker:2.0.0 +export VARTRACKER_CONTAINER_IMAGE=ghcr.io/your-org/vartracker:2.2.1 export VARTRACKER_CONTAINER_DIGEST=sha256:... ``` @@ -253,8 +253,10 @@ LoFreq primer-overlap rescue: - `bam` and `end-to-end` therefore run LoFreq with `--no-default-filter`, then apply the normal `lofreq filter` step so standard LoFreq PASS calls are unchanged. - With the default `--lofreq-primer-rescue auto`, the rescue step runs only when `--primer-bed` is supplied. In other words, `auto` means "use primer rescue when an amplicon primer scheme has been explicitly provided." - In `end-to-end` mode, the same `--primer-bed` is used for `samtools ampliconclip` and for rescue. In `bam` mode, vartracker does not clip the input BAMs; the primer BED is used only to identify primer-overlap sites for rescue. -- Rescue candidates must be single-ALT SNPs that overlap a primer interval, fail the default LoFreq filter, and pass conservative near-fixed thresholds (`AF>=0.95`, `DP>=100`, `DP4 alt count>=95`, `QUAL>=100`, `DP4 ref count<=20`) with one-sided alternate-strand support. Indels, multi-ALT records, lower-frequency variants, and non-primer-overlap variants are not rescued by this rule. +- Rescue candidates must be single-ALT SNPs that overlap a primer interval, fail LoFreq's default strand-bias filtering, and pass conservative near-fixed thresholds (`AF>=0.95`, `DP>=100`, `DP4 alt count>=95`, `QUAL>=100`, `DP4 ref count<=20`, minor ALT strand fraction `<=0.05`). Indels, multi-ALT records, lower-frequency variants, non-primer-overlap variants, and variants filtered for non-strand-bias reasons are not rescued by this rule. +- The raw LoFreq calls are retained as `_variants.raw.vcf.gz` and listed in the updated spreadsheet as `raw_vcf`. - Rescued variants are marked with `FILTER=RESCUED_PRIMER_OVERLAP`, `INFO/PRIMER_OVERLAP`, and `INFO/RESCUED_BY=overlap_primer_interval`; per-sample details are written to `_variants.rescued.tsv` and listed in the updated spreadsheet as `lofreq_rescued_tsv`. +- Variants called by raw LoFreq but filtered out of the final VCF are written to `_variants.filtered_out.tsv` with the LoFreq filter reason and core metrics. This is useful for auditing high-frequency calls that fail strand-bias or other LoFreq filters. - Use `--lofreq-primer-rescue off` to disable rescue even when a primer BED is supplied, or `--lofreq-primer-rescue on` to require rescue and fail if `--primer-bed` is missing. The rescue thresholds can be adjusted with the `--lofreq-rescue-*` options. Example amplicon run with primer rescue: @@ -284,10 +286,11 @@ Mode-specific expectations: - **End-to-end mode** requires `reads1` (and optionally `reads2`); remaining fields are generated. The `bam` and `end-to-end` workflows also write two consensus FASTA columns to -the updated Snakemake spreadsheet, plus the LoFreq rescue audit column: +the updated Snakemake spreadsheet, plus LoFreq audit columns: `consensus` for a simple consensus, `iupac_consensus` for an IUPAC-aware -consensus, and `lofreq_rescued_tsv` for the per-sample primer-overlap rescue -table. SNPs below +consensus, `raw_vcf` for raw LoFreq calls, `lofreq_rescued_tsv` for the +per-sample primer-overlap rescue table, and `lofreq_filtered_out_tsv` for raw +LoFreq records excluded from the final VCF. SNPs below `--consensus-snp-min-af` are ignored, SNPs from `--consensus-snp-min-af` up to `--consensus-snp-thresh` stay as reference bases in the simple consensus and become REF+ALT ambiguity codes in the IUPAC consensus, and SNPs at or above @@ -509,7 +512,9 @@ vartracker produces several output files: - **results.csv**: Comprehensive variant analysis with all metrics - **results_metadata.json**: Output schema version and results metadata +- **`_variants.raw.vcf.gz`** (`bam`/`end-to-end`): Raw LoFreq calls before default filtering and primer-overlap rescue - **`_variants.rescued.tsv`** (`bam`/`end-to-end`): LoFreq primer-overlap rescue audit table, empty when rescue is disabled or no variants are rescued +- **`_variants.filtered_out.tsv`** (`bam`/`end-to-end`): Raw LoFreq calls excluded from the final VCF, including filter reason and call metrics - **new_mutations.csv**: Mutations not present in the first sample - **persistent_new_mutations.csv**: New mutations that persist to the final sample - **cumulative_mutations.pdf**: Plot showing mutation accumulation over time @@ -566,7 +571,7 @@ The pipeline performs the following analysis: When using vartracker, please cite the software release you used. Citation metadata is provided in `CITATION.cff`, and GitHub releases are archived on Zenodo. -- Foster, C. (2026). *vartracker* (Version 2.2.0). Zenodo. https://doi.org/10.5281/zenodo.18452274 +- Foster, C. (2026). *vartracker* (Version 2.2.1). Zenodo. https://doi.org/10.5281/zenodo.18452274 Note: the DOI above is the Zenodo concept DOI for all versions; a version-specific DOI is minted by Zenodo after each GitHub release. diff --git a/pyproject.toml b/pyproject.toml index 1b2c97a..790e3cc 100755 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "vartracker" -version = "2.2.0" +version = "2.2.1" authors = [ {name = "Dr Charles Foster"}, ] diff --git a/tests/test_analysis.py b/tests/test_analysis.py index d7b886e..a97c007 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -133,18 +133,75 @@ def test_prepare_variant_heatmap_matrix_orders_variants_by_genome(): expected_long_label = f"nsp2:{head}+{middle}{tail}\n(G1946GT)" expected_index = [ - "nsp2:T809=\n(C1059T)", + "nsp2:T629=\n(C1059T)", expected_long_label, "S:D215G\n(A22206G)", ] assert list(matrix.index) == expected_index - assert matrix.loc["nsp2:T809=\n(C1059T)", "P0"] == 0.0 - assert matrix.loc["nsp2:T809=\n(C1059T)", "P1"] == 0.5 + assert matrix.loc["nsp2:T629=\n(C1059T)", "P0"] == 0.0 + assert matrix.loc["nsp2:T629=\n(C1059T)", "P1"] == 0.5 assert matrix.loc[expected_long_label, "P1"] == 0.8 assert matrix.loc["S:D215G\n(A22206G)", "P1"] == 1.0 +def test_prepare_variant_heatmap_matrix_normalises_starred_synonymous_label(): + table = pd.DataFrame( + [ + { + "gene": "ORF1ab", + "amino_acid_consequence": "924F", + "nsp_aa_change": "", + "type_of_change": "*synonymous", + "type_of_variant": "snp", + "alt_freq": "0.5 / 0.0", + "samples": "P0 / P1", + "variant": "C3037T", + "start": 3037, + }, + { + "gene": "ORF1ab", + "amino_acid_consequence": "924F", + "nsp_aa_change": "", + "type_of_change": "synonymous", + "type_of_variant": "snp", + "alt_freq": "0.0 / 0.6", + "samples": "P0 / P1", + "variant": "C3037T", + "start": 3037, + }, + ] + ) + + matrix = _prepare_variant_heatmap_matrix(table, ["P0", "P1"], 0.2, 0.3) + + assert list(matrix.index) == ["nsp3_PLpro:F106=\n(C3037T)"] + assert matrix.loc["nsp3_PLpro:F106=\n(C3037T)", "P0"] == 0.5 + assert matrix.loc["nsp3_PLpro:F106=\n(C3037T)", "P1"] == 0.6 + + +def test_prepare_variant_heatmap_matrix_repairs_stale_stop_gained_nsp_label(): + table = pd.DataFrame( + [ + { + "gene": "ORF1ab", + "amino_acid_consequence": "L889*", + "nsp_aa_change": "nsp3_PLpro:71L", + "type_of_change": "stop_gained", + "type_of_variant": "snp", + "alt_freq": "0.052", + "samples": "P0", + "variant": "T2931A", + "start": 2931, + }, + ] + ) + + matrix = _prepare_variant_heatmap_matrix(table, ["P0"], 0.0, 0.0) + + assert list(matrix.index) == ["nsp3_PLpro:L71*\n(T2931A)"] + + def test_prepare_variant_heatmap_matrix_excludes_selected_consequence_types(): table = pd.DataFrame( [ @@ -329,7 +386,7 @@ def test_process_joint_variants_only_adds_single_joint_prefix(tmp_path): "aa1_weight": "", "aa2_weight": "", "weight_difference": "", - "type_of_change": "joint_joint_frameshift", + "type_of_change": "joint_*frameshift", }, { "start": 101, diff --git a/tests/test_analysis_launcher.py b/tests/test_analysis_launcher.py index 977f832..9bdfb6c 100644 --- a/tests/test_analysis_launcher.py +++ b/tests/test_analysis_launcher.py @@ -50,6 +50,7 @@ def test_validate_lofreq_primer_rescue_on_requires_primer_bed(): min_alt_count=95, min_qual=100, max_ref_count=20, + max_minor_alt_fraction=0.05, ) @@ -63,4 +64,17 @@ def test_validate_lofreq_primer_rescue_rejects_invalid_threshold(): min_alt_count=95, min_qual=100, max_ref_count=20, + max_minor_alt_fraction=0.05, + ) + + with pytest.raises(ValueError, match="between 0 and 1"): + _validate_lofreq_primer_rescue( + "auto", + None, + min_af=0.95, + min_dp=100, + min_alt_count=95, + min_qual=100, + max_ref_count=20, + max_minor_alt_fraction=1.5, ) diff --git a/tests/test_constants.py b/tests/test_constants.py new file mode 100644 index 0000000..f6866a0 --- /dev/null +++ b/tests/test_constants.py @@ -0,0 +1,12 @@ +"""Tests for reference-coordinate formatting helpers.""" + +from __future__ import annotations + +from vartracker.constants import reformat_csq_notation + + +def test_reformat_csq_notation_preserves_orf1ab_stop_gained_nsp_change(): + reformatted, nsp_change = reformat_csq_notation("ORF1ab", "889L>889*") + + assert reformatted == "L889*" + assert nsp_change == "nsp3_PLpro:L71*" diff --git a/tests/test_lofreq_primer_rescue.py b/tests/test_lofreq_primer_rescue.py index 6822785..17476ed 100644 --- a/tests/test_lofreq_primer_rescue.py +++ b/tests/test_lofreq_primer_rescue.py @@ -1,7 +1,10 @@ import subprocess from pathlib import Path -from vartracker.lofreq_primer_rescue import rescue_lofreq_primer_variants +from vartracker.lofreq_primer_rescue import ( + lofreq_filter_with_audit, + rescue_lofreq_primer_variants, +) def test_rescue_lofreq_primer_variants_keeps_pass_and_rescues_overlap( @@ -15,7 +18,10 @@ def test_rescue_lofreq_primer_variants_keeps_pass_and_rescues_overlap( '##INFO=\n' "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" "chr1\t10\t.\tA\tG\t200\t.\tAF=0.98;DP=120;DP4=1,0,118,0\n" - "chr1\t30\t.\tC\tT\t150\t.\tAF=0.40;DP=120;DP4=30,30,30,30\n", + "chr1\t15\t.\tT\tC\t210\t.\tAF=0.98;DP=125;DP4=0,5,114,5\n" + "chr1\t18\t.\tG\tT\t220\t.\tAF=0.98;DP=125;DP4=0,4,115,5\n" + "chr1\t30\t.\tC\tT\t150\t.\tAF=0.40;DP=120;DP4=30,30,30,30\n" + "chr1\t50\t.\tT\tC\t160\t.\tAF=0.96;DP=130;DP4=2,2,60,65\n", encoding="utf-8", ) primers = tmp_path / "primers.bed" @@ -27,11 +33,14 @@ def test_rescue_lofreq_primer_variants_keeps_pass_and_rescues_overlap( '##INFO=\n' "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" "chr1\t10\t.\tA\tG\t200\tstrandbias\tAF=0.98;DP=120;DP4=1,0,118,0\n" + "chr1\t15\t.\tT\tC\t210\tsb_fdr\tAF=0.98;DP=125;DP4=0,5,114,5\n" + "chr1\t18\t.\tG\tT\t220\tlowqual\tAF=0.98;DP=125;DP4=0,4,115,5\n" "chr1\t30\t.\tC\tT\t150\tPASS\tAF=0.40;DP=120;DP4=30,30,30,30\n" + "chr1\t50\t.\tT\tC\t160\tsb_fdr\tAF=0.96;DP=130;DP4=2,2,60,65\n" ) def fake_run(cmd, check): - assert cmd[:4] == ["lofreq", "filter", "-i", str(raw_vcf)] + assert cmd[:5] == ["lofreq", "filter", "--print-all", "-i", str(raw_vcf)] assert check is True output_path = Path(cmd[cmd.index("-o") + 1]) output_path.write_text(default_filtered_vcf, encoding="utf-8") @@ -41,21 +50,26 @@ def fake_run(cmd, check): output_vcf = tmp_path / "final.vcf" rescued_tsv = tmp_path / "rescued.tsv" + filtered_out_tsv = tmp_path / "filtered_out.tsv" result = rescue_lofreq_primer_variants( raw_vcf=raw_vcf, primers_bed=primers, output_vcf=output_vcf, rescued_tsv=rescued_tsv, + filtered_out_tsv=filtered_out_tsv, ) assert result.normal_passed == 1 - assert result.rescued == 1 - assert result.discarded == 0 + assert result.rescued == 2 + assert result.discarded == 2 + assert result.filtered_out == 2 output_text = output_vcf.read_text(encoding="utf-8") assert "##FILTER=\n' + '##INFO=\n' + '##INFO=\n' + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + "chr1\t10\t.\tA\tG\t200\t.\tAF=0.98;DP=120;DP4=1,0,118,0\n" + "chr1\t50\t.\tT\tC\t160\t.\tAF=0.96;DP=130;DP4=2,2,60,65\n", + encoding="utf-8", + ) + default_filtered_vcf = ( + "##fileformat=VCFv4.2\n" + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + "chr1\t10\t.\tA\tG\t200\tPASS\tAF=0.98;DP=120;DP4=1,0,118,0\n" + "chr1\t50\t.\tT\tC\t160\tsb_fdr\tAF=0.96;DP=130;DP4=2,2,60,65\n" + ) + + def fake_run(cmd, check): + assert cmd[:5] == ["lofreq", "filter", "--print-all", "-i", str(raw_vcf)] + assert check is True + output_path = Path(cmd[cmd.index("-o") + 1]) + output_path.write_text(default_filtered_vcf, encoding="utf-8") + return subprocess.CompletedProcess(cmd, 0) + + monkeypatch.setattr("vartracker.lofreq_primer_rescue.subprocess.run", fake_run) + + output_vcf = tmp_path / "final.vcf" + filtered_out_tsv = tmp_path / "filtered_out.tsv" + result = lofreq_filter_with_audit( + raw_vcf=raw_vcf, + output_vcf=output_vcf, + filtered_out_tsv=filtered_out_tsv, + ) + + assert result.normal_passed == 1 + assert result.filtered_out == 1 + output_text = output_vcf.read_text(encoding="utf-8") + assert "chr1\t10\t.\tA\tG\t200\tPASS" in output_text + assert "chr1\t50\t.\tT\tC" not in output_text + filtered_out_lines = filtered_out_tsv.read_text(encoding="utf-8").splitlines() + assert filtered_out_lines[0] == "variant\treason_filtered\tmetrics" + assert filtered_out_lines[1].startswith("T50C\tsb_fdr\tAF=0.96") diff --git a/tests/test_main.py b/tests/test_main.py index b2ae61b..3885c39 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -54,12 +54,15 @@ def test_bam_parser_accepts_lofreq_primer_rescue_options(): "off", "--lofreq-rescue-min-af", "0.9", + "--lofreq-rescue-max-minor-alt-fraction", + "0.1", ] ) assert args.primer_bed == "primers.bed" assert args.lofreq_primer_rescue == "off" assert args.lofreq_rescue_min_af == 0.9 + assert args.lofreq_rescue_max_minor_alt_fraction == 0.1 def test_drop_exact_duplicate_result_rows_removes_only_exact_duplicates(capsys): @@ -1072,6 +1075,7 @@ def fake_vcf(args): assert recorded["workflow_kwargs"]["ampliconclip_tolerance"] == 2 assert recorded["workflow_kwargs"]["lofreq_primer_rescue"] == "auto" assert recorded["workflow_kwargs"]["lofreq_rescue_min_af"] == 0.95 + assert recorded["workflow_kwargs"]["lofreq_rescue_max_minor_alt_fraction"] == 0.05 assert recorded["vcf_input"] == str(updated_csv) assert modes_checked == ["e2e"] diff --git a/tests/test_snakemake_workflow.py b/tests/test_snakemake_workflow.py index b100cb9..a29e3b3 100644 --- a/tests/test_snakemake_workflow.py +++ b/tests/test_snakemake_workflow.py @@ -33,7 +33,23 @@ def test_snakemake_rules_write_logs_under_outdir(): assert "--tolerance {params.tolerance}" in snakefile assert "--no-default-filter" in snakefile assert "_variants.raw.vcf.gz" in snakefile + assert "_variants.raw.vcf.gz.tbi" in snakefile + assert ( + 'vcf_raw = temp(f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz")' + not in snakefile + ) + assert ( + 'vcf_raw = f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz"' in snakefile + ) + assert ( + 'vcf_raw_tbi = f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz.tbi"' + in snakefile + ) + assert "_variants.filtered_out.tsv" in snakefile assert "LOFREQ_PRIMER_RESCUE_ENABLED" in snakefile + assert "LOFREQ_RESCUE_MAX_MINOR_ALT_FRACTION" in snakefile + assert "max_minor_alt_fraction=LOFREQ_RESCUE_MAX_MINOR_ALT_FRACTION" in snakefile + assert "lofreq_filter_with_audit" in snakefile assert "rescue_lofreq_primer_variants" in snakefile assert "_variants.rescued.tsv" in snakefile assert "_validate_primer_bed_reference(PRIMER_BED, REF)" in snakefile @@ -41,4 +57,6 @@ def test_snakemake_rules_write_logs_under_outdir(): assert "_iupac_consensus.fasta" in snakefile assert "df['consensus']" in snakefile assert "df['iupac_consensus']" in snakefile + assert "df['raw_vcf']" in snakefile assert "df['lofreq_rescued_tsv']" in snakefile + assert "df['lofreq_filtered_out_tsv']" in snakefile diff --git a/tests/test_vcf_processing.py b/tests/test_vcf_processing.py index d4f589a..6816f40 100644 --- a/tests/test_vcf_processing.py +++ b/tests/test_vcf_processing.py @@ -401,8 +401,8 @@ def test_process_vcf_splits_sample_specific_bcsq_annotations(tmp_path): "start": 5, "amino_acid_consequence": "K2T", "bcsq_aa_notation": "2K>2T", - "presence_absence": "N / Y", - "alt_freq": ". / 0.100", + "presence_absence": "Y / Y", + "alt_freq": "0.100 / 0.100", }, ] ) @@ -410,6 +410,85 @@ def test_process_vcf_splits_sample_specific_bcsq_annotations(tmp_path): pd.testing.assert_frame_equal(observed, expected) +def test_process_vcf_merges_starred_and_unstarred_equivalent_bcsq(tmp_path): + vcf_path = tmp_path / "annotated.vcf" + vcf_path.write_text( + "##fileformat=VCFv4.2\n" + "##contig=\n" + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\ts2\n" + "chr1\t4\t.\tA\tG\t.\tPASS\tDP=100;AF=0.1;" + "BCSQ=*synonymous|GENE1|tx|protein_coding|+|2K|4A>G," + "synonymous|GENE1|tx|protein_coding|+|2K|4A>G" + "\tGT:DP:AF:BCSQ\t1:100:0.1:1\t1:100:0.2:4\n", + encoding="utf-8", + ) + + cov1 = tmp_path / "s1.depth.txt" + cov2 = tmp_path / "s2.depth.txt" + _write_depth_file(cov1) + _write_depth_file(cov2) + + table = process_vcf(str(vcf_path), [str(cov1), str(cov2)], 10, ["s1", "s2"]) + + assert len(table) == 1 + row = table.iloc[0] + assert row["variant"] == "A4G" + assert row["type_of_change"] == "synonymous" + assert row["presence_absence"] == "Y / Y" + assert row["alt_freq"] == "0.100 / 0.200" + assert row["variant_status"] == "original" + assert row["persistence_status"] == "original_retained" + + +def test_process_vcf_keeps_single_site_variant_present_when_sample_has_joint_csq( + tmp_path, +): + vcf_path = tmp_path / "annotated.vcf" + vcf_path.write_text( + "##fileformat=VCFv4.2\n" + "##contig=\n" + '##INFO=\n' + '##INFO=\n' + '##INFO=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + '##FORMAT=\n' + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\ts2\n" + "chr1\t4\t.\tA\tG\t.\tPASS\tDP=100;AF=0.2;" + "BCSQ=missense|GENE1|tx|protein_coding|+|2K>2A|4A>G," + "frameshift|GENE1|tx|protein_coding|+|2KAAAAAAAAAAAA>2K|4A>G+5A>C" + "\tGT:DP:AF:BCSQ\t1:100:0.1:1\t1:100:0.2:4\n", + encoding="utf-8", + ) + + cov1 = tmp_path / "s1.depth.txt" + cov2 = tmp_path / "s2.depth.txt" + _write_depth_file(cov1) + _write_depth_file(cov2) + + table = process_vcf(str(vcf_path), [str(cov1), str(cov2)], 10, ["s1", "s2"]) + simple = table[table["bcsq_nt_notation"].eq("4A>G")].iloc[0] + joint = table[table["bcsq_nt_notation"].str.contains("\\+", regex=True)].iloc[0] + + assert simple["type_of_change"] == "missense" + assert simple["presence_absence"] == "Y / Y" + assert simple["alt_freq"] == "0.100 / 0.200" + assert simple["variant_status"] == "original" + assert simple["persistence_status"] == "original_retained" + + assert joint["type_of_change"] == "frameshift" + assert joint["presence_absence"] == "N / Y" + assert joint["alt_freq"] == ". / 0.200" + + @pytest.mark.skipif(shutil.which("bcftools") is None, reason="bcftools not available") def test_merge_then_annotate_preserves_joint_annotations_across_samples(tmp_path): ref, gff = _write_minimal_reference_bundle(tmp_path) @@ -525,7 +604,7 @@ def test_merge_then_annotate_preserves_joint_annotations_across_samples(tmp_path "start": 5, "variant": "A5C", "amino_acid_consequence": "K2T", - "presence_absence": "N / Y", + "presence_absence": "Y / Y", "type_of_change": "missense", "joint_variant": False, }, diff --git a/vartracker/Snakefile b/vartracker/Snakefile index 1da3c2e..d5700ed 100644 --- a/vartracker/Snakefile +++ b/vartracker/Snakefile @@ -4,12 +4,14 @@ try: from vartracker.consensus import consensus_genotype_for_variant from vartracker.lofreq_primer_rescue import ( PrimerRescueThresholds, + lofreq_filter_with_audit, rescue_lofreq_primer_variants, ) except ModuleNotFoundError: from consensus import consensus_genotype_for_variant from lofreq_primer_rescue import ( PrimerRescueThresholds, + lofreq_filter_with_audit, rescue_lofreq_primer_variants, ) @@ -33,6 +35,9 @@ LOFREQ_RESCUE_MIN_DP = int(config.get("lofreq_rescue_min_dp", 100)) LOFREQ_RESCUE_MIN_ALT_COUNT = int(config.get("lofreq_rescue_min_alt_count", 95)) LOFREQ_RESCUE_MIN_QUAL = float(config.get("lofreq_rescue_min_qual", 100.0)) LOFREQ_RESCUE_MAX_REF_COUNT = int(config.get("lofreq_rescue_max_ref_count", 20)) +LOFREQ_RESCUE_MAX_MINOR_ALT_FRACTION = float( + config.get("lofreq_rescue_max_minor_alt_fraction", 0.05) +) if LOFREQ_PRIMER_RESCUE not in {"auto", "on", "off"}: raise ValueError("lofreq_primer_rescue must be one of: auto, on, off") @@ -120,6 +125,10 @@ rule all: input: expand(f"{OUTDIR}/{{sample}}/{{sample}}_variants.vcf.gz", sample=SAMPLES.keys()), expand(f"{OUTDIR}/{{sample}}/{{sample}}_variants.vcf.gz.csi", sample=SAMPLES.keys()), + expand(f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz", sample=SAMPLES.keys()), + expand(f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz.tbi", sample=SAMPLES.keys()), + expand(f"{OUTDIR}/{{sample}}/{{sample}}_variants.rescued.tsv", sample=SAMPLES.keys()), + expand(f"{OUTDIR}/{{sample}}/{{sample}}_variants.filtered_out.tsv", sample=SAMPLES.keys()), expand(f"{OUTDIR}/{{sample}}/{{sample}}_depth.txt", sample=SAMPLES.keys()), expand(f"{OUTDIR}/{{sample}}/{{sample}}_consensus.fasta", sample=SAMPLES.keys()), expand(f"{OUTDIR}/{{sample}}/{{sample}}_iupac_consensus.fasta", sample=SAMPLES.keys()), @@ -279,9 +288,11 @@ rule lofreq_call: bam = f"{OUTDIR}/{{sample}}/{{sample}}_aligned.indelqual.bam", ref = REF output: - vcf_raw = temp(f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz"), + vcf_raw = f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz", + vcf_raw_tbi = f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz.tbi", vcf_filtered = temp(f"{OUTDIR}/{{sample}}/{{sample}}_variants.filtered.vcf"), rescued_tsv = f"{OUTDIR}/{{sample}}/{{sample}}_variants.rescued.tsv", + filtered_out_tsv = f"{OUTDIR}/{{sample}}/{{sample}}_variants.filtered_out.tsv", vcf = f"{OUTDIR}/{{sample}}/{{sample}}_variants.vcf.gz", csi = f"{OUTDIR}/{{sample}}/{{sample}}_variants.vcf.gz.csi" log: @@ -304,12 +315,14 @@ rule lofreq_call: min_alt_count=LOFREQ_RESCUE_MIN_ALT_COUNT, min_qual=LOFREQ_RESCUE_MIN_QUAL, max_ref_count=LOFREQ_RESCUE_MAX_REF_COUNT, + max_minor_alt_fraction=LOFREQ_RESCUE_MAX_MINOR_ALT_FRACTION, ) result = rescue_lofreq_primer_variants( raw_vcf=str(output.vcf_raw), primers_bed=PRIMER_BED, output_vcf=str(output.vcf_filtered), rescued_tsv=str(output.rescued_tsv), + filtered_out_tsv=str(output.filtered_out_tsv), tmp_dir=str(log_path.parent), thresholds=thresholds, ) @@ -317,17 +330,24 @@ rule lofreq_call: handle.write( "lofreq primer rescue: " f"normal_passed={result.normal_passed} " - f"rescued={result.rescued} discarded={result.discarded}\n" + f"rescued={result.rescued} discarded={result.discarded} " + f"filtered_out={result.filtered_out}\n" ) else: - shell( - "lofreq filter -i {output.vcf_raw} " - "-o {output.vcf_filtered} 2>> {log}" + result = lofreq_filter_with_audit( + raw_vcf=str(output.vcf_raw), + output_vcf=str(output.vcf_filtered), + filtered_out_tsv=str(output.filtered_out_tsv), + tmp_dir=str(log_path.parent), ) with open(output.rescued_tsv, "w", encoding="utf-8") as handle: handle.write("variant\treason_filtered\treason_rescued\tmetrics\n") with log_path.open("a", encoding="utf-8") as handle: - handle.write("lofreq primer rescue: disabled\n") + handle.write( + "lofreq primer rescue: disabled " + f"normal_passed={result.normal_passed} " + f"filtered_out={result.filtered_out}\n" + ) shell("bgzip -c {output.vcf_filtered} > {output.vcf} 2>> {log}") shell("bcftools index -f {output.vcf} 2>> {log}") @@ -597,12 +617,18 @@ rule update_csv: df['vcf'] = df['sample_name'].apply( lambda x: os.path.abspath(f"{params.outdir}/{x}/{x}_variants.vcf.gz") ) + df['raw_vcf'] = df['sample_name'].apply( + lambda x: os.path.abspath(f"{params.outdir}/{x}/{x}_variants.raw.vcf.gz") + ) df['coverage'] = df['sample_name'].apply( lambda x: os.path.abspath(f"{params.outdir}/{x}/{x}_depth.txt") ) df['lofreq_rescued_tsv'] = df['sample_name'].apply( lambda x: os.path.abspath(f"{params.outdir}/{x}/{x}_variants.rescued.tsv") ) + df['lofreq_filtered_out_tsv'] = df['sample_name'].apply( + lambda x: os.path.abspath(f"{params.outdir}/{x}/{x}_variants.filtered_out.tsv") + ) df['consensus'] = df['sample_name'].apply( lambda x: os.path.abspath(f"{params.outdir}/{x}/{x}_consensus.fasta") ) diff --git a/vartracker/_version.py b/vartracker/_version.py index 4e0c032..11e633a 100644 --- a/vartracker/_version.py +++ b/vartracker/_version.py @@ -9,7 +9,7 @@ # NOTE: When bumping the project version remember to update this fallback value # alongside the version declared in pyproject.toml. -_FALLBACK_VERSION = "2.2.0" +_FALLBACK_VERSION = "2.2.1" try: __version__ = metadata.version("vartracker") diff --git a/vartracker/analysis.py b/vartracker/analysis.py index 9a395ad..e895d93 100644 --- a/vartracker/analysis.py +++ b/vartracker/analysis.py @@ -34,7 +34,7 @@ def _ensure_joint_prefix(change_type: object) -> str: text = str(change_type or "").strip() if not text: return "joint" - normalised = re.sub(r"^(joint_)+", "", text) + normalised = re.sub(r"^(joint_)+", "", text).lstrip("*") return f"joint_{normalised}" if normalised else "joint" @@ -493,20 +493,35 @@ def _build_gene_order_map( return gene_order_map, ordered -def _map_orf1ab_position_to_nsp(aa_position: int) -> str: - """Map an ORF1ab amino acid position to its NSP name.""" +def _map_orf1ab_position_to_nsp_info(aa_position: int) -> tuple[str, int]: + """Map an ORF1ab amino acid position to its NSP name and local position.""" nsp_products = cast(Sequence[str], NSPS.get("product", [])) aa_starts = cast(Sequence[int], NSPS.get("aa_start", [])) - for name, start in zip(nsp_products, aa_starts): - length = NSP_LENGTHS.get(str(name)) - if length is None: - continue - end = start + length - 1 + for idx, (name, start) in enumerate(zip(nsp_products, aa_starts)): + end = aa_starts[idx + 1] - 1 if idx + 1 < len(aa_starts) else 999999 if start <= aa_position <= end: - return name + return str(name), aa_position - start + 1 + + return "ORF1ab", aa_position + + +def _map_orf1ab_position_to_nsp(aa_position: int) -> str: + """Map an ORF1ab amino acid position to its NSP name.""" + return _map_orf1ab_position_to_nsp_info(aa_position)[0] + - return "ORF1ab" +def _format_orf1ab_change_as_nsp(amino_change: object) -> tuple[str, str] | None: + """Convert an ORF1ab amino-acid label to NSP-local notation when possible.""" + text = str(amino_change or "").strip() + match = re.fullmatch(r"([A-Za-z*]*)(\d+)([A-Za-z*]*)", text) + if not match: + return None + ref = match.group(1) + aa_position = int(match.group(2)) + alt = match.group(3) + gene_label, local_position = _map_orf1ab_position_to_nsp_info(aa_position) + return gene_label, f"{ref}{local_position}{alt}" def _extract_numeric_position(value: str) -> Optional[int]: @@ -547,7 +562,8 @@ def _resolve_variant_labels(row) -> Tuple[str, str, str]: gene = getattr(row, "gene", "") amino_change = getattr(row, "amino_acid_consequence", "") nsp_change = getattr(row, "nsp_aa_change", "") - change_type = str(getattr(row, "type_of_change", "")) + change_type = str(getattr(row, "type_of_change", "")).lstrip("*") + change_type = re.sub(r"^(joint_)+", "", change_type) gene_label = gene aa_label = amino_change @@ -559,10 +575,14 @@ def _resolve_variant_labels(row) -> Tuple[str, str, str]: gene_label = gene_part if aa_part: aa_label = aa_part + if "*" in str(amino_change) and "*" not in str(aa_label): + mapped_change = _format_orf1ab_change_as_nsp(amino_change) + if mapped_change is not None: + gene_label, aa_label = mapped_change else: - position = _extract_numeric_position(amino_change) - if position is not None: - gene_label = _map_orf1ab_position_to_nsp(position) + mapped_change = _format_orf1ab_change_as_nsp(amino_change) + if mapped_change is not None: + gene_label, aa_label = mapped_change if not aa_label or str(aa_label) in {"", "None"}: if isinstance(nsp_change, str) and ":" in nsp_change: @@ -710,7 +730,7 @@ def _prepare_variant_heatmap_matrix( ) records: List[Dict[str, Union[str, float, int]]] = [] - seen_labels = set() + record_index_by_label: dict[str, int] = {} qc_maps: dict[str, dict[str, str]] = {} for row in table.itertuples(index=False): @@ -729,10 +749,9 @@ def _prepare_variant_heatmap_matrix( gene_label, display_label, base_label = _resolve_variant_labels(row) - if base_label in seen_labels: - continue - - change_type = str(getattr(row, "type_of_change", "")).strip().lower() + change_type = ( + str(getattr(row, "type_of_change", "")).strip().lower().lstrip("*") + ) if not include_joint and change_type.startswith("joint"): continue if included_patterns and not any( @@ -842,8 +861,33 @@ def _prepare_variant_heatmap_matrix( for sample, value in zip(ordered_samples, row_values): record[sample] = value + if display_label in record_index_by_label: + existing = records[record_index_by_label[display_label]] + existing["gene_order"] = min( + int(existing.get("gene_order", record["gene_order"])), + int(record["gene_order"]), + ) + existing["start"] = min( + int(existing.get("start", record["start"])), + int(record["start"]), + ) + existing_qc = qc_maps.setdefault(display_label, {}) + for sample, value in zip(ordered_samples, row_values): + current_value = existing.get(sample, 0.0) + current_frequency = ( + float(current_value) + if isinstance(current_value, (int, float)) + else _coerce_frequency(str(current_value)) + ) + if value > current_frequency: + existing[sample] = value + existing_qc[sample] = row_qc_map.get(sample, "") + elif sample not in existing_qc: + existing_qc[sample] = row_qc_map.get(sample, "") + continue + + record_index_by_label[display_label] = len(records) records.append(record) - seen_labels.add(base_label) qc_maps[display_label] = { sample: row_qc_map.get(sample, "") for sample in ordered_samples } diff --git a/vartracker/analysis_launcher.py b/vartracker/analysis_launcher.py index 791fe68..28d3e15 100644 --- a/vartracker/analysis_launcher.py +++ b/vartracker/analysis_launcher.py @@ -75,6 +75,7 @@ def _validate_lofreq_primer_rescue( min_alt_count: int, min_qual: float, max_ref_count: int, + max_minor_alt_fraction: float, ) -> None: if mode not in {"auto", "on", "off"}: raise ValueError("lofreq_primer_rescue must be one of: auto, on, off") @@ -90,6 +91,8 @@ def _validate_lofreq_primer_rescue( raise ValueError("lofreq_rescue_min_qual must be >= 0") if max_ref_count < 0: raise ValueError("lofreq_rescue_max_ref_count must be >= 0") + if not 0 <= max_minor_alt_fraction <= 1: + raise ValueError("lofreq_rescue_max_minor_alt_fraction must be between 0 and 1") def run_workflow( @@ -114,6 +117,7 @@ def run_workflow( lofreq_rescue_min_alt_count: int = 95, lofreq_rescue_min_qual: float = 100.0, lofreq_rescue_max_ref_count: int = 20, + lofreq_rescue_max_minor_alt_fraction: float = 0.05, ) -> Optional[str]: """Run the lofreq variant calling workflow via the Snakemake API. @@ -146,6 +150,7 @@ def run_workflow( lofreq_rescue_min_alt_count, lofreq_rescue_min_qual, lofreq_rescue_max_ref_count, + lofreq_rescue_max_minor_alt_fraction, ) Path(outdir).mkdir(parents=True, exist_ok=True) @@ -167,6 +172,7 @@ def run_workflow( "lofreq_rescue_min_alt_count": lofreq_rescue_min_alt_count, "lofreq_rescue_min_qual": lofreq_rescue_min_qual, "lofreq_rescue_max_ref_count": lofreq_rescue_max_ref_count, + "lofreq_rescue_max_minor_alt_fraction": (lofreq_rescue_max_minor_alt_fraction), } if primer_bed_path: config_dict["primer_bed"] = primer_bed_path diff --git a/vartracker/constants.py b/vartracker/constants.py index 172314e..3fb663f 100644 --- a/vartracker/constants.py +++ b/vartracker/constants.py @@ -145,24 +145,16 @@ def bcf_orf1ab_to_nsp(mutation): Raises: ValueError: If mutation format is invalid """ - # Parse out important bits from mutation + # Parse out important bits from mutation. Stop codons are represented as "*" + # and must be retained in both ORF1ab and NSP notation. no_gene = re.sub(".*:", "", mutation) - parsed = re.findall(r"([A-Z]+)", no_gene) - - if len(parsed) == 2: - ref = parsed[0] - alt = parsed[1] - elif len(parsed) == 1: - ref = "" - alt = parsed[0] - else: + parsed = re.fullmatch(r"([A-Z*]*)(\d+)([A-Z*]*)", no_gene) + if not parsed: raise ValueError(f"Invalid mutation format: {mutation}") - pos_match = re.search(r"(\d+)", no_gene) - if not pos_match: - raise ValueError(f"No position found in mutation: {mutation}") - - pos = int(pos_match.group()) + ref = parsed.group(1) + pos = int(parsed.group(2)) + alt = parsed.group(3) # Find corresponding nsp idx = bisect_right(NSPS["aa_start"], pos) - 1 @@ -264,9 +256,18 @@ def reformat_csq_notation(gene, string): if splitter == -1: return (string, "") - ref = re.sub("[0-9]", "", string[:splitter]) - pos = re.sub("[A-Za-z]", "", string[:splitter]) - alt = re.sub("[0-9]", "", string[splitter + 1 :]) + left = string[:splitter] + right = string[splitter + 1 :] + left_match = re.fullmatch(r"(\d+)([A-Za-z*]+)", left) + right_match = re.fullmatch(r"(\d+)([A-Za-z*]+)", right) + if left_match and right_match: + ref = left_match.group(2) + pos = left_match.group(1) + alt = right_match.group(2) + else: + ref = re.sub("[0-9]", "", left) + pos = re.sub("[A-Za-z*]", "", left) + alt = re.sub("[0-9]", "", right) reformatted = ref + pos + alt if gene == "ORF1ab": diff --git a/vartracker/lofreq_primer_rescue.py b/vartracker/lofreq_primer_rescue.py index d0a80cb..0cbaaec 100644 --- a/vartracker/lofreq_primer_rescue.py +++ b/vartracker/lofreq_primer_rescue.py @@ -35,6 +35,7 @@ class PrimerRescueThresholds: min_alt_count: int = 95 min_qual: float = 100.0 max_ref_count: int = 20 + max_minor_alt_fraction: float = 0.05 @dataclass(frozen=True) @@ -45,6 +46,17 @@ class PrimerRescueResult: rescued: int discarded: int rescued_tsv: str + filtered_out: int = 0 + filtered_out_tsv: str | None = None + + +@dataclass(frozen=True) +class LofreqFilterAuditResult: + """Summary of default LoFreq filtering with an audit table.""" + + normal_passed: int + filtered_out: int + filtered_out_tsv: str def _open_text(path: str | Path, mode: str = "rt") -> TextIO: @@ -126,6 +138,13 @@ def _filter_is_pass(fields: list[str]) -> bool: return fields[6] in ("PASS", ".") +def _is_strand_bias_filter(reason: str) -> bool: + return any( + item.lower() in {"sb_fdr", "strandbias", "strand_bias"} + for item in reason.split(";") + ) + + def read_default_filter_results( path: str | Path, ) -> tuple[set[tuple[str, str, str, str]], dict[tuple[str, str, str, str], str]]: @@ -150,9 +169,13 @@ def rescue_metrics( fields: list[str], primers: dict[str, list[tuple[int, int]]], thresholds: PrimerRescueThresholds, + reason_filtered: str, ) -> str | None: """Return rescue metrics for a candidate record, or ``None`` if it fails.""" + if not _is_strand_bias_filter(reason_filtered): + return None + chrom = fields[0] pos = int(fields[1]) ref = fields[3] @@ -206,17 +229,52 @@ def rescue_metrics( return None if ref_count > thresholds.max_ref_count: return None - if not (alt_fwd == 0 or alt_rev == 0): + minor_alt_fraction = min(alt_fwd, alt_rev) / alt_count + if minor_alt_fraction > thresholds.max_minor_alt_fraction: return None dp4_af = alt_count / dp4_total return ( f"AF={af:g},DP={dp},QUAL={qual:g}," f"DP4={ref_fwd}/{ref_rev}/{alt_fwd}/{alt_rev}," - f"alt_count={alt_count},ref_count={ref_count},dp4_af={dp4_af:g}" + f"alt_count={alt_count},ref_count={ref_count},dp4_af={dp4_af:g}," + f"minor_alt_fraction={minor_alt_fraction:g}" ) +def variant_metrics(fields: list[str]) -> str: + """Return concise metrics for a raw LoFreq record.""" + if len(fields) < 8: + return "." + + info = parse_info(fields[7]) + metrics = [] + if "AF" in info: + try: + metrics.append(f"AF={_first_float(info['AF']):g}") + except ValueError: + metrics.append(f"AF={info['AF']}") + if "DP" in info: + try: + metrics.append(f"DP={_first_int(info['DP'])}") + except ValueError: + metrics.append(f"DP={info['DP']}") + if fields[5] != ".": + metrics.append(f"QUAL={fields[5]}") + counts = _dp4_counts(info.get("DP4", "")) + if counts is not None: + ref_fwd, ref_rev, alt_fwd, alt_rev = counts + ref_count = ref_fwd + ref_rev + alt_count = alt_fwd + alt_rev + total = ref_count + alt_count + metrics.append(f"DP4={ref_fwd}/{ref_rev}/{alt_fwd}/{alt_rev}") + metrics.append(f"alt_count={alt_count}") + metrics.append(f"ref_count={ref_count}") + if total: + metrics.append(f"dp4_af={alt_count / total:g}") + return ",".join(metrics) if metrics else "." + + def _add_info_flag(info_text: str, flag: str) -> str: if info_text in ("", "."): return flag @@ -284,71 +342,186 @@ def _write_headers(headers: list[str], output: TextIO) -> None: output.write(line) -def write_final_vcf( +def _run_lofreq_filter_print_all( + raw_vcf: str | Path, + tmp_dir: str | Path | None = None, + lofreq: str = "lofreq", +) -> tuple[set[tuple[str, str, str, str]], dict[tuple[str, str, str, str], str]]: + with tempfile.TemporaryDirectory(dir=tmp_dir) as run_tmp_dir: + filtered_vcf = os.path.join(run_tmp_dir, "lofreq.default_filtered.vcf") + subprocess.run( + [lofreq, "filter", "--print-all", "-i", str(raw_vcf), "-o", filtered_vcf], + check=True, + ) + return read_default_filter_results(filtered_vcf) + + +def write_default_filtered_vcf( raw_vcf: str | Path, output_vcf: str | Path, - rescued_tsv_path: str | Path, + filtered_out_tsv_path: str | Path, pass_keys: set[tuple[str, str, str, str]], filter_reasons: dict[tuple[str, str, str, str], str], - primers: dict[str, list[tuple[int, int]]], - thresholds: PrimerRescueThresholds, -) -> PrimerRescueResult: - """Write the default-filtered VCF plus any rescued primer-overlap records.""" +) -> LofreqFilterAuditResult: + """Write normal LoFreq PASS records plus a filtered-out audit table.""" normal_passed = 0 - rescued = 0 - discarded = 0 + filtered_out = 0 headers: list[str] = [] with ( _open_text(raw_vcf) as raw, _open_text(output_vcf, "wt") as output, - Path(rescued_tsv_path).open("w", encoding="utf-8") as rescued_tsv, + Path(filtered_out_tsv_path).open("w", encoding="utf-8") as filtered_out_tsv, ): - rescued_tsv.write("variant\treason_filtered\treason_rescued\tmetrics\n") - + filtered_out_tsv.write("variant\treason_filtered\tmetrics\n") for line in raw: if line.startswith("#"): headers.append(line) continue if headers: - _write_headers(headers, output) + output.writelines(headers) headers = [] fields = line.rstrip("\n").split("\t") if len(fields) < 8: - discarded += 1 continue key = _record_key(fields) if key in pass_keys: output.write("\t".join(_pass_fields(fields)) + "\n") normal_passed += 1 - else: - metrics = rescue_metrics(fields, primers, thresholds) - if metrics is None: + continue + + reason_filtered = filter_reasons.get( + key, "not_passed_default_lofreq_filter" + ) + filtered_out_tsv.write( + f"{_variant_name(fields)}\t{reason_filtered}\t" + f"{variant_metrics(fields)}\n" + ) + filtered_out += 1 + + if headers: + output.writelines(headers) + + return LofreqFilterAuditResult( + normal_passed=normal_passed, + filtered_out=filtered_out, + filtered_out_tsv=str(filtered_out_tsv_path), + ) + + +def write_final_vcf( + raw_vcf: str | Path, + output_vcf: str | Path, + rescued_tsv_path: str | Path, + filtered_out_tsv_path: str | Path | None, + pass_keys: set[tuple[str, str, str, str]], + filter_reasons: dict[tuple[str, str, str, str], str], + primers: dict[str, list[tuple[int, int]]], + thresholds: PrimerRescueThresholds, +) -> PrimerRescueResult: + """Write the default-filtered VCF plus any rescued primer-overlap records.""" + + normal_passed = 0 + rescued = 0 + discarded = 0 + filtered_out = 0 + headers: list[str] = [] + + with ( + _open_text(raw_vcf) as raw, + _open_text(output_vcf, "wt") as output, + Path(rescued_tsv_path).open("w", encoding="utf-8") as rescued_tsv, + ): + rescued_tsv.write("variant\treason_filtered\treason_rescued\tmetrics\n") + filtered_out_tsv = ( + Path(filtered_out_tsv_path).open("w", encoding="utf-8") + if filtered_out_tsv_path is not None + else None + ) + if filtered_out_tsv is not None: + filtered_out_tsv.write("variant\treason_filtered\tmetrics\n") + + try: + for line in raw: + if line.startswith("#"): + headers.append(line) + continue + + if headers: + _write_headers(headers, output) + headers = [] + + fields = line.rstrip("\n").split("\t") + if len(fields) < 8: discarded += 1 continue - output.write("\t".join(_rescued_fields(fields)) + "\n") - reason_filtered = filter_reasons.get( - key, "not_passed_default_lofreq_filter" - ) - rescued_tsv.write( - f"{_variant_name(fields)}\t{reason_filtered}\t" - f"overlap_primer_interval\t{metrics}\n" - ) - rescued += 1 + key = _record_key(fields) + if key in pass_keys: + output.write("\t".join(_pass_fields(fields)) + "\n") + normal_passed += 1 + else: + reason_filtered = filter_reasons.get( + key, "not_passed_default_lofreq_filter" + ) + metrics = rescue_metrics( + fields, primers, thresholds, reason_filtered + ) + if metrics is None: + discarded += 1 + if filtered_out_tsv is not None: + filtered_out_tsv.write( + f"{_variant_name(fields)}\t{reason_filtered}\t" + f"{variant_metrics(fields)}\n" + ) + filtered_out += 1 + continue + + output.write("\t".join(_rescued_fields(fields)) + "\n") + rescued_tsv.write( + f"{_variant_name(fields)}\t{reason_filtered}\t" + f"overlap_primer_interval\t{metrics}\n" + ) + rescued += 1 - if headers: - _write_headers(headers, output) + if headers: + _write_headers(headers, output) + finally: + if filtered_out_tsv is not None: + filtered_out_tsv.close() return PrimerRescueResult( normal_passed=normal_passed, rescued=rescued, discarded=discarded, rescued_tsv=str(rescued_tsv_path), + filtered_out=filtered_out, + filtered_out_tsv=( + str(filtered_out_tsv_path) if filtered_out_tsv_path is not None else None + ), + ) + + +def lofreq_filter_with_audit( + raw_vcf: str | Path, + output_vcf: str | Path, + filtered_out_tsv: str | Path, + tmp_dir: str | Path | None = None, + lofreq: str = "lofreq", +) -> LofreqFilterAuditResult: + """Apply default LoFreq filtering and write filtered-out variant details.""" + + pass_keys, filter_reasons = _run_lofreq_filter_print_all(raw_vcf, tmp_dir, lofreq) + return write_default_filtered_vcf( + raw_vcf, + output_vcf, + filtered_out_tsv, + pass_keys, + filter_reasons, ) @@ -357,6 +530,7 @@ def rescue_lofreq_primer_variants( primers_bed: str | Path, output_vcf: str | Path, rescued_tsv: str | Path | None = None, + filtered_out_tsv: str | Path | None = None, tmp_dir: str | Path | None = None, lofreq: str = "lofreq", thresholds: PrimerRescueThresholds | None = None, @@ -366,21 +540,18 @@ def rescue_lofreq_primer_variants( thresholds = thresholds or PrimerRescueThresholds() if rescued_tsv is None: rescued_tsv = f"{output_vcf}.rescued.tsv" + if filtered_out_tsv is None: + filtered_out_tsv = f"{output_vcf}.filtered_out.tsv" primers = load_primers(primers_bed) - with tempfile.TemporaryDirectory(dir=tmp_dir) as run_tmp_dir: - filtered_vcf = os.path.join(run_tmp_dir, "lofreq.default_filtered.vcf") - subprocess.run( - [lofreq, "filter", "-i", str(raw_vcf), "-o", filtered_vcf], - check=True, - ) - pass_keys, filter_reasons = read_default_filter_results(filtered_vcf) + pass_keys, filter_reasons = _run_lofreq_filter_print_all(raw_vcf, tmp_dir, lofreq) return write_final_vcf( raw_vcf, output_vcf, rescued_tsv, + filtered_out_tsv, pass_keys, filter_reasons, primers, @@ -404,6 +575,10 @@ def _parse_args(argv: list[str] | None = None) -> argparse.Namespace: "--rescued-tsv", help="TSV of rescued variants; default: OUTPUT_VCF.rescued.tsv", ) + parser.add_argument( + "--filtered-out-tsv", + help="TSV of raw variants retained by LoFreq calling but filtered from final VCF", + ) parser.add_argument("--tmp-dir", help="Optional temporary directory") parser.add_argument("--lofreq", default="lofreq", help="lofreq executable") parser.add_argument( @@ -436,6 +611,15 @@ def _parse_args(argv: list[str] | None = None) -> argparse.Namespace: default=PrimerRescueThresholds.max_ref_count, help="Maximum DP4 ref count for rescue candidates only (default: 20)", ) + parser.add_argument( + "--max-minor-alt-fraction", + type=float, + default=PrimerRescueThresholds.max_minor_alt_fraction, + help=( + "Maximum minor ALT strand fraction for rescue candidates only " + "(default: 0.05)" + ), + ) return parser.parse_args(argv) @@ -447,19 +631,21 @@ def main(argv: list[str] | None = None) -> int: min_alt_count=args.min_alt_count, min_qual=args.min_qual, max_ref_count=args.max_ref_count, + max_minor_alt_fraction=args.max_minor_alt_fraction, ) result = rescue_lofreq_primer_variants( raw_vcf=args.raw_vcf, primers_bed=args.primers_bed, output_vcf=args.output_vcf, rescued_tsv=args.rescued_tsv, + filtered_out_tsv=args.filtered_out_tsv, tmp_dir=args.tmp_dir, lofreq=args.lofreq, thresholds=thresholds, ) print( f"normal_passed={result.normal_passed} rescued={result.rescued} " - f"discarded={result.discarded}", + f"discarded={result.discarded} filtered_out={result.filtered_out}", file=sys.stderr, ) return 0 diff --git a/vartracker/main.py b/vartracker/main.py index c9fb2c3..ddbc2b1 100644 --- a/vartracker/main.py +++ b/vartracker/main.py @@ -891,6 +891,15 @@ def _add_lofreq_primer_rescue_arguments(group: argparse._ArgumentGroup) -> None: default=PrimerRescueThresholds.max_ref_count, help="Maximum DP4 reference count for primer rescue candidates only (default: 20)", ) + group.add_argument( + "--lofreq-rescue-max-minor-alt-fraction", + type=float, + default=PrimerRescueThresholds.max_minor_alt_fraction, + help=( + "Maximum minor ALT strand fraction for primer rescue candidates only " + "(default: 0.05)" + ), + ) def _move_action_group_after( @@ -2245,6 +2254,9 @@ def _run_e2e_command(args): lofreq_rescue_min_alt_count=args.lofreq_rescue_min_alt_count, lofreq_rescue_min_qual=args.lofreq_rescue_min_qual, lofreq_rescue_max_ref_count=args.lofreq_rescue_max_ref_count, + lofreq_rescue_max_minor_alt_fraction=( + args.lofreq_rescue_max_minor_alt_fraction + ), ) if rulegraph_path: @@ -2397,6 +2409,9 @@ def _run_bam_command(args): lofreq_rescue_min_alt_count=args.lofreq_rescue_min_alt_count, lofreq_rescue_min_qual=args.lofreq_rescue_min_qual, lofreq_rescue_max_ref_count=args.lofreq_rescue_max_ref_count, + lofreq_rescue_max_minor_alt_fraction=( + args.lofreq_rescue_max_minor_alt_fraction + ), ) if rulegraph_path: diff --git a/vartracker/vcf_processing.py b/vartracker/vcf_processing.py index bf1c3e0..0c5f257 100644 --- a/vartracker/vcf_processing.py +++ b/vartracker/vcf_processing.py @@ -873,18 +873,54 @@ def _decode_sample_bcsq_annotations(v, samples, annotations): return decoded -def _mask_allele_frequencies_for_annotation( - annotation, allele_freqs, samples, sample_bcsq_map +def _normalise_bcsq_change_type(change_type): + """Remove bcftools' leading compound-context marker from a consequence type.""" + return str(change_type or "").lstrip("*") + + +def _normalise_bcsq_annotation(annotation): + """Return a BCSQ annotation key that ignores leading consequence markers.""" + parts = str(annotation).split("|", 1) + if not parts: + return str(annotation) + parts[0] = _normalise_bcsq_change_type(parts[0]) + return "|".join(parts) + + +def _group_equivalent_bcsq_annotations(annotations): + """Group BCSQ annotations that represent the same consequence.""" + groups = {} + ordered_groups = [] + for annotation in annotations: + key = _normalise_bcsq_annotation(annotation) + if key not in groups: + groups[key] = [] + ordered_groups.append(groups[key]) + groups[key].append(annotation) + return ordered_groups + + +def _representative_bcsq_annotation(annotations): + """Choose the least-decorated BCSQ annotation from an equivalent group.""" + for annotation in annotations: + change_type = str(annotation).split("|", 1)[0] + if not change_type.startswith("*"): + return annotation + return annotations[0] + + +def _mask_allele_frequencies_for_annotation_group( + annotations, allele_freqs, samples, sample_bcsq_map ): - """Keep allele frequencies only for samples where the annotation applies.""" + """Keep allele frequencies for samples matching any equivalent annotation.""" + annotation_set = set(annotations) masked = [] for sample, allele_freq in zip(samples, allele_freqs): if allele_freq == ".": masked.append(".") continue - masked.append( - allele_freq if annotation in sample_bcsq_map.get(sample, []) else "." - ) + sample_annotations = set(sample_bcsq_map.get(sample, [])) + masked.append(allele_freq if annotation_set & sample_annotations else ".") return masked @@ -980,6 +1016,29 @@ def _annotation_alt_for_record(v, anno): return None +def _annotation_is_single_record_allele(v, anno, annotation_alt) -> bool: + """Return True when a BCSQ annotation describes only the current allele.""" + if len(anno) <= 6 or (len(anno) == 1 and str(anno[0]).startswith("@")): + return False + + dna_change = str(anno[6] or "").strip() + if not dna_change or "+" in dna_change: + return False + + match = re.fullmatch(r"(\d+)([^>]+)>([^>]+)", dna_change) + if not match: + return False + + pos = int(match.group(1)) + ref = match.group(2) + alt = match.group(3) + if pos != int(v.POS) or ref != str(v.REF): + return False + if annotation_alt is not None and alt != str(annotation_alt): + return False + return alt in {str(value) for value in (v.ALT or [])} + + def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): """ Process VCF file and extract variant information. @@ -1047,11 +1106,13 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): # Process annotations if "BCSQ" in info: annotations = v.INFO["BCSQ"].split(",") + annotation_groups = _group_equivalent_bcsq_annotations(annotations) sample_bcsq_map = _decode_sample_bcsq_annotations(v, samples, annotations) produced_annotation_specific_row = False if sample_bcsq_map: - for annot in annotations: + for annotation_group in annotation_groups: + annot = _representative_bcsq_annotation(annotation_group) anno = annot.split("|") annotation_alt = _annotation_alt_for_record(v, anno) annotation_allele_freqs = ( @@ -1059,9 +1120,17 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): if annotation_alt is not None else allele_freqs ) - masked_allele_freqs = _mask_allele_frequencies_for_annotation( - annot, annotation_allele_freqs, samples, sample_bcsq_map - ) + if _annotation_is_single_record_allele(v, anno, annotation_alt): + masked_allele_freqs = annotation_allele_freqs + else: + masked_allele_freqs = ( + _mask_allele_frequencies_for_annotation_group( + annotation_group, + annotation_allele_freqs, + samples, + sample_bcsq_map, + ) + ) masked_trajectory = _summarise_sample_trajectory( masked_allele_freqs, samples ) @@ -1088,7 +1157,8 @@ def process_vcf(vcf_file, covs, min_depth, sample_names_override=None): produced_annotation_specific_row = True if not produced_annotation_specific_row: - for annot in annotations: + for annotation_group in annotation_groups: + annot = _representative_bcsq_annotation(annotation_group) anno = annot.split("|") annotation_alt = _annotation_alt_for_record(v, anno) annotation_allele_freqs = ( @@ -1179,6 +1249,7 @@ def _process_annotation( ): """Process a single annotation from bcftools csq.""" selected_alt = alt_allele or (v.ALT[0] if v.ALT else "") + type_of_change = _normalise_bcsq_change_type(anno[0] if anno else "") if len(anno) == 1 and anno[0].startswith("@"): # Joint variant annotation return { @@ -1194,7 +1265,7 @@ def _process_annotation( "bcsq_nt_notation": anno[0], "bcsq_aa_notation": anno[0], "type_of_variant": v.var_type, - "type_of_change": anno[0], + "type_of_change": type_of_change, "variant_status": variant_status, "persistence_status": persistent_status, "presence_absence": " / ".join(presence_absence), @@ -1240,7 +1311,7 @@ def _process_annotation( "bcsq_nt_notation": anno[6] if len(anno) > 5 else "", "bcsq_aa_notation": anno[5] if len(anno) > 5 else "", "type_of_variant": v.var_type, - "type_of_change": anno[0], + "type_of_change": type_of_change, "variant_status": variant_status, "persistence_status": persistent_status, "presence_absence": " / ".join(presence_absence),