From 6fae6e71a48aa8f465ec9e77459648ba4a789630 Mon Sep 17 00:00:00 2001 From: MikeWLloyd Date: Mon, 10 Aug 2026 15:06:47 -0400 Subject: [PATCH 1/3] fix #109 report low coverage edge case --- src/optitype/io/readers.py | 10 +++--- src/optitype/pipeline.py | 68 ++++++++++++++++++++++++++++++++++++++ tests/test_cli.py | 26 +++++++++++++++ tests/test_io.py | 28 ++++++++++++++++ tests/test_pipeline.py | 28 ++++++++++++++++ 5 files changed, 156 insertions(+), 4 deletions(-) diff --git a/src/optitype/io/readers.py b/src/optitype/io/readers.py index bf39e66..9ed1fa1 100644 --- a/src/optitype/io/readers.py +++ b/src/optitype/io/readers.py @@ -190,11 +190,13 @@ def pysam_to_dataframe(samfile: str) -> tuple[pd.DataFrame, pd.DataFrame]: logger.debug("\n %s %d reads loaded. Creating dataframe...", elapsed(), len(hits)) - pos_df = pd.DataFrame.from_dict(hits, orient="index") - pos_df.columns = sam.references[:] + pos_df = pd.DataFrame.from_dict(hits, orient="index", columns=sam.references[:]) - details_df = pd.DataFrame.from_dict(read_details, orient="index") - details_df.columns = ["mismatches", "read_length"] + details_df = pd.DataFrame.from_dict( + read_details, + orient="index", + columns=["mismatches", "read_length"], + ) if hit_counter > 0: logger.debug( diff --git a/src/optitype/pipeline.py b/src/optitype/pipeline.py index 1c8a744..bffa38a 100644 --- a/src/optitype/pipeline.py +++ b/src/optitype/pipeline.py @@ -68,6 +68,29 @@ class PipelineResult: output_plot: str +class LowCoverageError(RuntimeError): + """Raised when HLA typing cannot proceed due to insufficient mapped HLA reads.""" + + +def _write_low_coverage_result(out_csv: str, message: str) -> None: + """Write a structured failure TSV so callers get an explicit output artifact.""" + failed = pd.DataFrame([ + { + "A1": None, + "A2": None, + "B1": None, + "B2": None, + "C1": None, + "C2": None, + "Reads": 0, + "Objective": None, + "Status": "FAILED", + "Message": message, + } + ]) + failed.to_csv(out_csv, sep="\t", index=False) + + def get_num_threads(configured_threads: int) -> int: """Get the number of threads to use, capped by available CPUs.""" try: @@ -324,6 +347,15 @@ def run_pipeline( pos2, read_details2 = pysam_to_dataframe(bam_paths[1]) binary2 = np.sign(pos2) + if pos.empty and pos2.empty: + msg = ( + "HLA typing was not possible due to low coverage: no mapped HLA reads " + "were found in either input read file." + ) + logger.error(msg) + _write_low_coverage_result(out_csv, msg) + raise LowCoverageError(f"{msg} See {out_csv}.") + if not bam_input and config.delete_bam: os.remove(bam_paths[0]) os.remove(bam_paths[1]) @@ -361,15 +393,42 @@ def cut_last_char(x): return x[:-1] else: pos, read_details = pysam_to_dataframe(bam_paths[0]) + if pos.empty: + msg = ( + "HLA typing was not possible due to low coverage: no mapped HLA reads " + "were found in the input file." + ) + logger.error(msg) + _write_low_coverage_result(out_csv, msg) + raise LowCoverageError(f"{msg} See {out_csv}.") + if not bam_input and config.delete_bam: os.remove(bam_paths[0]) binary = np.sign(pos) + if binary.shape[0] == 0: + msg = ( + "HLA typing was not possible due to low coverage: no usable mapped HLA reads " + "remain after read pairing/selection." + ) + logger.error(msg) + _write_low_coverage_result(out_csv, msg) + raise LowCoverageError(f"{msg} See {out_csv}.") + # Filter to frequent alleles alleles_to_keep = [col for col in binary.columns if _is_frequent(col, table)] binary = binary[alleles_to_keep] + if binary.shape[1] == 0: + msg = ( + "HLA typing was not possible due to low coverage: no informative HLA allele " + "hits remained after filtering." + ) + logger.error(msg) + _write_low_coverage_result(out_csv, msg) + raise LowCoverageError(f"{msg} See {out_csv}.") + logger.debug("%s Temporary pruning of identical rows and columns", elapsed()) unique_col, _representing = ht.prune_identical_alleles(binary, report_groups=True) @@ -381,6 +440,15 @@ def cut_last_char(x): return x[:-1] minimal_alleles = ht.prune_overshadowed_alleles(temp_pruned) + if len(minimal_alleles) == 0: + msg = ( + "HLA typing was not possible due to low coverage: no minimal informative " + "alleles could be retained." + ) + logger.error(msg) + _write_low_coverage_result(out_csv, msg) + raise LowCoverageError(f"{msg} See {out_csv}.") + logger.debug("%s Keeping only the minimal number of required alleles %s", elapsed(), minimal_alleles.shape) binary = binary[minimal_alleles] diff --git a/tests/test_cli.py b/tests/test_cli.py index 6ad3031..08df001 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -6,6 +6,7 @@ from click.testing import CliRunner from optitype.cli import main +from optitype.pipeline import LowCoverageError def test_cli_help(): @@ -75,3 +76,28 @@ def test_cli_run_yara_missing_binary(tmp_path): ) assert result.exit_code != 0 assert "yara_mapper" in result.output + + +def test_cli_run_low_coverage_reports_reason_and_result_path(tmp_path): + """CLI should surface low-coverage reason and result TSV path to users.""" + fq = tmp_path / "reads.fq" + fq.write_text("@read1\nACGT\n+\nIIII\n") + + outdir = tmp_path / "out" + prefix = "lowcov" + out_csv = outdir / f"{prefix}_result.tsv" + out_csv.parent.mkdir(parents=True, exist_ok=True) + out_csv.write_text("Status\tMessage\nFAILED\tlow coverage\n") + + msg = f"HLA typing was not possible due to low coverage. See {out_csv}." + + runner = CliRunner() + with patch("optitype.pipeline.run_pipeline", side_effect=LowCoverageError(msg)): + result = runner.invoke( + main, + ["run", "-i", str(fq), "--dna", "-o", str(outdir), "--prefix", prefix], + ) + + assert result.exit_code != 0 + assert "low coverage" in result.output.lower() + assert str(out_csv) in result.output diff --git a/tests/test_io.py b/tests/test_io.py index d1b1a77..417a48a 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -4,6 +4,7 @@ import pandas as pd from optitype.io.data import get_data_path, load_reference_data, get_reference_fasta +from optitype.io import readers def test_get_data_path(): @@ -48,3 +49,30 @@ def test_get_reference_fasta_invalid(): """Test that invalid seq_type raises ValueError.""" with pytest.raises(ValueError): get_reference_fasta("invalid") + + +def test_pysam_to_dataframe_empty_alignments(monkeypatch): + """Empty BAM/SAM input should return empty DataFrames with valid columns.""" + + class DummySam: + header = {"PG": [{"ID": "yara", "CL": ""}]} + nreferences = 3 + references = ["HLA:A*01:01", "HLA:B*07:02", "HLA:C*07:02"] + + def __iter__(self): + return iter(()) + + class DummyPysam: + @staticmethod + def AlignmentFile(_samfile, _mode): + return DummySam() + + monkeypatch.setattr(readers, "PYSAM_AVAILABLE", True) + monkeypatch.setattr(readers, "pysam", DummyPysam(), raising=False) + + pos_df, details_df = readers.pysam_to_dataframe("empty.bam") + + assert pos_df.empty + assert list(pos_df.columns) == ["HLA:A*01:01", "HLA:B*07:02", "HLA:C*07:02"] + assert details_df.empty + assert list(details_df.columns) == ["mismatches", "read_length"] diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py index 934cd7e..6adbe32 100644 --- a/tests/test_pipeline.py +++ b/tests/test_pipeline.py @@ -12,6 +12,7 @@ from optitype.io.data import load_reference_data from optitype.pipeline import ( PipelineConfig, + LowCoverageError, run_pipeline, load_config, load_mapper_config, @@ -174,6 +175,33 @@ def test_invalid_enumerate(self, tmp_path): with pytest.raises(ValueError, match="enumerate"): run_pipeline([str(f)], "dna", str(tmp_path), enumerate_count=0) + def test_low_coverage_writes_failure_tsv(self, tmp_path): + out_prefix = "lowcov" + out_csv = tmp_path / f"{out_prefix}_result.tsv" + + with patch( + "optitype.pipeline.pysam_to_dataframe", + return_value=( + pd.DataFrame(columns=["HLA:A*01:01"]), + pd.DataFrame(columns=["mismatches", "read_length"]), + ), + ): + with pytest.raises(LowCoverageError, match="low coverage"): + run_pipeline( + input_files=["sample.bam"], + seq_type="dna", + output_dir=str(tmp_path), + prefix=out_prefix, + config=PipelineConfig(), + verbose=False, + ) + + assert out_csv.exists() + failed = pd.read_csv(out_csv, sep="\t") + assert failed.loc[0, "Status"] == "FAILED" + assert "low coverage" in failed.loc[0, "Message"] + assert int(failed.loc[0, "Reads"]) == 0 + class TestMapperConfig: """Tests for mapper config loading.""" From 5f828afacf99a0e4a46a26e008713c959a3ea0a0 Mon Sep 17 00:00:00 2001 From: MikeWLloyd Date: Tue, 11 Aug 2026 09:29:42 -0400 Subject: [PATCH 2/3] sys.exit 3 for alt catching --- src/optitype/cli.py | 5 ++++- tests/test_cli.py | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/optitype/cli.py b/src/optitype/cli.py index 973835c..049efba 100644 --- a/src/optitype/cli.py +++ b/src/optitype/cli.py @@ -173,7 +173,7 @@ def run( # With custom settings optitype run -i reads.fq --dna -o results/ --solver cbc --threads 8 """ - from optitype.pipeline import PipelineConfig, load_config, load_mapper_config, run_pipeline + from optitype.pipeline import LowCoverageError, PipelineConfig, load_config, load_mapper_config, run_pipeline # Resolve razers3 path (--razers3 option or search PATH) if razers3 is None: @@ -237,6 +237,9 @@ def run( click.echo(f"Results written to: {result.output_csv}") click.echo(f"Coverage plot: {result.output_plot}") + except LowCoverageError as e: + click.echo(str(e), err=True) + raise SystemExit(3) from None except Exception as e: raise click.ClickException(str(e)) from None diff --git a/tests/test_cli.py b/tests/test_cli.py index 08df001..a32cdfc 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -98,6 +98,6 @@ def test_cli_run_low_coverage_reports_reason_and_result_path(tmp_path): ["run", "-i", str(fq), "--dna", "-o", str(outdir), "--prefix", prefix], ) - assert result.exit_code != 0 + assert result.exit_code == 3 assert "low coverage" in result.output.lower() assert str(out_csv) in result.output From 918295e86840493ced8330a0012153b300a2dabe Mon Sep 17 00:00:00 2001 From: MikeWLloyd Date: Tue, 11 Aug 2026 10:14:55 -0400 Subject: [PATCH 3/3] change from exit to warn to allow nextflow to capture output on 'fail' --- src/optitype/cli.py | 4 ++-- tests/test_cli.py | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/optitype/cli.py b/src/optitype/cli.py index 049efba..9a11b06 100644 --- a/src/optitype/cli.py +++ b/src/optitype/cli.py @@ -238,8 +238,8 @@ def run( click.echo(f"Coverage plot: {result.output_plot}") except LowCoverageError as e: - click.echo(str(e), err=True) - raise SystemExit(3) from None + click.secho(f"WARNING: {e}", fg="yellow", err=True) + return except Exception as e: raise click.ClickException(str(e)) from None diff --git a/tests/test_cli.py b/tests/test_cli.py index a32cdfc..dd8085c 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -98,6 +98,7 @@ def test_cli_run_low_coverage_reports_reason_and_result_path(tmp_path): ["run", "-i", str(fq), "--dna", "-o", str(outdir), "--prefix", prefix], ) - assert result.exit_code == 3 + assert result.exit_code == 0 + assert "warning" in result.output.lower() assert "low coverage" in result.output.lower() assert str(out_csv) in result.output