- Applied black formatting.

Author: BenjaminIsaac0111

scripts/batch_qc_stats.py

@@ -4,46 +4,75 @@
 import matplotlib.pyplot as plt
 from scipy.sparse import csr_matrix
 
+
 def calculate_qc_stats(h5ad_path, min_umis=500, min_genes=200, max_mt=0.15):
     with h5py.File(h5ad_path, "r") as f:
-        barcodes = [b.decode("utf-8") if isinstance(b, bytes) else b for b in f["obs"]["_index"][:]]
-        gene_names = [g.decode("utf-8") if isinstance(g, bytes) else g for g in f["var"]["_index"][:]]
+        barcodes = [
+            b.decode("utf-8") if isinstance(b, bytes) else b
+            for b in f["obs"]["_index"][:]
+        ]
+        gene_names = [
+            g.decode("utf-8") if isinstance(g, bytes) else g
+            for g in f["var"]["_index"][:]
+        ]
         X_group = f["X"]
         if isinstance(X_group, h5py.Group):
-            X = csr_matrix((X_group["data"][:], X_group["indices"][:], X_group["indptr"][:]), shape=(len(barcodes), len(gene_names)))
+            X = csr_matrix(
+                (X_group["data"][:], X_group["indices"][:], X_group["indptr"][:]),
+                shape=(len(barcodes), len(gene_names)),
+            )
         else:
             X = f["X"][:]
 
     n_counts = np.array(X.sum(axis=1)).flatten()
     n_genes = np.array((X > 0).sum(axis=1)).flatten()
-    mt_genes = [i for i, name in enumerate(gene_names) if "mt-" in name.lower() or "mt:" in name.lower()]
+    mt_genes = [
+        i
+        for i, name in enumerate(gene_names)
+        if "mt-" in name.lower() or "mt:" in name.lower()
+    ]
     if mt_genes:
         mt_counts = np.array(X[:, mt_genes].sum(axis=1)).flatten()
         pct_counts_mt = mt_counts / (n_counts + 1e-9)
     else:
         pct_counts_mt = np.zeros_like(n_counts)
 
-    keep_mask = (n_counts >= min_umis) & (n_genes >= min_genes) & (pct_counts_mt <= max_mt)
-    return len(barcodes), np.sum(keep_mask), np.sum(n_counts < min_umis), np.sum(n_genes < min_genes), np.sum(pct_counts_mt > max_mt)
+    keep_mask = (
+        (n_counts >= min_umis) & (n_genes >= min_genes) & (pct_counts_mt <= max_mt)
+    )
+    return (
+        len(barcodes),
+        np.sum(keep_mask),
+        np.sum(n_counts < min_umis),
+        np.sum(n_genes < min_genes),
+        np.sum(pct_counts_mt > max_mt),
+    )
+
 
 if __name__ == "__main__":
     st_dir = r"A:\hest_data\st"
-    
+
     # Get all h5ad files
-    samples = [f.replace(".h5ad", "") for f in os.listdir(st_dir) if f.endswith(".h5ad")]
+    samples = [
+        f.replace(".h5ad", "") for f in os.listdir(st_dir) if f.endswith(".h5ad")
+    ]
     samples.sort()
-    
+
     print(f"Analyzing {len(samples)} samples...")
-    print(f"{'Sample':<15} | {'Total':<6} | {'Kept':<6} | {'% Kept':<8} | {'Low UMI':<8} | {'Low Gene':<8} | {'High MT':<8}")
+    print(
+        f"{'Sample':<15} | {'Total':<6} | {'Kept':<6} | {'% Kept':<8} | {'Low UMI':<8} | {'Low Gene':<8} | {'High MT':<8}"
+    )
     print("-" * 80)
-    
+
     all_results = []
     for s in samples:
         path = os.path.join(st_dir, f"{s}.h5ad")
         try:
             total, kept, l_umi, l_gene, h_mt = calculate_qc_stats(path)
             pct = kept / total if total > 0 else 0
-            print(f"{s:<15} | {total:<6} | {kept:<6} | {pct:7.1%} | {l_umi:<8} | {l_gene:<8} | {h_mt:<8}")
+            print(
+                f"{s:<15} | {total:<6} | {kept:<6} | {pct:7.1%} | {l_umi:<8} | {l_gene:<8} | {h_mt:<8}"
+            )
             all_results.append((total, kept, pct))
         except Exception as e:
             print(f"Error processing {s}: {e}")

scripts/diagnose_qc.py

@@ -4,24 +4,33 @@
 import matplotlib.pyplot as plt
 from scipy.sparse import csr_matrix
 
+
 def diagnose_qc(h5ad_path, output_path, min_umis=500, min_genes=200, max_mt=0.15):
     print(f"Loading {h5ad_path}...")
-    
+
     with h5py.File(h5ad_path, "r") as f:
         # Load barcodes and gene names
-        barcodes = [b.decode("utf-8") if isinstance(b, bytes) else b for b in f["obs"]["_index"][:]]
-        gene_names = [g.decode("utf-8") if isinstance(g, bytes) else g for g in f["var"]["_index"][:]]
-        
+        barcodes = [
+            b.decode("utf-8") if isinstance(b, bytes) else b
+            for b in f["obs"]["_index"][:]
+        ]
+        gene_names = [
+            g.decode("utf-8") if isinstance(g, bytes) else g
+            for g in f["var"]["_index"][:]
+        ]
+
         # Load expression matrix X
         X_group = f["X"]
         if isinstance(X_group, h5py.Group):
             data = X_group["data"][:]
             indices = X_group["indices"][:]
             indptr = X_group["indptr"][:]
-            X = csr_matrix((data, indices, indptr), shape=(len(barcodes), len(gene_names)))
+            X = csr_matrix(
+                (data, indices, indptr), shape=(len(barcodes), len(gene_names))
+            )
         else:
             X = f["X"][:]
-            
+
         # Spatial coordinates (usually in obsm/spatial)
         if "obsm" in f and "spatial" in f["obsm"]:
             coords = f["obsm"]["spatial"][:]
@@ -33,10 +42,14 @@ def diagnose_qc(h5ad_path, output_path, min_umis=500, min_genes=200, max_mt=0.15
     # Calculate QC metrics
     n_counts = np.array(X.sum(axis=1)).flatten()
     n_genes = np.array((X > 0).sum(axis=1)).flatten()
-    
+
     # MT fraction
     # Robust MT detection: check for mt- or MT- anywhere, but prioritize common patterns
-    mt_genes = [i for i, name in enumerate(gene_names) if "mt-" in name.lower() or "mt:" in name.lower()]
+    mt_genes = [
+        i
+        for i, name in enumerate(gene_names)
+        if "mt-" in name.lower() or "mt:" in name.lower()
+    ]
     if mt_genes:
         print(f"Found {len(mt_genes)} mitochondrial genes.")
         mt_counts = np.array(X[:, mt_genes].sum(axis=1)).flatten()
@@ -49,12 +62,12 @@ def diagnose_qc(h5ad_path, output_path, min_umis=500, min_genes=200, max_mt=0.15
     pass_umi = n_counts >= min_umis
     pass_gene = n_genes >= min_genes
     pass_mt = pct_counts_mt <= max_mt
-    
+
     keep_mask = pass_umi & pass_gene & pass_mt
-    
+
     total_spots = len(barcodes)
     kept_spots = np.sum(keep_mask)
-    
+
     print(f"Total spots: {total_spots}")
     print(f"Kept spots: {kept_spots} ({kept_spots/total_spots:.1%})")
     print(f"Filtered: {total_spots - kept_spots}")
@@ -64,31 +77,41 @@ def diagnose_qc(h5ad_path, output_path, min_umis=500, min_genes=200, max_mt=0.15
 
     # Plotting
     fig, axes = plt.subplots(2, 2, figsize=(15, 12))
-    
+
     # 1. UMI vs Genes
     axes[0, 0].scatter(n_counts, n_genes, c=keep_mask, cmap="RdYlGn", alpha=0.5, s=10)
-    axes[0, 0].axvline(min_umis, color="red", linestyle="--", label=f"Min UMI={min_umis}")
-    axes[0, 0].axhline(min_genes, color="blue", linestyle="--", label=f"Min Genes={min_genes}")
+    axes[0, 0].axvline(
+        min_umis, color="red", linestyle="--", label=f"Min UMI={min_umis}"
+    )
+    axes[0, 0].axhline(
+        min_genes, color="blue", linestyle="--", label=f"Min Genes={min_genes}"
+    )
     axes[0, 0].set_xlabel("Total UMI counts")
     axes[0, 0].set_ylabel("Number of detected genes")
     axes[0, 0].set_title("QC: UMI vs Genes")
     axes[0, 0].legend()
-    
+
     # 2. MT Fraction distribution
     axes[0, 1].hist(pct_counts_mt, bins=50, color="gray", alpha=0.7)
-    axes[0, 1].axvline(max_mt, color="red", linestyle="--", label=f"Max MT={max_mt:.0%}")
+    axes[0, 1].axvline(
+        max_mt, color="red", linestyle="--", label=f"Max MT={max_mt:.0%}"
+    )
     axes[0, 1].set_xlabel("Mitochondrial Fraction")
     axes[0, 1].set_ylabel("Count")
     axes[0, 1].set_title("QC: MT Fraction Distribution")
     axes[0, 1].legend()
-    
+
     # 3. Spatial: Before (Total)
-    axes[1, 0].scatter(coords[:, 0], coords[:, 1], c="lightgray", s=15, label="All Spots")
-    axes[1, 0].scatter(coords[keep_mask, 0], coords[keep_mask, 1], c="green", s=15, label="Pass QC")
+    axes[1, 0].scatter(
+        coords[:, 0], coords[:, 1], c="lightgray", s=15, label="All Spots"
+    )
+    axes[1, 0].scatter(
+        coords[keep_mask, 0], coords[keep_mask, 1], c="green", s=15, label="Pass QC"
+    )
     axes[1, 0].set_title("Spatial Distribution: Kept vs Filtered")
     axes[1, 0].set_aspect("equal")
     axes[1, 0].legend()
-    
+
     # 4. Summary Table (as text)
     stats_text = (
         f"Sample: {os.path.basename(h5ad_path)}\n\n"
@@ -102,21 +125,23 @@ def diagnose_qc(h5ad_path, output_path, min_umis=500, min_genes=200, max_mt=0.15
     )
     axes[1, 1].text(0.1, 0.5, stats_text, fontsize=14, family="monospace")
     axes[1, 1].axis("off")
-    
+
     plt.tight_layout()
     plt.savefig(output_path)
     print(f"Plot saved to {output_path}")
 
+
 if __name__ == "__main__":
     import argparse
+
     parser = argparse.ArgumentParser()
     parser.add_argument("--sample", type=str, default="MEND29")
     args = parser.parse_args()
-    
+
     sample_id = args.sample
     h5ad_file = f"A:\\hest_data\\st\\{sample_id}.h5ad"
     output_file = f"qc_diagnosis_{sample_id}.png"
-    
+
     if not os.path.exists(h5ad_file):
         print(f"Error: {h5ad_file} not found.")
     else:

scripts/find_outliers.py

@@ -1,5 +1,6 @@
 import os
 import sys
+
 sys.path.append(r"z:\Projects\SpatialTranscriptFormer\scripts")
 from batch_qc_stats import calculate_qc_stats
 import numpy as np
@@ -11,22 +12,30 @@
 results = []
 for s in samples:
     try:
-        total, kept, l_umi, l_gene, h_mt = calculate_qc_stats(os.path.join(st_dir, f"{s}.h5ad"))
-        results.append({
-            "sample": s,
-            "total": total,
-            "kept": kept,
-            "pct": kept / total if total > 0 else 0,
-            "low_umi": l_umi,
-            "low_gene": l_gene,
-            "high_mt": h_mt
-        })
+        total, kept, l_umi, l_gene, h_mt = calculate_qc_stats(
+            os.path.join(st_dir, f"{s}.h5ad")
+        )
+        results.append(
+            {
+                "sample": s,
+                "total": total,
+                "kept": kept,
+                "pct": kept / total if total > 0 else 0,
+                "low_umi": l_umi,
+                "low_gene": l_gene,
+                "high_mt": h_mt,
+            }
+        )
     except Exception as e:
         print(f"Error {s}: {e}")
 
 results.sort(key=lambda x: x["pct"])
 
-print(f"{'Sample':<15} | {'Kept %':<8} | {'Filtered':<10} | {'Low UMI':<8} | {'Low Gene':<8} | {'High MT':<8}")
+print(
+    f"{'Sample':<15} | {'Kept %':<8} | {'Filtered':<10} | {'Low UMI':<8} | {'Low Gene':<8} | {'High MT':<8}"
+)
 print("-" * 75)
 for r in results[:15]:
-    print(f"{r['sample']:<15} | {r['pct']:7.1%} | {r['total']-r['kept']:<10} | {r['low_umi']:<8} | {r['low_gene']:<8} | {r['high_mt']:<8}")
+    print(
+        f"{r['sample']:<15} | {r['pct']:7.1%} | {r['total']-r['kept']:<10} | {r['low_umi']:<8} | {r['low_gene']:<8} | {r['high_mt']:<8}"
+    )

scripts/monitor.py

@@ -2,6 +2,7 @@
 """
 Real-time Training Monitor Entrypoint for SpatialTranscriptFormer.
 """
+
 import argparse
 import logging
 from spatial_transcript_former.dashboard.app import init_app, app

src/spatial_transcript_former/checkpoint.py

@@ -13,7 +13,6 @@
 
 import torch
 
-
 # Keys serialized into config.json.  These correspond to
 # SpatialTranscriptFormer.__init__ arguments (minus runtime-only
 # arguments like ``pathway_init`` and ``pretrained``).

src/spatial_transcript_former/data/pathways.py

@@ -18,6 +18,7 @@
     "c2_cgp": "https://data.broadinstitute.org/gsea-msigdb/msigdb/release/2024.1.Hs/c2.cgp.v2024.1.Hs.symbols.gmt",
 }
 
+
 def download_msigdb_gmt(url: str, filename: str, cache_dir: str = ".cache") -> str:
     """
     Download an MSigDB GMT file if not already cached.

src/spatial_transcript_former/data/spatial_stats.py

@@ -71,7 +71,7 @@ def morans_i(x: np.ndarray, W: csr_matrix) -> float:
     x_mean = x.mean()
     z = x - x_mean
 
-    denominator = np.sum(z ** 2)
+    denominator = np.sum(z**2)
     if denominator < 1e-12:
         return 0.0  # Constant expression → no spatial pattern
 

src/spatial_transcript_former/models/interaction.py

@@ -348,9 +348,7 @@ def forward(
                     x_layer = x_layer + layer.dropout1(attn_output)
                     x_layer = x_layer + layer._ff_block(layer.norm2(x_layer))
                 else:
-                    x_layer = layer.norm1(
-                        x_layer + layer.dropout1(attn_output)
-                    )
+                    x_layer = layer.norm1(x_layer + layer.dropout1(attn_output))
                     x_layer = layer.norm2(x_layer + layer._ff_block(x_layer))
             out = x_layer
         else:

src/spatial_transcript_former/predict.py

@@ -22,7 +22,6 @@
 import torch.nn as nn
 from torchvision import transforms
 
-
 # ═══════════════════════════════════════════════════════════════════════
 #  FeatureExtractor
 # ═══════════════════════════════════════════════════════════════════════

src/spatial_transcript_former/recipes/hest/build_vocab.py

@@ -44,8 +44,12 @@ def scan_h5ad_files(data_dir):
 
 
 def calculate_global_genes(
-    data_dir, ids, num_genes=1000, target_pathways=None,
-    svg_weight=0.0, svg_k=6,
+    data_dir,
+    ids,
+    num_genes=1000,
+    target_pathways=None,
+    svg_weight=0.0,
+    svg_k=6,
 ):
     st_dir = os.path.join(data_dir, "st")
     if not ids:
@@ -158,15 +162,13 @@ def calculate_global_genes(
         # Hybrid score: weighted sum of ranks (lower = better)
         alpha = svg_weight
         hybrid_score = {
-            g: (1 - alpha) * expr_rank[g] + alpha * mi_rank[g]
-            for g in all_genes
+            g: (1 - alpha) * expr_rank[g] + alpha * mi_rank[g] for g in all_genes
         }
         sorted_all_genes = sorted(all_genes, key=lambda g: hybrid_score[g])
 
         # Build stats list with Moran's I column
         sorted_all = [
-            (g, gene_totals[g], gene_morans_avg.get(g, 0.0))
-            for g in sorted_all_genes
+            (g, gene_totals[g], gene_morans_avg.get(g, 0.0)) for g in sorted_all_genes
         ]
         print(
             f"Hybrid ranking: expression weight={(1 - alpha):.1f}, "
@@ -221,7 +223,7 @@ def main():
         type=float,
         default=0.0,
         help="Weight for spatial variability (Moran's I) in gene ranking. "
-             "0.0=expression-only (default), 1.0=SVG-only, 0.5=balanced.",
+        "0.0=expression-only (default), 1.0=SVG-only, 0.5=balanced.",
     )
     parser.add_argument(
         "--svg-k",
@@ -242,8 +244,12 @@ def main():
         sys.exit(1)
 
     top_genes, all_stats = calculate_global_genes(
-        args.data_dir, ids, args.num_genes, target_pathways=args.pathways,
-        svg_weight=args.svg_weight, svg_k=args.svg_k,
+        args.data_dir,
+        ids,
+        args.num_genes,
+        target_pathways=args.pathways,
+        svg_weight=args.svg_weight,
+        svg_k=args.svg_k,
     )
 
     print(f"Saving top {len(top_genes)} genes to {output_path}")