diff --git a/.gitattributes b/.gitattributes index ef815a9..36bfc38 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,2 +1,3 @@ cpp/**/btllib/** linguist-vendored -tests/**/*.py linguist-vendored \ No newline at end of file +tests/**/*.py linguist-vendored +*.npz binary \ No newline at end of file diff --git a/src/seqwin/_version.py b/src/seqwin/_version.py index 222c11c..b703f5c 100644 --- a/src/seqwin/_version.py +++ b/src/seqwin/_version.py @@ -1 +1 @@ -__version__ = '0.4.0' \ No newline at end of file +__version__ = '0.4.1' \ No newline at end of file diff --git a/src/seqwin/cli.py b/src/seqwin/cli.py index 4f7a19f..59b68d6 100644 --- a/src/seqwin/cli.py +++ b/src/seqwin/cli.py @@ -140,6 +140,9 @@ def main( help='Do not evaluate signature sequences with BLAST.', rich_help_panel='Signature options' ), + no_filter: bool = typer.Option( + False, '--no-filter', show_default=False, hidden=True + ), # blast_neg_only: bool = typer.Option( # False, '--fast-blast', is_flag=True, flag_value=True, show_default=False, # help='Only evaluate (BLAST) against non-target assemblies.' @@ -227,6 +230,7 @@ def main( min_len=min_len, max_len=max_len, run_blast=not no_blast, + no_filter=no_filter, #blast_neg_only=blast_neg_only, seed=seed, n_cpu=n_cpu, diff --git a/src/seqwin/config.py b/src/seqwin/config.py index 1ea5874..dcdf19d 100644 --- a/src/seqwin/config.py +++ b/src/seqwin/config.py @@ -94,9 +94,9 @@ class Config(BaseModel): min_len (int): Minimum length of output signatures. [200] max_len (int | None): Estimated maximum length of output signatures. If None, no explicit limit is applied (capped by `max_nodes_cap`). [None] run_blast (bool): If True, evaluate signature sequences with BLAST. [True] + no_filter (bool): If True, skip filtering k-mers and save the raw k-mer graph. [False] blast_neg_only (bool, deprecated): If True, only include non-target assemblies in the BLAST database. [False] - no_filter (bool): If True, skip filtering k-mers (debug only). [False] penalty_th_cap (float): If `penalty_th` is None, penalty threshold cannot be higher than this value. [0.2] edge_w_th_mul (float): Multiplier for determining the threshold for low-weight edges. [0.3] min_nodes_floor (int): Lowest possible value for `min_nodes` (see `RunState`), regardless of `min_len`. [3] @@ -138,10 +138,10 @@ class Config(BaseModel): min_len: int = 200 max_len: int | None = None run_blast: bool = True + no_filter: bool = False blast_neg_only: bool = False # NOTE: need to fix when this is turned on # Graph filtering options (not included in CLI) - no_filter: bool = False penalty_th_cap: float = 0.2 edge_w_th_mul: float = 0.3 min_nodes_floor: int = 3 @@ -255,6 +255,7 @@ class WorkingDir: config (str): The `Config` instance saved as JSON. ['config.json'] assemblies_dir (str): Directory for downloaded assemblies. ['assemblies'] assemblies_csv (str): The `Assemblies` instance (`assemblies.py`) saved as CSV. ['assemblies.csv'] + graph (str): NumPy arrays representing the k-mer graph (`kmers`, `idx`, `nodes` and `edges`). ['graph.npz'] mash (str): Mash sketch of all assemblies (.msh will be added by Mash). ['sketches'] blast_dir (str): Directory for the BLAST database. ['blastdb'] blast_log (str): Console output of the makeblastdb command, inside `blast_dir`. ['makeblastdb.log'] @@ -266,6 +267,7 @@ class WorkingDir: config: str = 'config.json' assemblies_dir: str = 'assemblies' assemblies_csv: str = 'assemblies.csv' + graph: str = 'graph.npz' mash: str = 'sketches' blast_dir: str = 'blastdb' blast_log: str = 'makeblastdb.log' diff --git a/src/seqwin/core.py b/src/seqwin/core.py index 9a7d8be..cc8867d 100644 --- a/src/seqwin/core.py +++ b/src/seqwin/core.py @@ -6,6 +6,7 @@ Dependencies: ------------- +- numpy - pandas - .assemblies - .kmers @@ -32,6 +33,7 @@ logger = logging.getLogger(__name__) +import numpy as np import pandas as pd from .assemblies import Assemblies, get_assemblies @@ -126,25 +128,35 @@ def run(self) -> None: assemblies = self.assemblies overwrite = config.overwrite + no_filter = config.no_filter working_dir = state.working_dir kmers, jaccard = get_kmers(assemblies, config, state) - if kmers.subgraphs is None: - # if config.no_filter is True (debug only) - markers = None + if no_filter: + graph_path = working_dir / WORKINGDIR.graph + file_to_write(graph_path, overwrite) + np.savez( + graph_path, + allow_pickle=False, + kmers=kmers.kmers, + idx=kmers.idx, + nodes=kmers.nodes, + edges=kmers.edges, + ) + logger.info(f'Filtering is turned off. Raw minimizer graph is saved as {graph_path}') else: markers = get_markers(kmers, assemblies, config, state) - self.kmers = kmers - self.mash = jaccard - self.markers = markers + self.kmers = kmers + self.mash = jaccard + self.markers = markers - # save run instance - results_path = working_dir / WORKINGDIR.results - file_to_write(results_path, overwrite) - results_path.write_bytes(pickle.dumps(self)) - logger.info(f'Run instance (includes all run data) saved as {results_path}') + # save run instance + results_path = working_dir / WORKINGDIR.results + file_to_write(results_path, overwrite) + results_path.write_bytes(pickle.dumps(self)) + logger.info(f'Run instance (includes all run data) saved as {results_path}') def run(config: Config) -> Seqwin: diff --git a/src/seqwin/kmers.py b/src/seqwin/kmers.py index bf1f4aa..1721819 100644 --- a/src/seqwin/kmers.py +++ b/src/seqwin/kmers.py @@ -271,17 +271,7 @@ def get_kmers( kmers = KmerGraph(assemblies, kmerlen, windowsize, n_cpu) if no_filter: - # skip kmers.filter(), debug only - kmers_path = working_dir / 'kmers.npz' - np.savez( - kmers_path, - kmers=kmers.kmers, - idx=kmers.idx, - nodes=kmers.nodes, - edges=kmers.edges, - allow_pickle=False - ) - log_and_raise(SystemExit, f'Filtering is turned off. Raw k-mers are saved to {kmers_path}') + return kmers, None # calculate filter params # 1. calculate penalty threshold diff --git a/tests/smoke/conftest.py b/tests/smoke/conftest.py index d45d98d..9813ea7 100644 --- a/tests/smoke/conftest.py +++ b/tests/smoke/conftest.py @@ -38,3 +38,8 @@ def non_targets_txt(fixtures_dir: Path) -> Path: @pytest.fixture(scope='session') def expected_fasta(fixtures_dir: Path) -> str: return read_text(fixtures_dir / 'expected' / 'signatures.fasta') + + +@pytest.fixture(scope='session') +def expected_graph(fixtures_dir: Path) -> Path: + return fixtures_dir / 'expected' / 'graph.npz' diff --git a/tests/smoke/fixtures/expected/graph.npz b/tests/smoke/fixtures/expected/graph.npz new file mode 100644 index 0000000..49cfa34 Binary files /dev/null and b/tests/smoke/fixtures/expected/graph.npz differ diff --git a/tests/smoke/test_outputs.py b/tests/smoke/test_outputs.py index af26309..c146b61 100644 --- a/tests/smoke/test_outputs.py +++ b/tests/smoke/test_outputs.py @@ -1,5 +1,6 @@ from pathlib import Path +import numpy as np from typer.testing import CliRunner from seqwin import cli @@ -18,6 +19,45 @@ ) +def _assert_graph_matches_expected(actual_path: Path, expected_path: Path) -> None: + actual = np.load(actual_path, allow_pickle=False) + expected = np.load(expected_path, allow_pickle=False) + + actual_keys = set(actual.files) + expected_keys = set(expected.files) + required_keys = {'kmers', 'idx', 'nodes', 'edges'} + + assert actual_keys == expected_keys, f'graph arrays mismatch: actual={sorted(actual_keys)}, expected={sorted(expected_keys)}' + assert actual_keys == required_keys, f'graph arrays must be exactly {sorted(required_keys)}, got {sorted(actual_keys)}' + + for name in sorted(required_keys): + actual_array = actual[name] + expected_array = expected[name] + assert actual_array.dtype == expected_array.dtype, f'{name} dtype mismatch: actual={actual_array.dtype}, expected={expected_array.dtype}' + assert actual_array.shape == expected_array.shape, f'{name} shape mismatch: actual={actual_array.shape}, expected={expected_array.shape}' + + np.testing.assert_array_equal(actual['idx'], expected['idx'], err_msg='idx array values mismatch') + np.testing.assert_array_equal(actual['kmers'], expected['kmers'], err_msg='kmers array values mismatch') + np.testing.assert_array_equal(actual['edges'], expected['edges'], err_msg='edges array values mismatch') + + nodes_dtype = actual['nodes'].dtype + for field_name in nodes_dtype.names or (): + if field_name == 'penalty': + np.testing.assert_allclose( + actual['nodes'][field_name], + expected['nodes'][field_name], + rtol=0, + atol=1e-12, + err_msg='nodes.penalty values mismatch', + ) + else: + np.testing.assert_array_equal( + actual['nodes'][field_name], + expected['nodes'][field_name], + err_msg=f'nodes.{field_name} values mismatch', + ) + + def _run_cli(*args: str) -> Path: result = runner.invoke(cli.app, args) assert result.exit_code == 0, result.stdout @@ -66,3 +106,17 @@ def test_multithreading_matches_expected(tmp_path: Path, targets_txt: Path, non_ ) assert read_text(out_dir / WORKINGDIR.markers_fasta) == expected_fasta + + +def test_graph_matches_expected(tmp_path: Path, targets_txt: Path, non_targets_txt: Path, expected_graph: Path) -> None: + out_dir = _run_cli( + '--tar-paths', str(targets_txt), + '--neg-paths', str(non_targets_txt), + '--prefix', str(tmp_path), + '--threads', '1', + '--title', 'no-filter', + '--no-filter', + *_shared_config, + ) + + _assert_graph_matches_expected(out_dir / WORKINGDIR.graph, expected_graph)