v1.1.0

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`:

zipstrain/pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "zipstrain"
-version = "1.0.1"
+version = "1.1.0"
 description = ""
 authors = [
     {name = "ParsaGhadermazi",email = "[email protected]"}

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,
             )

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,22 +272,31 @@ 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
         chrom,
         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}")
 

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:
         """

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 <profile_2_file> \
     <stb-file-arg> \
     --min-cov <min_cov> \
+    <min-snp-freq-arg> \
     --min-gene-compare-len <min-gene-compare-len> \
     <ani-method-arg> \
     <calculate-arg> \
@@ -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/<sample> 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()