Merge pull request #182 from jkmckenna/0.2.5

0.2.5

Author: jkmckenna

README.md

@@ -12,22 +12,3 @@ The following CLI tools need to be installed and configured before using the inf
 1) [Dorado](https://github.com/nanoporetech/dorado) -> Basecalling, alignment, demultiplexing. Required for Nanopore SMF experiments, but not Illumina SMF experiments.
 2) [Minimap2](https://github.com/lh3/minimap2) -> Aligner if not using dorado. Support for other aligners could eventually be added if needed.
 3) [Modkit](https://github.com/nanoporetech/modkit) -> Extracting read level methylation metrics from the MM/ML tags in BAM files. Only required for direct modification detection SMF protocols.
-
-## Announcements
-
-### 12/02/25 - Version 0.2.3 is available through PyPI
-Version 0.2.3 provides the core smftools functionality through several command line commands (load, preprocess, spatial, hmm).
-
-### 11/05/25 - Version 0.2.1 is available through PyPI
-Version 0.2.1 makes the core workflow (smftools load) a command line tool that takes in an experiment_config.csv file for input/output and parameter management.
-
-### 05/29/25 - Version 0.1.6 is available through PyPI.
-Informatics, preprocessing, tools, plotting modules have core functionality that is approaching stability on MacOS(Intel/Silicon) and Linux(Ubuntu). I will work on improving documentation/tutorials shortly. The base PyTorch/Scikit-Learn ML-infrastructure is going through some organizational changes to work with PyTorch Lightning, Hydra, and WanDB to facilitate organizational scaling, multi-device usage, and logging.
-
-### 10/01/24 - More recent versions are being updated frequently. Installation from source over PyPI is recommended!
-
-### 09/09/24 - The version 0.1.1 package ([smftools-0.1.1](https://pypi.org/project/smftools/)) is installable through pypi!
-The informatics module has been bumped to alpha-phase status. This module can deal with POD5s and unaligned BAMS from nanopore conversion and direct SMF experiments, as well as FASTQs from Illumina conversion SMF experiments. Primary output from this module is an AnnData object containing all relevant SMF data, which is compatible with all downstream smftools modules. The other modules are still in pre-alpha phase. Preprocessing, Tools, and Plotting modules should be promoted to alpha-phase within the next month or so.
-
-### 08/30/24 - The version 0.1.0 package ([smftools-0.1.0](https://pypi.org/project/smftools/)) is installable through pypi!
-Currently, this package (smftools-0.1.0) is going through rapid improvement (dependency handling accross Linux and Mac OS, testing, documentation, debugging) and is still too early in development for widespread use. The underlying functionality was originally developed as a collection of scripts for single molecule footprinting (SMF) experiments in our lab, but is being packaged/developed to facilitate the expansion of SMF to any lab that is interested in performing these styles of experiments/analyses. The alpha-phase package is expected to be available within a couple months, so stay tuned!

docs/source/basic_usage.md

@@ -13,7 +13,7 @@ This command takes a user passed config file handling:
     - Experiment info (SMF modality, sequencer type, barcoding kit if nanopore, sample sheet with metadata mapping)
     - Options to override default workflow parameters from smftools/config. Params are handled from default.yaml -> modality_type.yaml -> user passed config.csv.
 
-![](docs/source/_static/smftools_informatics_diagram.png)
+![](_static/smftools_informatics_diagram.png)
 
 ## Preprocess Usage
 
@@ -23,7 +23,7 @@ This command performs preprocessing on the anndata object. It automatically runs
 smftools preprocess "/Path_to_experiment_config.csv"
 ```
 
-![](docs/source/_static/smftools_preprocessing_diagram.png)
+![](_static/smftools_preprocessing_diagram.png)
 
 ## Spatial Usage
 

docs/source/release-notes/0.1.0.md

@@ -1,4 +1,5 @@
 (v0.1.0)=
 
 0.1.0 2024-08-30
-smftools initial release. Pre-Alpha phase. J McKenna
+smftools initial release. Pre-Alpha phase. J McKenna
+Currently, this package (smftools-0.1.0) is going through rapid improvement (dependency handling accross Linux and Mac OS, testing, documentation, debugging) and is still too early in development for widespread use. The underlying functionality was originally developed as a collection of scripts for single molecule footprinting (SMF) experiments in our lab, but is being packaged/developed to facilitate the expansion of SMF to any lab that is interested in performing these styles of experiments/analyses. The alpha-phase package is expected to be available within a couple months, so stay tuned!

docs/source/release-notes/0.1.1.md

@@ -0,0 +1,4 @@
+(v0.1.1)=
+
+0.1.1 2024-09-09
+The informatics module has been bumped to alpha-phase status. This module can deal with POD5s and unaligned BAMS from nanopore conversion and direct SMF experiments, as well as FASTQs from Illumina conversion SMF experiments. Primary output from this module is an AnnData object containing all relevant SMF data, which is compatible with all downstream smftools modules. The other modules are still in pre-alpha phase. Preprocessing, Tools, and Plotting modules should be promoted to alpha-phase within the next month or so.

docs/source/release-notes/0.1.6.md

@@ -0,0 +1,4 @@
+(v0.1.6)=
+
+0.1.6 2025-05-29
+Informatics, preprocessing, tools, plotting modules have core functionality that is approaching stability on MacOS (Intel/Silicon) and Linux (Ubuntu). Documentation/tutorials are still being improved. The base PyTorch/Scikit-Learn ML-infrastructure is going through organizational changes to work with PyTorch Lightning, Hydra, and WanDB to facilitate scaling, multi-device usage, and logging.

docs/source/release-notes/0.2.1.md

@@ -0,0 +1,4 @@
+(v0.2.1)=
+
+0.2.1 2025-11-05
+Version 0.2.1 makes the core workflow (smftools load) a command line tool that takes in an experiment_config.csv file for input/output and parameter management.

docs/source/release-notes/0.2.3.md

@@ -0,0 +1,4 @@
+(v0.2.3)=
+
+0.2.3 2025-12-02
+Version 0.2.3 provides the core smftools functionality through several command line commands (load, preprocess, spatial, hmm).

src/smftools/hmm/HMM.py

@@ -10,6 +10,9 @@
 import torch.nn as nn
 from scipy.sparse import issparse
 
+from smftools.logging_utils import get_logger
+
+logger = get_logger(__name__)
 # =============================================================================
 # Registry / Factory
 # =============================================================================
@@ -1228,7 +1231,11 @@ def fit_em(
             self._normalize_emission()
 
             if verbose:
-                print(f"[SingleBernoulliHMM.fit] iter={it} ll_proxy={hist[-1]:.6f}")
+                logger.info(
+                    "[SingleBernoulliHMM.fit] iter=%s ll_proxy=%.6f",
+                    it,
+                    hist[-1],
+                )
 
             if len(hist) > 1 and abs(hist[-1] - hist[-2]) < float(tol):
                 break
@@ -1450,7 +1457,11 @@ def fit_em(
             self._normalize_emission()
 
             if verbose:
-                print(f"[MultiBernoulliHMM.fit] iter={it} ll_proxy={hist[-1]:.6f}")
+                logger.info(
+                    "[MultiBernoulliHMM.fit] iter=%s ll_proxy=%.6f",
+                    it,
+                    hist[-1],
+                )
 
             if len(hist) > 1 and abs(hist[-1] - hist[-2]) < float(tol):
                 break
@@ -1783,7 +1794,11 @@ def fit_em(
             self._normalize_trans_by_bin()
 
             if verbose:
-                print(f"[DistanceBinnedSingle.fit] iter={it} ll_proxy={hist[-1]:.6f}")
+                logger.info(
+                    "[DistanceBinnedSingle.fit] iter=%s ll_proxy=%.6f",
+                    it,
+                    hist[-1],
+                )
 
             if len(hist) > 1 and abs(hist[-1] - hist[-2]) < float(tol):
                 break

src/smftools/hmm/call_hmm_peaks.py

@@ -3,6 +3,10 @@
 from pathlib import Path
 from typing import Any, Dict, Optional, Sequence, Union
 
+from smftools.logging_utils import get_logger
+
+logger = get_logger(__name__)
+
 
 def call_hmm_peaks(
     adata,
@@ -107,8 +111,10 @@ def call_hmm_peaks(
                 candidates = [feature_key]
 
             if not candidates:
-                print(
-                    f"[call_hmm_peaks] WARNING: no layers found matching '{feature_key}' in ref '{ref}'. Skipping."
+                logger.warning(
+                    "[call_hmm_peaks] No layers found matching '%s' in ref '%s'. Skipping.",
+                    feature_key,
+                    ref,
                 )
                 continue
 
@@ -121,17 +127,22 @@ def call_hmm_peaks(
 
             for layer_name in candidates:
                 if layer_name not in adata.layers:
-                    print(
-                        f"[call_hmm_peaks] WARNING: layer '{layer_name}' not in adata.layers; skipping."
+                    logger.warning(
+                        "[call_hmm_peaks] Layer '%s' not in adata.layers; skipping.",
+                        layer_name,
                     )
                     continue
 
                 # Dense layer data
                 L = adata.layers[layer_name]
                 L = L.toarray() if issparse(L) else np.asarray(L)
                 if L.shape != (adata.n_obs, adata.n_vars):
-                    print(
-                        f"[call_hmm_peaks] WARNING: layer '{layer_name}' has shape {L.shape}, expected ({adata.n_obs}, {adata.n_vars}); skipping."
+                    logger.warning(
+                        "[call_hmm_peaks] Layer '%s' has shape %s, expected (%s, %s); skipping.",
+                        layer_name,
+                        L.shape,
+                        adata.n_obs,
+                        adata.n_vars,
                     )
                     continue
 
@@ -154,7 +165,11 @@ def call_hmm_peaks(
                     peak_metric, prominence=peak_prom, distance=min_distance
                 )
                 if peak_indices.size == 0:
-                    print(f"[call_hmm_peaks] No peaks for layer '{layer_name}' in ref '{ref}'.")
+                    logger.info(
+                        "[call_hmm_peaks] No peaks for layer '%s' in ref '%s'.",
+                        layer_name,
+                        ref,
+                    )
                     continue
 
                 peak_centers = coordinates[peak_indices]
@@ -185,7 +200,7 @@ def call_hmm_peaks(
                     safe_layer = str(layer_name).replace("/", "_")
                     fname = output_dir / f"{tag}_{safe_layer}_{safe_ref}_peaks.png"
                     fig.savefig(fname, bbox_inches="tight", dpi=200)
-                    print(f"[call_hmm_peaks] Saved plot to {fname}")
+                    logger.info("[call_hmm_peaks] Saved plot to %s", fname)
                     plt.close(fig)
                 else:
                     fig.tight_layout()
@@ -285,8 +300,11 @@ def call_hmm_peaks(
                 else:
                     adata.var[any_col] = False
 
-                print(
-                    f"[call_hmm_peaks] Annotated {len(peak_centers)} peaks for layer '{layer_name}' in ref '{ref}'."
+                logger.info(
+                    "[call_hmm_peaks] Annotated %s peaks for layer '%s' in ref '%s'.",
+                    len(peak_centers),
+                    layer_name,
+                    ref,
                 )
 
     # global any-peak across all layers/refs

src/smftools/hmm/display_hmm.py

@@ -1,19 +1,24 @@
+from smftools.logging_utils import get_logger
+
+logger = get_logger(__name__)
+
+
 def display_hmm(hmm, state_labels=["Non-Methylated", "Methylated"], obs_labels=["0", "1"]):
     import torch
 
-    print("\n**HMM Model Overview**")
-    print(hmm)
+    logger.info("**HMM Model Overview**")
+    logger.info("%s", hmm)
 
-    print("\n**Transition Matrix**")
+    logger.info("**Transition Matrix**")
     transition_matrix = torch.exp(hmm.edges).detach().cpu().numpy()
     for i, row in enumerate(transition_matrix):
         label = state_labels[i] if state_labels else f"State {i}"
         formatted_row = ", ".join(f"{p:.6f}" for p in row)
-        print(f"{label}: [{formatted_row}]")
+        logger.info("%s: [%s]", label, formatted_row)
 
-    print("\n**Emission Probabilities**")
+    logger.info("**Emission Probabilities**")
     for i, dist in enumerate(hmm.distributions):
         label = state_labels[i] if state_labels else f"State {i}"
         probs = dist.probs.detach().cpu().numpy()
         formatted_emissions = {obs_labels[j]: probs[j] for j in range(len(probs))}
-        print(f"{label}: {formatted_emissions}")
+        logger.info("%s: %s", label, formatted_emissions)