feat: Collection of Minor Changes after review

.github/workflows/python_package.yml

@@ -5,7 +5,7 @@ on:
   workflow_call:
 
 env:
-  UV_VERSION: 0.9.5
+  UV_VERSION: 0.11.28
 
 jobs:
   Tests:

.github/workflows/workflow_linting.yml

@@ -4,7 +4,7 @@ on:
   push:
 
 env:
-  PIXI_VERSION: v0.65.0
+  PIXI_VERSION: v0.72.2
 
 jobs:
   Tests:
@@ -21,4 +21,4 @@ jobs:
           cache: true
 
       - name: Run linting
-        run: pixi run lint --check
+        run: pixi run --locked lint --check

README.md

@@ -30,7 +30,7 @@ To install the GCS connector for Hail/Spark, run
 pixi run setup-gcs
 ```
 
-Log in before running any Hail-dependent tools:
+Log in before running any Hail-dependent tools, requires [the gcloud cli be installed](https://docs.cloud.google.com/sdk/docs/install-sdk):
 
 ```bash
 gcloud auth application-default login
@@ -51,7 +51,7 @@ To install both connectors at once, run `pixi run setup-cloud`.
 
 ### Running the workflow
 
-The workflow does not bundle a default `configfile:` — pass one explicitly with
+The workflow does not set a default `configfile:` — pass one explicitly with
 `--configfile`. Two ready-made configs are provided under `workflows/config/`:
 
 - `config_gcs.yml` — reads all cloud inputs from GCS (`gs://gcp-public-data--gnomad/`,

divref/divref/alias.py

@@ -1,4 +1,2 @@
-from typing import TypeAlias
-
-HailPath: TypeAlias = str
+type HailPath = str
 """Type alias for filesystem paths accepted by Hail: local, S3 (s3a://), GCS (gs://), or HDFS (hdfs://)."""  # noqa: E501

divref/divref/defaults.py

@@ -22,8 +22,11 @@
 )
 """HGDP+1KG sample metadata (AWS S3 alternative)."""
 
-POPULATIONS: list[str] = ["afr", "amr", "eas", "sas", "nfe"]
+# The five continental genetic-ancestry groups DivRef selects from the gnomAD HGDP+1KG legend
+# (a subset of gnomAD's full ancestry set, which also includes e.g. fin, asj, ami, mid, oth).
+# Update this list only to track a deliberate change to the DivRef population definition.
+POPULATIONS: Final[list[str]] = ["afr", "amr", "eas", "sas", "nfe"]
 """Default HGDP+1KG populations."""
 
-REFERENCE_GENOME: str = "GRCh38"
+REFERENCE_GENOME: Final[str] = "GRCh38"
 """Default reference genome assembly."""

divref/divref/duckdb_index.py

@@ -30,6 +30,24 @@ class TablePair(Metric):
     sites_table_path: Path
 
 
+def sequences_table_exists(conn: duckdb.DuckDBPyConnection) -> bool:
+    """
+    Return whether the `sequences` table exists in the connected DuckDB index.
+
+    Args:
+        conn: Open DuckDB connection to a DivRef index.
+
+    Returns:
+        True if a `sequences` table is present.
+    """
+    return (
+        conn.execute(
+            "SELECT 1 FROM information_schema.tables WHERE table_name = 'sequences'"
+        ).fetchone()
+        is not None
+    )
+
+
 def read_pops_legend(table_path: Path) -> list[str]:
     """
     Read a Hail table's population legend from its globals.

divref/divref/hail.py

@@ -18,7 +18,7 @@ def _export_gcs_credentials(gcs_credentials_path: Path | None) -> None:
     """
     if gcs_credentials_path is None:
         raise ValueError("gcs_credentials_path is required when use_s3 is False.")
-    if not gcs_credentials_path.exists():
+    if not gcs_credentials_path.is_file():
         raise FileNotFoundError(
             f"GCS credentials file not found at {gcs_credentials_path}. Run "
             "`gcloud auth application-default login` or pass a valid --gcs-credentials-path."
@@ -83,19 +83,19 @@ def hail_init(
     if not use_s3:
         _export_gcs_credentials(gcs_credentials_path)
 
-    jars_dir = os.path.join(pyspark.__path__[0], "jars")
+    jars_dir = Path(pyspark.__path__[0]) / "jars"
     cloud_jars: list[str] = []
     spark_conf: dict[str, str] = {}
 
     if use_s3:
-        hadoop_aws_jar = os.path.join(jars_dir, "hadoop-aws.jar")
-        aws_sdk_bundle_jar = os.path.join(jars_dir, "aws-java-sdk-bundle.jar")
+        hadoop_aws_jar = jars_dir / "hadoop-aws.jar"
+        aws_sdk_bundle_jar = jars_dir / "aws-java-sdk-bundle.jar"
         for jar in (hadoop_aws_jar, aws_sdk_bundle_jar):
-            if not os.path.exists(jar):
+            if not jar.is_file():
                 raise FileNotFoundError(
                     f"S3 connector JAR not found at {jar}. Run 'pixi run setup-s3' to download it."
                 )
-        cloud_jars.extend([hadoop_aws_jar, aws_sdk_bundle_jar])
+        cloud_jars.extend([str(hadoop_aws_jar), str(aws_sdk_bundle_jar)])
         spark_conf.update({
             "spark.hadoop.fs.s3a.impl": "org.apache.hadoop.fs.s3a.S3AFileSystem",
             # All workflow inputs live on public Open Data buckets that allow anonymous
@@ -112,13 +112,13 @@ def hail_init(
             "spark.hadoop.fs.s3a.threads.max": "64",
         })
     else:
-        gcs_jar = os.path.join(jars_dir, "gcs-connector.jar")
-        if not os.path.exists(gcs_jar):
+        gcs_jar = jars_dir / "gcs-connector.jar"
+        if not gcs_jar.is_file():
             raise FileNotFoundError(
                 f"GCS connector JAR not found at {gcs_jar}. "
                 "Run 'pixi run setup-gcs' to download it."
             )
-        cloud_jars.append(gcs_jar)
+        cloud_jars.append(str(gcs_jar))
         spark_conf.update({
             "spark.hadoop.fs.gs.impl": "com.google.cloud.hadoop.fs.gcs.GoogleHadoopFileSystem",
             "spark.hadoop.fs.AbstractFileSystem.gs.impl": "com.google.cloud.hadoop.fs.gcs.GoogleHadoopFS",  # noqa: E501

divref/divref/haplotype.py

@@ -1,28 +1,9 @@
 """Shared utilities for Hail-based DivRef pipeline tools."""
 
-from typing import Hashable
-from typing import TypeVar
-
 import hail as hl
 
 from divref import defaults
 
-_V = TypeVar("_V", bound=Hashable)
-"""Type variable for hashable dictionary values used in to_hashable_items."""
-
-
-def to_hashable_items(d: dict[str, _V]) -> tuple[tuple[str, _V], ...]:
-    """
-    Convert a dictionary to a sorted tuple of items for use as a hashable key.
-
-    Args:
-        d: Dictionary with hashable values to convert.
-
-    Returns:
-        Sorted tuple of (key, value) pairs.
-    """
-    return tuple(sorted(d.items()))
-
 
 def _max_reference_end(variants: hl.Expression) -> hl.Expression:
     """
@@ -83,7 +64,13 @@ def get_haplo_sequence(
         raise ValueError(
             "get_haplo_sequence requires at least one variant; received an empty sequence"
         )
-    sorted_variants = hl.sorted(variants, key=lambda x: (x.locus.position, hl.len(x.alleles[0])))
+    # Tertiary tiebreak on alt length so a substitution (len(alt) == 1) composes before a pure
+    # insertion (len(alt) > 1) at the same position. A remaining tie (same position, ref length,
+    # and alt length) falls back to input order, but that requires two alleles at one site, which
+    # cannot occur on a single phased haplotype.
+    sorted_variants = hl.sorted(
+        variants, key=lambda x: (x.locus.position, hl.len(x.alleles[0]), hl.len(x.alleles[1]))
+    )
     min_pos = sorted_variants[0].locus.position
     max_ref_end = _max_reference_end(sorted_variants)
     full_context = hl.get_sequence(
@@ -168,7 +155,9 @@ def haplo_coordinates(
     """
     # Same sort key as get_haplo_sequence so the two stay in lockstep; only the minimum position
     # (`[0].locus.position`) is read here, and `_max_reference_end` is order-independent.
-    sorted_variants = hl.sorted(variants, key=lambda x: (x.locus.position, hl.len(x.alleles[0])))
+    sorted_variants = hl.sorted(
+        variants, key=lambda x: (x.locus.position, hl.len(x.alleles[0]), hl.len(x.alleles[1]))
+    )
     min_variant = sorted_variants[0]
     return hl.struct(
         start=min_variant.locus.position - 1 - window_size,

divref/divref/haplotype_compat.py

@@ -12,7 +12,9 @@
 # A component variant: (contig, pos, ref, alt). pos is 1-based.
 Variant = tuple[str, int, str, str]
 
-# Reason labels in precedence order (first match wins in classify_pair) and stable output order.
+# The complete set of incompatibility reason labels classify_pair can return. This tuple is
+# documentation only: classify_pair selects a reason via its own if/elif order (not this tuple),
+# and compatibility_flag emits flags alphabetically (not in this order).
 # A SNP co-located with an indel is NOT listed: the SNP substitutes the shared base and the indel
 # acts on/after it, so the overlap composes to a well-defined haplotype (see classify_pair). An
 # insertion anchored inside a deletion IS a conflict (`insertion_in_deletion`): its anchor base is
@@ -69,11 +71,21 @@ def parse_variants_string(s: str) -> list[Variant]:
 
     Returns:
         Variants sorted by ascending position, then by reference-allele length.
+
+    Raises:
+        ValueError: If a token is not `contig:position:ref:alt` or its position is not an integer.
     """
     out: list[Variant] = []
     for token in s.split(","):
-        contig, pos, ref, alt = token.split(":")
-        out.append((contig, int(pos), ref, alt))
+        fields = token.split(":")
+        if len(fields) != 4:
+            raise ValueError(f"expected 'contig:pos:ref:alt', got {token!r} in {s!r}")
+        contig, pos, ref, alt = fields
+        try:
+            position = int(pos)
+        except ValueError as e:
+            raise ValueError(f"non-integer position in variant token {token!r} in {s!r}") from e
+        out.append((contig, position, ref, alt))
     out.sort(key=lambda v: (v[1], len(v[2])))
     return out
 
@@ -126,7 +138,7 @@ def _same_position_reason(v1: Variant, v2: Variant, k1: str, k2: str) -> str:
 
 def classify_pair(v1: Variant, v2: Variant) -> str | None:
     """
-    Classify an adjacent variant pair into an incompatibility reason.
+    Classify a variant pair (`v1` at or before `v2` by position) into an incompatibility reason.
 
     Returns `None` when the pair is compatible: either it does not overlap (`variant_distance >= 0`)
     or it is a SNP co-located with an indel, which composes to a well-defined haplotype (the SNP

divref/divref/main.py

@@ -46,7 +46,7 @@ def run() -> None:
     """Set up logging, then hand over to defopt for running command line tools."""
     setup_logging()
     logger = logging.getLogger("divref")
-    logger.info("Executing: " + " ".join(sys.argv))
+    logger.info("Executing: %s", " ".join(sys.argv))
     defopt.run(
         funcs=_tools,
         argv=sys.argv[1:],

divref/divref/tools/append_contig_to_duckdb_index.py

@@ -17,6 +17,7 @@
 from divref import defaults
 from divref.duckdb_index import TablePair
 from divref.duckdb_index import read_and_validate_pops_legends
+from divref.duckdb_index import sequences_table_exists
 from divref.haplotype import get_haplo_sequence
 from divref.haplotype import haplo_coordinates
 from divref.haplotype_compat import compatibility_flag
@@ -65,10 +66,7 @@ def read_window_size(conn: duckdb.DuckDBPyConnection) -> int:
 
 def sequences_row_count(conn: duckdb.DuckDBPyConnection) -> int:
     """Current number of rows in `sequences`, or 0 if the table does not exist yet."""
-    exists = conn.execute(
-        "SELECT 1 FROM information_schema.tables WHERE table_name = 'sequences'"
-    ).fetchone()
-    if exists is None:
+    if not sequences_table_exists(conn):
         return 0
     row = conn.execute("SELECT COUNT(*) FROM sequences").fetchone()
     if row is None:
@@ -78,10 +76,7 @@ def sequences_row_count(conn: duckdb.DuckDBPyConnection) -> int:
 
 def _contig_already_appended(conn: duckdb.DuckDBPyConnection, contig: str) -> bool:
     """Whether the `sequences` table already holds any rows for `contig` (False if no table yet)."""
-    exists = conn.execute(
-        "SELECT 1 FROM information_schema.tables WHERE table_name = 'sequences'"
-    ).fetchone()
-    if exists is None:
+    if not sequences_table_exists(conn):
         return False
     row = conn.execute("SELECT 1 FROM sequences WHERE contig = ? LIMIT 1", [contig]).fetchone()
     return row is not None
@@ -279,6 +274,7 @@ def build_contig_sequences_table(
     seq_ht = seq_ht.annotate(
         sequence=get_haplo_sequence(window_size, seq_ht.variants),
         contig=seq_ht.variants[0].locus.contig,
+        # `start`/`end` are 0-based half-open (from `haplo_coordinates`; variant loci are 1-based).
         start=coords.start,
         end=coords.end,
     )
@@ -456,7 +452,10 @@ def iter_dataframe_chunks(
             yield df
 
 
-@dataclass(frozen=True, kw_only=True)
+# eq=False so the frozen dataclass keeps a (default, identity-based) __hash__ despite its list
+# fields; without it, frozen+eq would synthesize a __hash__ that raises TypeError on unhashable
+# lists. The struct is only ever field-accessed, never compared or hashed by value.
+@dataclass(frozen=True, kw_only=True, eq=False)
 class _RemapArrays:
     """
     The four pop-legend remap arrays used to build a contig's sequences table.
@@ -579,10 +578,7 @@ def _stream_tsv_into_sequences(
     Returns:
         The number of rows appended for this contig.
     """
-    table_exists = conn.execute(
-        "SELECT 1 FROM information_schema.tables WHERE table_name = 'sequences'"
-    ).fetchone()
-    if table_exists is None:
+    if not sequences_table_exists(conn):
         # `haplotype_filter` is appended last in both the schema-defining empty frame and every
         # inserted chunk, so column order stays consistent across the CREATE and the INSERTs.
         empty_df = _add_compatibility_flag(
@@ -600,7 +596,9 @@ def _stream_tsv_into_sequences(
     ):
         chunk_df = _add_compatibility_flag(chunk)
         conn.register("chunk_df", chunk_df)
-        conn.execute("INSERT INTO sequences SELECT * FROM chunk_df")
+        # BY NAME matches columns by name, not position, so a future change to the export/scan
+        # column order surfaces as a bind error rather than silently misaligning columns.
+        conn.execute("INSERT INTO sequences BY NAME SELECT * FROM chunk_df")
         conn.unregister("chunk_df")
         appended_rows += chunk_df.height
     return appended_rows
@@ -728,7 +726,8 @@ def append_contig_to_duckdb_index(
             f"Spark executor memory must be at least 1GB. Saw {spark_executor_memory_gb}GB."
         )
 
-    out_duckdb_file: Path = Path(f"{str(output_base)}.haplotypes_gnomad_merge.index.duckdb")
+    out_duckdb_file: Path = Path(f"{output_base}.haplotypes_gnomad_merge.index.duckdb")
+    # readable (not writable): enforces the index file already exists, i.e. init_duckdb_index ran.
     assert_path_is_readable(out_duckdb_file)
 
     table_pairs: list[TablePair] = list(TablePair.read(in_table_pairs_tsv))