From dd48d1ed2d114a94e76d855aa2c258d4ba4d978b Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Fri, 24 Jul 2026 17:38:17 +0000 Subject: [PATCH] v1.1.0 --- changelog.md | 7 +++ zipstrain/pyproject.toml | 2 +- zipstrain/src/zipstrain/cli.py | 9 +++- zipstrain/src/zipstrain/compare.py | 57 ++++++++++++++++++++---- zipstrain/src/zipstrain/database.py | 1 + zipstrain/src/zipstrain/task_manager.py | 58 +++++++++++++++++++++++-- 6 files changed, 120 insertions(+), 14 deletions(-) diff --git a/changelog.md b/changelog.md index b2e722a..ff090ba 100644 --- a/changelog.md +++ b/changelog.md @@ -2,6 +2,13 @@ Entries are brief by design and describe changes relative to the previous released version. +## 1.1.0 + +Compared with `1.0.1`: + +- Compare: popANI now applies a **minor-allele frequency floor** before counting shared alleles (`--min-snp-freq`, default `0.01`). A base is a "present allele" only if its frequency is ≥ the floor or it is the consensus base, so low-frequency sequencing errors no longer make two deeply-sequenced samples share an allele at (nearly) every position. This fixes popANI being inflated toward 100% at high coverage. Set `--min-snp-freq 0` for the previous raw-count behaviour. Applied consistently in the polars and duckdb engines; `GenomeComparisonConfig` gains a `min_snp_freq` field. (Distinct from the profile-side `--min-freq`: this floor is applied at compare time, so it also corrects popANI on already-written profiles.) +- Profiling: fixed a crash when the output directory is on a CIFS/SMB mount (e.g. PetaLibrary `/pl`). The post-run reorganize step used `shutil.move`, which falls back to permission-preserving copies (`copystat`/`chmod`) that raise `Operation not permitted` on such mounts, aborting `zipstrain profile` *after* all profiles were written (leaving them stranded in `batch_*/`). It now uses a CIFS-safe move (atomic rename, else metadata-free copy + symlink recreation). + ## 1.0.1 Compared with `1.0.0`: diff --git a/zipstrain/pyproject.toml b/zipstrain/pyproject.toml index b541d90..0e60666 100644 --- a/zipstrain/pyproject.toml +++ b/zipstrain/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "zipstrain" -version = "1.0.1" +version = "1.1.0" description = "" authors = [ {name = "ParsaGhadermazi",email = "54489047+ParsaGhadermazi@users.noreply.github.com"} diff --git a/zipstrain/src/zipstrain/cli.py b/zipstrain/src/zipstrain/cli.py index 60f6dd0..57fecf1 100644 --- a/zipstrain/src/zipstrain/cli.py +++ b/zipstrain/src/zipstrain/cli.py @@ -1368,7 +1368,8 @@ def get_gene_range_table(gene_file, output_file): @click.option('--duckdb-memory-limit', default=None, help="DuckDB memory limit (e.g., 2GB, 1024MB).") @click.option('--duckdb-temp-directory', default=None, help="Directory DuckDB can use for spill files.") @click.option('--duckdb-threads', type=int, default=None, help="Number of DuckDB worker threads.") -def single_compare_genome(profile_location_1, profile_location_2, stb_file, min_cov, min_gene_compare_len, output_file, genome, ani_method, calculate, engine, duckdb_memory_limit, duckdb_temp_directory, duckdb_threads): +@click.option('--min-snp-freq', type=float, default=cp.DEFAULT_MIN_SNP_FREQ, show_default=True, help="Minor-allele frequency floor for popANI: a base counts as a present allele only if its frequency is >= this value (or it is the consensus base). Removes low-frequency sequencing errors so popANI is not inflated at high coverage. Set 0 to disable (legacy raw-count behaviour).") +def single_compare_genome(profile_location_1, profile_location_2, stb_file, min_cov, min_gene_compare_len, output_file, genome, ani_method, calculate, engine, duckdb_memory_limit, duckdb_temp_directory, duckdb_threads, min_snp_freq): """ Compare two profile parquets and calculate genome-level comparison statistics. @@ -1424,6 +1425,7 @@ def single_compare_genome(profile_location_1, profile_location_2, stb_file, min_ memory_limit=duckdb_memory_limit, temp_directory=duckdb_temp_directory, threads=duckdb_threads, + min_snp_freq=min_snp_freq, ) ut.rewrite_parquet_with_metadata(output_file, compare_metadata) return @@ -1441,6 +1443,7 @@ def single_compare_genome(profile_location_1, profile_location_2, stb_file, min_ engine="polars", stb_file=stb_file, calculate=calculations, + min_snp_freq=min_snp_freq, ).with_columns( sample_1=pl.lit(profile_1_name), sample_2=pl.lit(profile_2_name), @@ -1986,6 +1989,7 @@ def _on_step(message: str) -> None: @click.option("--compare-genes", is_flag=True, default=False, show_default=True, help="Compare genes instead of genomes.") @click.option("--scope", default=None, help="Comparison scope. Defaults to 'all' for genomes and 'all:all' for genes.") @click.option("--min-cov", default=5, show_default=True, help="Minimum coverage to consider a position.") +@click.option("--min-snp-freq", type=float, default=cp.DEFAULT_MIN_SNP_FREQ, show_default=True, help="Minor-allele frequency floor for popANI: a base counts as a present allele only if its frequency is >= this value (or it is the consensus base). Removes low-frequency sequencing errors so popANI is not inflated at high coverage. Set 0 to disable (legacy raw-count behaviour).") @click.option("--min-gene-compare-len", default=100, show_default=True, help="Minimum gene length to consider for comparison.") @click.option("--stb-file", default=None, help="Scaffold-to-genome mapping file. Required for --method matrix.") @click.option("--comp-db-file", default=None, help="Optional existing comparison parquet to resume/extend (standard method). Auto-detected from --run-dir if omitted.") @@ -2008,7 +2012,7 @@ def _on_step(message: str) -> None: @click.option("--container-address", default=None, help="Optional container image/address override. Defaults to the current ZipStrain version tag for docker/apptainer.") @click.option("--no-csv", is_flag=True, default=False, show_default=True, help="Do not write a companion .csv next to the comparison parquet.") @click.option("--force-csv", is_flag=True, default=False, show_default=True, help="Write the companion .csv even when the estimated size exceeds 100 MB.") -def compare(profile_db, run_dir, method, compare_genes, scope, min_cov, min_gene_compare_len, stb_file, comp_db_file, allow_mismatch, ani_method, engine, calculate, bed_file, gene_range_table, backend, memory_limit_gb, duckdb_memory_limit, duckdb_threads, max_concurrent_batches, poll_interval, task_per_batch, execution_mode, slurm_config, container_engine, container_address, no_csv, force_csv): +def compare(profile_db, run_dir, method, compare_genes, scope, min_cov, min_snp_freq, min_gene_compare_len, stb_file, comp_db_file, allow_mismatch, ani_method, engine, calculate, bed_file, gene_range_table, backend, memory_limit_gb, duckdb_memory_limit, duckdb_threads, max_concurrent_batches, poll_interval, task_per_batch, execution_mode, slurm_config, container_engine, container_address, no_csv, force_csv): """ Compare profiled samples at the genome level (default) or gene level (--compare-genes). @@ -2097,6 +2101,7 @@ def compare(profile_db, run_dir, method, compare_genes, scope, min_cov, min_gene min_cov=min_cov, min_gene_compare_len=min_gene_compare_len, stb_file_loc=stb_file, + min_snp_freq=min_snp_freq, ), comp_db_loc=comp_db_file, ) diff --git a/zipstrain/src/zipstrain/compare.py b/zipstrain/src/zipstrain/compare.py index e254de7..e2fdc2d 100644 --- a/zipstrain/src/zipstrain/compare.py +++ b/zipstrain/src/zipstrain/compare.py @@ -194,6 +194,29 @@ def _duckdb_configure_connection( conn.execute(f"SET threads={int(threads)}") +DEFAULT_MIN_SNP_FREQ = 0.01 + + +def _apply_snp_floor(lf: "pl.LazyFrame", min_snp_freq: float) -> "pl.LazyFrame": + """Zero minor-allele counts below ``min_snp_freq`` frequency; the consensus base is + always kept. Removes low-frequency sequencing-error alleles so that popANI (which + counts a position as "shared" when the two samples share any present allele) is not + inflated toward 100% at high coverage. ``min_snp_freq <= 0`` disables the floor + (legacy raw-count behaviour).""" + if not min_snp_freq or min_snp_freq <= 0: + return lf + bases = ("A", "T", "C", "G") + cov = pl.col("A") + pl.col("T") + pl.col("C") + pl.col("G") + mx = pl.max_horizontal(*[pl.col(b) for b in bases]) + return lf.with_columns([ + pl.when((pl.col(b) >= min_snp_freq * cov) | (pl.col(b) >= mx)) + .then(pl.col(b)) + .otherwise(0) + .alias(b) + for b in bases + ]) + + def _duckdb_ani_expression(ani_method: str) -> str: if ani_method == "popani": # Cast at expression time (post-join) to avoid upfront casting in filtered scans. @@ -249,11 +272,23 @@ def _duckdb_shared_query( genome_scope: str, ani_method: str, gene_scope: str = "all", + min_snp_freq: float = DEFAULT_MIN_SNP_FREQ, ) -> str: ani_expr = _duckdb_ani_expression(ani_method) con_expr = _duckdb_ani_expression("conani") genome_scope_sql = _duckdb_quote_sql_string(genome_scope) gene_scope_sql = _duckdb_quote_sql_string(gene_scope) + if min_snp_freq and min_snp_freq > 0: + # A base is a "present allele" only if its frequency is >= min_snp_freq or it is + # the consensus (max) base; low-frequency sequencing-error alleles are zeroed so + # popANI is not inflated at high coverage. WHERE still filters on raw coverage. + base_cols = ",\n ".join( + f"CASE WHEN CAST({b} AS DOUBLE) >= {float(min_snp_freq)} * (A + T + C + G) " + f"OR {b} >= GREATEST(A, T, C, G) THEN {b} ELSE 0 END AS {b}" + for b in ("A", "T", "C", "G") + ) + else: + base_cols = "A,\n T,\n C,\n G" return f""" WITH p1f AS ( SELECT @@ -261,10 +296,7 @@ def _duckdb_shared_query( pos, gene, genome, - A, - T, - C, - G + {base_cols} FROM {p1_source} p1 WHERE (A + T + C + G) >= {min_cov} AND ('{genome_scope_sql}' = 'all' OR genome = '{genome_scope_sql}') @@ -276,10 +308,7 @@ def _duckdb_shared_query( pos, gene, genome, - A, - T, - C, - G + {base_cols} FROM {p2_source} p2 WHERE (A + T + C + G) >= {min_cov} AND ('{genome_scope_sql}' = 'all' OR genome = '{genome_scope_sql}') @@ -558,6 +587,7 @@ def _shared_loci_polars( genome_scope: str = "all", gene_scope: str = "all", ani_method: str = "popani", + min_snp_freq: float = DEFAULT_MIN_SNP_FREQ, ) -> pl.LazyFrame: ani_expr = getattr(PolarsANIExpressions(), ani_method)() con_expr = PolarsANIExpressions().conani() @@ -566,6 +596,8 @@ def _shared_loci_polars( mpile2=mpile2, min_cov=min_cov, ) + p1 = _apply_snp_floor(p1, min_snp_freq) + p2 = _apply_snp_floor(p2, min_snp_freq) if _profile_is_coordinate_sorted(mpile1) and _profile_is_coordinate_sorted(mpile2): p1 = p1.set_sorted(["chrom", "pos"]) p2 = p2.set_sorted(["chrom", "pos"]) @@ -719,6 +751,7 @@ def duckdb_compare_genomes_to_parquet( memory_limit: Optional[str] = None, temp_directory: Optional[Union[str, Path]] = None, threads: Optional[int] = None, + min_snp_freq: float = DEFAULT_MIN_SNP_FREQ, ) -> None: """Run genome comparison in DuckDB and write final output directly to parquet. @@ -740,6 +773,7 @@ def duckdb_compare_genomes_to_parquet( min_cov=min_cov, genome_scope=genome_scope, ani_method=ani_method, + min_snp_freq=min_snp_freq, ) query = _duckdb_genome_compare_query( shared_query=shared_query, @@ -767,6 +801,7 @@ def duckdb_compare_genomes( memory_limit: Optional[str] = None, temp_directory: Optional[Union[str, Path]] = None, threads: Optional[int] = None, + min_snp_freq: float = DEFAULT_MIN_SNP_FREQ, ) -> pl.LazyFrame: """Run genome comparison in DuckDB and return selected metrics as a LazyFrame.""" con = duckdb.connect() @@ -785,6 +820,7 @@ def duckdb_compare_genomes( min_cov=min_cov, genome_scope=genome_scope, ani_method=ani_method, + min_snp_freq=min_snp_freq, ) query = _duckdb_genome_compare_query( shared_query=shared_query, @@ -989,6 +1025,7 @@ def compare_genomes_polars( ani_method: str = "popani", stb_file: Optional[Union[str, Path]] = None, calculate: Optional[Union[str, Iterable[str]]] = None, + min_snp_freq: float = DEFAULT_MIN_SNP_FREQ, ) -> pl.LazyFrame: """Compare two profiles fully in Polars and return genome-level statistics.""" calculations = parse_genome_calculations(calculate) @@ -998,6 +1035,7 @@ def compare_genomes_polars( min_cov=min_cov, genome_scope=genome_scope, ani_method=ani_method, + min_snp_freq=min_snp_freq, ) genome_comp_parts: list[pl.LazyFrame] = [] if "ani" in calculations: @@ -1129,6 +1167,7 @@ def compare_genomes( engine: Literal["polars", "duckdb"] = "polars", stb_file: Optional[Union[str, Path]] = None, calculate: Optional[Union[str, Iterable[str]]] = None, + min_snp_freq: float = DEFAULT_MIN_SNP_FREQ, ) -> pl.LazyFrame: """Compare two profiles with selectable execution engine.""" calculations = parse_genome_calculations(calculate) @@ -1142,6 +1181,7 @@ def compare_genomes( ani_method=ani_method, stb_file=stb_file, calculate=calculations, + min_snp_freq=min_snp_freq, ) if engine == "duckdb": return duckdb_compare_genomes( @@ -1156,6 +1196,7 @@ def compare_genomes( memory_limit=duckdb_memory_limit, temp_directory=duckdb_temp_directory, threads=duckdb_threads, + min_snp_freq=min_snp_freq, ) raise ValueError(f"Unsupported engine: {engine}") diff --git a/zipstrain/src/zipstrain/database.py b/zipstrain/src/zipstrain/database.py index 60c032c..bcf0576 100644 --- a/zipstrain/src/zipstrain/database.py +++ b/zipstrain/src/zipstrain/database.py @@ -224,6 +224,7 @@ class GenomeComparisonConfig(BaseModel): min_cov: int =Field(description="Minimum coverage a base on the reference fasta that must have in order to be compared.") min_gene_compare_len: int=Field(description="Minimum length of a gene that needs to be covered at min_cov to be considered for gene similarity calculations") stb_file_loc:str | None = Field(default=None, description="Optional location of the scaffold to bin file.") + min_snp_freq: float = Field(default=0.01, description="Minor-allele frequency floor for popANI presence: a base counts as a present allele only if its frequency is >= this value (or it is the consensus base). Removes low-frequency sequencing errors so popANI is not inflated at high coverage. Set 0 to disable.") def is_compatible(self, other: GenomeComparisonConfig) -> bool: """ diff --git a/zipstrain/src/zipstrain/task_manager.py b/zipstrain/src/zipstrain/task_manager.py index 05c4d53..370b01f 100644 --- a/zipstrain/src/zipstrain/task_manager.py +++ b/zipstrain/src/zipstrain/task_manager.py @@ -56,6 +56,7 @@ from rich.columns import Columns import polars as pl import psutil +import os import shutil import signal from datetime import datetime, timezone @@ -684,6 +685,7 @@ async def generate_tasks(self) -> list[Task]: "profile_2_file": FileInput(row["profile_location_2"]), "stb-file-arg": StringInput(stb_file_arg), "min_cov": IntInput(self.comp_config.min_cov), + "min-snp-freq-arg": StringInput(f"--min-snp-freq {self.comp_config.min_snp_freq}"), "min-gene-compare-len": IntInput(self.comp_config.min_gene_compare_len), "ani-method-arg": StringInput(ani_method_arg), "calculate-arg": StringInput(calculate_arg), @@ -1755,6 +1757,7 @@ class FastCompareTask(Task): --profile-location-2 \ \ --min-cov \ + \ --min-gene-compare-len \ \ \ @@ -1855,6 +1858,55 @@ class PrepareCompareGenomeRunOutputsSlurmBatch(SlurmBatch): COMPARE_FINAL_BATCH_DIRNAME = "Outputs" +def _cifs_safe_move(src, dst) -> None: + """Move ``src`` to ``dst`` without permission-preserving copies. + + ``shutil.move`` falls back to ``copytree``/``copy2`` when ``os.rename`` + fails, and those call ``copystat``/``chmod`` which raise + ``Operation not permitted`` on CIFS/SMB mounts such as PetaLibrary + (``/pl``). Here we try an atomic rename first, then fall back to a manual + recursive move that copies file *contents only* (``shutil.copyfile``, + no metadata) and recreates symlinks verbatim. + """ + src = os.fspath(src) + dst = os.fspath(dst) + if os.path.islink(src): + target = os.readlink(src) + if os.path.lexists(dst): + if os.path.isdir(dst) and not os.path.islink(dst): + shutil.rmtree(dst, ignore_errors=True) + else: + try: + os.remove(dst) + except OSError: + pass + os.symlink(target, dst) + try: + os.remove(src) + except OSError: + pass + return + try: + os.replace(src, dst) # fast atomic path on POSIX filesystems + return + except OSError: + pass + if os.path.isdir(src): + os.makedirs(dst, exist_ok=True) + for name in os.listdir(src): + _cifs_safe_move(os.path.join(src, name), os.path.join(dst, name)) + try: + os.rmdir(src) + except OSError: + shutil.rmtree(src, ignore_errors=True) + else: + shutil.copyfile(src, dst) # data only; no chmod/copystat (CIFS-safe) + try: + os.remove(src) + except OSError: + pass + + def _move_into(destination_dir: pathlib.Path, source: pathlib.Path) -> None: """Move ``source`` into ``destination_dir``, prefixing on name collision. @@ -1865,7 +1917,7 @@ def _move_into(destination_dir: pathlib.Path, source: pathlib.Path) -> None: target = destination_dir / source.name if target.exists(): target = destination_dir / f"{source.parent.name}_{source.name}" - shutil.move(str(source), str(target)) + _cifs_safe_move(source, target) def _reorganize_profile_run_output(run_dir: pathlib.Path) -> None: @@ -1901,7 +1953,7 @@ def _reorganize_profile_run_output(run_dir: pathlib.Path) -> None: if entry.is_dir(): # A per-sample task directory: lift it to run_dir/ and tidy. sample_dir = run_dir / entry.name - shutil.move(str(entry), str(sample_dir)) + _cifs_safe_move(entry, sample_dir) intermediate_dir = sample_dir / INTERMEDIATE_FILES_DIRNAME for item in sorted(sample_dir.iterdir()): if item.name == INTERMEDIATE_FILES_DIRNAME: @@ -1968,7 +2020,7 @@ def _reorganize_compare_run_output(run_dir: pathlib.Path) -> None: _move_into(log_dir, entry) # Whatever remains (task subdirectories) is intermediate output. if any(batch_dir.iterdir()): - shutil.move(str(batch_dir), str(intermediate_dir / batch_dir.name)) + _cifs_safe_move(batch_dir, intermediate_dir / batch_dir.name) else: try: batch_dir.rmdir()