Merge pull request #355 from jkmckenna/0.3.2

0.3.2

Author: jkmckenna

.github/workflows/ci.yml

@@ -2,9 +2,9 @@ name: CI
 
 on:
   push:
-    branches: ["main", "0.3.0"]
+    branches: ["main"]
   pull_request:
-    branches: ["main", "0.3.0"]
+    branches: ["main"]
 
 concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
@@ -49,7 +49,7 @@ jobs:
       - name: Lint with ruff
         run: ruff check --output-format=github .
 
-  smoke:
+  pytest:
     runs-on: ubuntu-latest
     strategy:
       fail-fast: false
@@ -74,9 +74,9 @@ jobs:
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
-          python -m pip install .[dev]
-      - name: Run smoke tests
-        run: pytest -m smoke -q
+          python -m pip install .[dev,torch,plotting]
+      - name: Run pytest
+        run: pytest -m "smoke" -q
 
   docs:
     runs-on: ubuntu-latest

AGENTS.md

@@ -1,9 +1,10 @@
 # AGENTS.md
 
-This file tells coding agents (including OpenAI Codex and Claude Code) how to work in this repo.
+This file tells coding agents (including OpenAI's Codex, Anthropic's Claude Code, and Google's Gemini) how to work in this repo.
 
-Coding agents can only read from AGENTS.md or Claude.md files.
-Agents can not edit AGENTS.md or Claude.md files.
+- For AGENTS.md or CLAUDE.md files:
+  - Agents can read from these files.
+  - Agents can never edit these files.
 
 ## Goals
 - Make minimal, correct changes.
@@ -12,15 +13,17 @@ Agents can not edit AGENTS.md or Claude.md files.
 - Generate production grade, scalable code.
 
 ## Prompt interface
-- When asked about a problem or task, first describe the plan to handle the task.
-- Keep taking prompts until the plan is validated.
+- When asked about a problem or task, first read all files relevent to the task's scope.
+- Describe the problem given the context.
+- Formulate a plan to address the problem within scope.
+- Refine the plan with user input.
 - Implement code after being told to proceed.
 
 ## Repo orientation
 - Read existing patterns before inventing new ones.
 - Don’t refactor broadly unless asked.
-- If you’re unsure about intended behavior, look for tests/docs first.
-- If behavior is not clear after tests/docs, look at the Click commands section in this file.
+- If you’re unsure about intended behavior, look for tests or docs first.
+- If behavior is not clear after reading tests and docs, look at the Click commands section in this file.
 - Ignore all files in any directory named "archived".
 - User defined parameters exist within src/smftools/config.
 - Parameters are herited from default.yaml -> MODALITY.yaml -> user_defined_config.csv

src/smftools/cli/chimeric_adata.py

@@ -14,6 +14,47 @@
 ZERO_HAMMING_DISTANCE_SPANS = "zero_hamming_distance_spans"
 
 
+def _max_positive_span_length(delta_row: "np.ndarray") -> int:
+    """Return the max contiguous run length where delta span values are > 0."""
+    import numpy as np
+
+    values = np.asarray(delta_row)
+    if values.ndim != 1 or values.size == 0:
+        return 0
+
+    positive_mask = values > 0
+    if not np.any(positive_mask):
+        return 0
+
+    transitions = np.diff(positive_mask.astype(np.int8))
+    starts = np.flatnonzero(transitions == 1) + 1
+    ends = np.flatnonzero(transitions == -1) + 1
+
+    if positive_mask[0]:
+        starts = np.r_[0, starts]
+    if positive_mask[-1]:
+        ends = np.r_[ends, positive_mask.size]
+
+    return int(np.max(ends - starts))
+
+
+def _compute_chimeric_by_mod_hamming_distance(
+    delta_layer: "np.ndarray",
+    span_threshold: int,
+) -> "np.ndarray":
+    """Flag reads with any delta-hamming span strictly larger than ``span_threshold``."""
+    import numpy as np
+
+    delta_values = np.asarray(delta_layer)
+    if delta_values.ndim != 2:
+        raise ValueError("delta_layer must be a 2D array with shape (n_obs, n_vars).")
+
+    flags = np.zeros(delta_values.shape[0], dtype=bool)
+    for obs_idx, row in enumerate(delta_values):
+        flags[obs_idx] = _max_positive_span_length(row) > span_threshold
+    return flags
+
+
 def _build_top_segments_obs_tuples(
     read_df: "pd.DataFrame",
     obs_names: "pd.Index",
@@ -295,7 +336,10 @@ def chimeric_adata_core(
             suffix = "_strand_FASTA_base"
             variant_seq1_label = seq1_col[: -len(suffix)] if seq1_col.endswith(suffix) else seq1_col
             variant_seq2_label = seq2_col[: -len(suffix)] if seq2_col.endswith(suffix) else seq2_col
-            logger.info("Detected variant call layer '%s'; will overlay on span clustermaps.", variant_call_layer_name)
+            logger.info(
+                "Detected variant call layer '%s'; will overlay on span clustermaps.",
+                variant_call_layer_name,
+            )
 
     # ============================================================
     # 1) Rolling NN distances + layer clustermaps
@@ -625,7 +669,9 @@ def chimeric_adata_core(
                         _vc = adata[mask].layers[variant_call_layer_name]
                         _vc = _vc.toarray() if hasattr(_vc, "toarray") else np.asarray(_vc)
                         _variant_call_df = pd.DataFrame(
-                            _vc, index=adata[mask].obs_names.astype(str), columns=adata.var_names,
+                            _vc,
+                            index=adata[mask].obs_names.astype(str),
+                            columns=adata.var_names,
                         )
 
                     subset = subset[:, site_mask].copy()
@@ -641,6 +687,9 @@ def chimeric_adata_core(
                             variant_call_data=_variant_call_df,
                             seq1_label=variant_seq1_label,
                             seq2_label=variant_seq2_label,
+                            ref1_marker_color=getattr(cfg, "variant_overlay_seq1_color", "white"),
+                            ref2_marker_color=getattr(cfg, "variant_overlay_seq2_color", "black"),
+                            variant_marker_size=getattr(cfg, "variant_overlay_marker_size", 4.0),
                             title=title,
                             save_name=out_png,
                         )
@@ -1000,7 +1049,9 @@ def chimeric_adata_core(
                         _vc = adata[sample_mask].layers[variant_call_layer_name]
                         _vc = _vc.toarray() if hasattr(_vc, "toarray") else np.asarray(_vc)
                         _cross_variant_call_df = pd.DataFrame(
-                            _vc, index=adata[sample_mask].obs_names.astype(str), columns=adata.var_names,
+                            _vc,
+                            index=adata[sample_mask].obs_names.astype(str),
+                            columns=adata.var_names,
                         )
 
                     # --- Plots ---
@@ -1073,6 +1124,15 @@ def chimeric_adata_core(
                                 variant_call_data=_cross_variant_call_df,
                                 seq1_label=variant_seq1_label,
                                 seq2_label=variant_seq2_label,
+                                ref1_marker_color=getattr(
+                                    cfg, "variant_overlay_seq1_color", "white"
+                                ),
+                                ref2_marker_color=getattr(
+                                    cfg, "variant_overlay_seq2_color", "black"
+                                ),
+                                variant_marker_size=getattr(
+                                    cfg, "variant_overlay_marker_size", 4.0
+                                ),
                                 title=title,
                                 save_name=out_png,
                             )
@@ -1146,12 +1206,31 @@ def chimeric_adata_core(
                 )
                 delta_layer = np.clip(within_layer - cross_layer, 0, None)
                 adata.layers[DELTA_ZERO_HAMMING_DISTANCE_SPANS] = delta_layer
+                threshold = getattr(cfg, "delta_hamming_chimeric_span_threshold", 200)
+                try:
+                    threshold = int(threshold)
+                except (TypeError, ValueError):
+                    logger.warning(
+                        "Invalid delta_hamming_chimeric_span_threshold=%s; using default 200.",
+                        threshold,
+                    )
+                    threshold = 200
+                if threshold < 0:
+                    logger.warning(
+                        "delta_hamming_chimeric_span_threshold=%s is negative; clamping to 0.",
+                        threshold,
+                    )
+                    threshold = 0
+                adata.obs["chimeric_by_mod_hamming_distance"] = (
+                    _compute_chimeric_by_mod_hamming_distance(delta_layer, threshold)
+                )
             else:
                 logger.warning(
                     "Cannot compute delta: missing %s or %s layer.",
                     ZERO_HAMMING_DISTANCE_SPANS,
                     CROSS_SAMPLE_ZERO_HAMMING_DISTANCE_SPANS,
                 )
+                adata.obs["chimeric_by_mod_hamming_distance"] = False
 
             if DELTA_ZERO_HAMMING_DISTANCE_SPANS in adata.layers:
                 for reference in references:
@@ -1300,7 +1379,9 @@ def chimeric_adata_core(
                     .astype("category")
                     .cat.categories.tolist()
                 )
-                references = adata.obs[cfg.reference_column].astype("category").cat.categories.tolist()
+                references = (
+                    adata.obs[cfg.reference_column].astype("category").cat.categories.tolist()
+                )
 
                 for reference in references:
                     ref_mask = adata.obs[cfg.reference_column] == reference
@@ -1332,9 +1413,7 @@ def chimeric_adata_core(
                         safe_sample = str(sample).replace(os.sep, "_")
                         safe_ref = str(reference).replace(os.sep, "_")
                         n_reads = int(sample_mask.sum())
-                        trio_title = (
-                            f"{sample} {reference} (n={n_reads})"
-                        )
+                        trio_title = f"{sample} {reference} (n={n_reads})"
                         out_png = span_trio_dir / f"{safe_sample}__{safe_ref}.png"
                         try:
                             plot_hamming_span_trio(
@@ -1345,6 +1424,15 @@ def chimeric_adata_core(
                                 variant_call_data=_variant_call_df,
                                 seq1_label=variant_seq1_label,
                                 seq2_label=variant_seq2_label,
+                                ref1_marker_color=getattr(
+                                    cfg, "variant_overlay_seq1_color", "white"
+                                ),
+                                ref2_marker_color=getattr(
+                                    cfg, "variant_overlay_seq2_color", "black"
+                                ),
+                                variant_marker_size=getattr(
+                                    cfg, "variant_overlay_marker_size", 4.0
+                                ),
                                 title=trio_title,
                                 save_name=out_png,
                             )

src/smftools/cli/hmm_adata.py

@@ -1089,6 +1089,10 @@ def hmm_adata_core(
                 deaminase=deaminase,
                 min_signal=0,
                 index_col_suffix=cfg.reindexed_var_suffix,
+                overlay_variant_calls=getattr(cfg, "overlay_variant_calls", False),
+                variant_overlay_seq1_color=getattr(cfg, "variant_overlay_seq1_color", "white"),
+                variant_overlay_seq2_color=getattr(cfg, "variant_overlay_seq2_color", "black"),
+                variant_overlay_marker_size=getattr(cfg, "variant_overlay_marker_size", 4.0),
             )
 
     hmm_length_dir = hmm_directory / "12b_hmm_length_clustermaps"
@@ -1145,6 +1149,10 @@ def hmm_adata_core(
                 min_signal=0,
                 index_col_suffix=cfg.reindexed_var_suffix,
                 length_feature_ranges=length_feature_ranges,
+                overlay_variant_calls=getattr(cfg, "overlay_variant_calls", False),
+                variant_overlay_seq1_color=getattr(cfg, "variant_overlay_seq1_color", "white"),
+                variant_overlay_seq2_color=getattr(cfg, "variant_overlay_seq2_color", "black"),
+                variant_overlay_marker_size=getattr(cfg, "variant_overlay_marker_size", 4.0),
             )
 
     hmm_dir = hmm_directory / "13_hmm_bulk_traces"

src/smftools/cli/variant_adata.py

@@ -221,13 +221,18 @@ def _find_converted_column(unconverted_col: str, var_columns) -> str | None:
                     prefix = f"{chrom}_"
                     end = f"_{strand}_{strand}{suffix}"
                     candidates = [
-                        c for c in var_columns
+                        c
+                        for c in var_columns
                         if c.startswith(prefix) and c.endswith(end) and c != unconverted_col
                     ]
                     if len(candidates) == 1:
                         return candidates[0]
                     if len(candidates) > 1:
-                        logger.info("Multiple converted column candidates for '%s': %s", unconverted_col, candidates)
+                        logger.info(
+                            "Multiple converted column candidates for '%s': %s",
+                            unconverted_col,
+                            candidates,
+                        )
                         return candidates[0]
                     break
             return None
@@ -252,6 +257,7 @@ def _find_converted_column(unconverted_col: str, var_columns) -> str | None:
             seq1_column=seq1_col,
             seq2_column=seq2_col,
             read_span_layer=cfg.mismatch_frequency_read_span_layer,
+            reference_col=cfg.reference_column,
         )
 
     ############################################### Plot mismatch base frequencies ###############################################
@@ -365,7 +371,38 @@ def _find_converted_column(unconverted_col: str, var_columns) -> str | None:
                     variant_segment_layer=segment_layer_name,
                     read_span_layer=cfg.mismatch_frequency_read_span_layer,
                     save_path=segment_dir,
+                    ref1_marker_color=getattr(cfg, "variant_overlay_seq1_color", "white"),
+                    ref2_marker_color=getattr(cfg, "variant_overlay_seq2_color", "black"),
+                    marker_size=getattr(cfg, "variant_overlay_marker_size", 4.0),
+                    show_position_axis=True,
+                )
+
+            segment_type_dir = (
+                variant_directory
+                / "deduplicated"
+                / "05_variant_segment_clustermaps_with_mismatch_type"
+            )
+            if segment_type_dir.exists():
+                logger.info(
+                    "Variant segment mismatch-type clustermaps already exist at %s; skipping.",
+                    segment_type_dir,
+                )
+            else:
+                make_dirs([segment_type_dir])
+                plot_variant_segment_clustermaps(
+                    adata,
+                    seq1_column=seq1_col,
+                    seq2_column=seq2_col,
+                    sample_col=cfg.sample_name_col_for_plotting,
+                    reference_col=cfg.reference_column,
+                    variant_segment_layer=segment_layer_name,
+                    read_span_layer=cfg.mismatch_frequency_read_span_layer,
+                    save_path=segment_type_dir,
+                    ref1_marker_color=getattr(cfg, "variant_overlay_seq1_color", "white"),
+                    ref2_marker_color=getattr(cfg, "variant_overlay_seq2_color", "black"),
+                    marker_size=getattr(cfg, "variant_overlay_marker_size", 4.0),
                     show_position_axis=True,
+                    mismatch_type_obs_col="chimeric_variant_sites_type",
                 )
 
     # ============================================================

src/smftools/config/default.yaml

@@ -188,6 +188,12 @@ clustermap_cmap_cpg: "coolwarm"
 clustermap_cmap_a: "coolwarm"
 spatial_clustermap_sortby: "gpc"
 
+# Clustermap variant params
+overlay_variant_calls: false
+variant_overlay_seq1_color: "black"
+variant_overlay_seq2_color: "white"
+variant_overlay_marker_size: 4.0
+
 # Spatial analysis - Rolling NN Hamming
 rolling_nn_layer: "nan0_0minus1"
 rolling_nn_plot_layer: "nan0_0minus1"
@@ -225,6 +231,7 @@ rolling_nn_zero_pairs_segment_histogram_bins: 30
 cross_sample_analysis: true
 cross_sample_grouping_col: null
 cross_sample_random_seed: 42
+delta_hamming_chimeric_span_threshold: 200
 
 # Latent Analysis - UMAP/Leiden params
 layer_for_umap_plotting: 'nan_half'
@@ -316,21 +323,13 @@ hmm_merge_layer_features:
   - ["all_accessible_features", 60]
 clustermap_cmap_hmm: "coolwarm"
 hmm_clustermap_feature_layers:
-  - all_accessible_features
   - all_accessible_features_merged
-  - small_accessible_patch
-  - mid_accessible_patch
-  - large_accessible_patch
-  - large_accessible_patch_merged
-  - nucleosome_depleted_region
   - nucleosome_depleted_region_merged
   - small_bound_stretch
   - medium_bound_stretch
   - putative_nucleosome
-  - large_bound_stretch
   - all_footprint_features
 hmm_clustermap_length_layers:
-  - all_accessible_features
   - all_accessible_features_merged
   - all_footprint_features
 hmm_clustermap_sortby: "hmm"