Feature/dge in browser

DGE_README.md

@@ -0,0 +1,93 @@
+# MDV Differential Gene Expression (DGE) Implementation
+
+This document describes the "Tier 1" browser-based Differential Gene Expression (DGE) implementation in MDV.
+
+## 1. Overview
+
+The MDV DGE engine (`mdv-dge`) provides high-performance, real-time statistical analysis of single-cell RNA-seq data directly in the browser. It is designed to work with hundreds of thousands of cells and tens of thousands of genes without requiring a dedicated high-memory backend for computation.
+
+## 2. Architecture
+
+The DGE feature is split into four layers:
+
+1.  **UI Layer (`DGEDialogReact.tsx`)**: A React/MUI dialog that allows users to select a grouping column (e.g., Leiden cluster, Disease state) and compare a target group against a reference group (or "rest").
+2.  **Integration Layer (`dgeIntegration.ts`)**: The "glue" that connects MDV's `DataStore` and `ChartManager` to the DGE engine. It handles data discovery via `rows_as_columns` links, coordinates batch loading of gene data, and writes results back as new columns.
+3.  **Calculation Engine (`DGEDimension.ts` & `dgeWorker.ts`)**: Orchestrates the DGE pipeline. It divides genes into batches and can optionally offload calculations to **Web Workers** to keep the UI responsive.
+4.  **Statistical Core (`dgeStats.ts`)**: Pure TypeScript implementation of the statistical algorithms (Welch's t-test, Welford's online variance, Benjamini-Hochberg correction).
+
+## 3. Statistical Methodology
+
+### Welch's t-test
+We use **Welch’s t-test** (unequal variances t-test) to compare gene expression between two groups. This is the same default method used by popular tools like **Scanpy**.
+
+*   **P-value Calculation**: Computed using the Student's t-distribution CDF (via the Regularized Incomplete Beta function).
+*   **Multiple Testing Correction**: P-values are adjusted using the **Benjamini-Hochberg (BH)** procedure to control the False Discovery Rate (FDR).
+
+### Online Variance (Welford's Algorithm)
+To handle large datasets without multiple passes through the data, we use **Welford’s online algorithm**. This allows us to calculate the mean and variance of gene expression in a single pass over the cell indices, minimizing memory pressure.
+
+### Adaptive Effect Size
+The engine automatically detects the normalization of the input data and selects the appropriate effect size metric:
+
+*   **Log1p-normalized data**: Uses **Log2 Fold Change (Log2FC)** with `expm1` back-transformation.
+*   **Linear / raw count data**: Uses **Log2 Fold Change (Log2FC)** directly on means (no `expm1`).
+*   **Z-scored data**: Uses **Mean Difference** (`meanTarget - meanReference`).
+
+Detection is automatic via a multi-gene probe (see gotcha below).
+
+### Gotcha: Data Type Detection for Sparse Raw Counts
+
+> **If the volcano plot shows many genes clamped at exactly ±10 effect size, or highly-expressed housekeeping genes have implausible log2FC values, check that data type detection is working correctly.**
+
+The engine probes the first 50 genes to decide if expression values are log1p-normalized (small non-negative, max ≤ 20), linear/raw counts (non-negative, max > 20), or z-scored (has negatives). For **sparse raw count data**, most genes have very low counts (1–3 per cell) which look identical to log1p-normalized values. If detection breaks at the first gene it finds, it may misclassify raw counts as log1p — causing `expm1()` to be applied to already-linear values, which exponentially inflates fold changes for highly-expressed genes.
+
+**The fix** (implemented in `dgeIntegration.ts`): the probe loop scans across **all** 50 sampled genes, tracking the global maximum value, and only decides after examining all of them. It exits early only when it finds a definitive signal (a negative value → z-scored, or a value > 20 → linear).
+
+Even with correct detection, genes expressed exclusively in one group (e.g., 72 UC cells, 0 Healthy cells) will produce mathematically extreme fold changes that are clamped to ±10. This is expected for very sparse genes and can be mitigated with a minimum cell count filter if needed.
+
+## 4. Performance & Scalability
+
+### Batching
+Data is requested from the backend in batches (default: **2000 genes per batch**). This strikes a balance between network overhead and browser memory usage.
+
+### Memory Optimization
+*   **SharedArrayBuffer**: Expression data is loaded into `SharedArrayBuffer` objects, allowing efficient sharing between the main thread and Web Workers.
+*   **Sparse Data Handling**: In sparse datasets, `NaN` values in MDV's dense buffers are treated as `0` during statistical accumulation, ensuring correct cell counts for variance calculations.
+
+### Browser Limits
+*   **Memory**: Each gene for 60k cells takes ~240KB. A batch of 2000 genes takes ~480MB.
+*   **Dataset Size**: While Tier 1 works well for up to 500k cells, datasets of **10M+ cells** require either sub-sampling or server-side execution (Tier 2) due to browser memory and transfer limits.
+
+## 5. Backend Integration: HTTP Streaming
+
+To prevent the backend from crashing when requesting thousands of columns at once, we implemented **HTTP Streaming** (Chunked Transfer Encoding) in the Flask server:
+
+*   **`yield_byte_data`**: The Python backend streams individual gene columns from HDF5 as they are read, instead of buffering the entire multi-gigabyte response in RAM.
+*   **No Content-Length**: This removes the risk of `ERR_CONTENT_LENGTH_MISMATCH` errors caused by large binary payloads.
+
+## 6. Configuration & Troubleshooting
+
+### Adjusting Batch Size
+If the UI feels sluggish or you encounter memory errors, you can adjust the default DGE batch size in the following centralized location:
+
+*   **File:** `src/lib/constants.ts`
+*   **Constant:** `DEFAULT_DGE_BATCH_SIZE`
+
+This single source of truth is used as the default fallback in both the integration layer (`dgeIntegration.ts`) and the core DGE engine (`DGEDimension.ts`).
+
+**After making changes** you must perform a hard browser refresh (Cmd+Shift+R or Ctrl+F5) to clear the cached JavaScript.
+
+#### Guidance on choosing a value
+
+| Batch Size | Effect |
+| :--- | :--- |
+| **500** | Safer for large datasets (>100k cells). Lower peak memory per batch (~120MB). More HTTP round-trips. |
+| **1000** | Good balance for datasets up to ~200k cells. |
+| **2000** *(default)* | Optimal for datasets up to ~60k cells. Reduces round-trips significantly. |
+| **5000+** | Only suitable for small datasets (<20k cells). Risk of browser memory pressure / GC pauses. |
+
+### Visualizing Results
+After a DGE run completes, a **Volcano Plot** is automatically generated. You can also manually create scatter plots or color other charts using the newly generated `dge_effect_size` and `dge_neg_log10_pval_adj` columns.
+
+### Caching
+DGE results are cached in memory. If you run DGE a second time on the same dataset, the browser will reuse the cached gene expression data, making the second run nearly instantaneous.

python/mdvtools/mdvproject.py

@@ -1755,6 +1755,7 @@ def get_byte_data(self, columns, group):
         byte_list = []
         gr = h5[group]
         if not isinstance(gr, h5py.Group):
+            h5.close()
             raise TypeError("Expected 'gr' to be of type h5py.Group.")
         for column in columns:
             sg = column.get("subgroup")
@@ -1771,6 +1772,27 @@ def get_byte_data(self, columns, group):
         h5.close()
         return b"".join(byte_list)
 
+    def yield_byte_data(self, columns, group):
+        """
+        Generator version of get_byte_data for HTTP streaming.
+        Yields bytes for each column one by one to avoid large memory allocations
+        and Content-Length mismatches.
+        """
+        with h5py.File(self.h5file, "r") as h5:
+            gr = h5[group]
+            if not isinstance(gr, h5py.Group):
+                raise TypeError("Expected 'gr' to be of type h5py.Group.")
+            for column in columns:
+                sg = column.get("subgroup")
+                if sg:
+                    sgindex = int(column["sgindex"])
+                    yield get_subgroup_bytes(
+                        gr[sg], sgindex, column.get("sgtype") == "sparse"
+                    )
+                else:
+                    data = gr[column["field"]]
+                    yield numpy.array(data).tobytes()
+
     def set_region_data(
         self,
         datasource,

python/mdvtools/server.py

@@ -5,7 +5,9 @@
     request,
     make_response,
     jsonify,
-    current_app
+    current_app,
+    Response,
+    stream_with_context
 )
 from mdvtools.server_utils import (
     send_file,
@@ -196,10 +198,12 @@ def get_data():
                 raise Exception(
                     "Request must contain JSON with 'columns' and 'data_source'"
                 )
-            bytes_ = project.get_byte_data(data["columns"], data["data_source"])
-            response = make_response(bytes_)
-            response.headers.set("Content-Type", "application/octet-stream")
-            return response
+
+            # Use streaming to avoid large memory allocations and Content-Length mismatches
+            return Response(
+                stream_with_context(project.yield_byte_data(data["columns"], data["data_source"])),
+                mimetype="application/octet-stream"
+            )
         except Exception as e:
             log(e)
             return "Problem handling request", 400

python/mdvtools/tests/generate_dge_reference.py

@@ -0,0 +1,142 @@
+"""Generate Scanpy DGE reference data for mdv-dge validation.
+
+Runs rank_genes_groups on the PBMC3k dataset for all leiden clusters vs rest,
+exports results as JSON fixtures for use in vitest browser-side DGE tests.
+
+Usage:
+    cd /app/python && poetry run python -m mdvtools.tests.generate_dge_reference
+"""
+
+import json
+import sys
+import numpy as np
+import scanpy as sc
+
+
+H5AD_PATH = "/app/mdv/16/scanpy-pbmc3k.h5ad"
+FIXTURE_DIR = "/app/src/tests/fixtures"
+TOP_N_SMALL = 50
+VALIDATION_GENES = 100
+
+
+def extract_results(adata, method: str) -> dict:
+    """Extract rank_genes_groups results into a serializable dict."""
+    result = sc.get.rank_genes_groups_df(adata, group=None)
+
+    groups = {}
+    for group_name in result["group"].unique():
+        gdf = result[result["group"] == group_name].copy()
+        gdf = gdf.sort_values("pvals")
+        groups[str(group_name)] = {
+            "genes": gdf["names"].tolist(),
+            "scores": [float(x) for x in gdf["scores"]],
+            "log2fc": [float(x) for x in gdf["logfoldchanges"]],
+            "pvals": [float(x) for x in gdf["pvals"]],
+            "pvals_adj": [float(x) for x in gdf["pvals_adj"]],
+        }
+
+    group_sizes = adata.obs["leiden"].value_counts().to_dict()
+
+    return {
+        "method": method,
+        "groupby": "leiden",
+        "reference": "rest",
+        "n_cells": int(adata.n_obs),
+        "n_genes": int(adata.n_vars),
+        "group_sizes": {str(k): int(v) for k, v in group_sizes.items()},
+        "groups": groups,
+    }
+
+
+def extract_small(full_results: dict) -> dict:
+    """Extract top N genes per group for fast unit tests."""
+    small = {k: v for k, v in full_results.items() if k != "groups"}
+    small["top_n"] = TOP_N_SMALL
+    small["groups"] = {}
+    for gname, gdata in full_results["groups"].items():
+        small["groups"][gname] = {
+            k: v[:TOP_N_SMALL] for k, v in gdata.items()
+        }
+    return small
+
+
+def extract_validation_data(adata) -> dict:
+    """Extract log1p-normalized expression data for a subset of genes + leiden assignments.
+
+    Uses adata.raw.X (log1p-normalized) since that's what rank_genes_groups uses
+    by default. This allows the JS test to recompute DGE from scratch and compare.
+    """
+    # Use raw data (log1p-normalized) to match rank_genes_groups default behavior
+    raw = adata.raw
+    raw_gene_names = raw.var_names.tolist()
+
+    # Select genes evenly spaced from the raw gene set, but only those that
+    # appear in the rank_genes_groups results (i.e. all raw genes)
+    selected_indices = np.linspace(0, len(raw_gene_names) - 1, VALIDATION_GENES, dtype=int)
+    selected_genes = [raw_gene_names[i] for i in selected_indices]
+
+    X = raw.X
+    if hasattr(X, "toarray"):
+        X = X.toarray()
+    expression = {}
+    for idx, gene in zip(selected_indices, selected_genes):
+        expression[gene] = [float(x) for x in X[:, idx]]
+
+    leiden = adata.obs["leiden"].tolist()
+
+    return {
+        "n_cells": int(adata.n_obs),
+        "n_genes_total": int(raw.n_vars),
+        "n_genes_selected": len(selected_genes),
+        "genes": selected_genes,
+        "leiden": leiden,
+        "expression": expression,
+    }
+
+
+def write_json(data: dict, path: str):
+    with open(path, "w") as f:
+        json.dump(data, f, allow_nan=False, separators=(",", ":"))
+    size_mb = len(json.dumps(data, separators=(",", ":"))) / 1024 / 1024
+    print(f"  Written {path} ({size_mb:.1f} MB)")
+
+
+def main():
+    print(f"Loading {H5AD_PATH}...")
+    adata = sc.read_h5ad(H5AD_PATH)
+    print(f"  {adata.n_obs} cells, {adata.n_vars} genes")
+    print(f"  Leiden clusters: {sorted(adata.obs['leiden'].unique())}")
+
+    # --- T-test ---
+    print("\nRunning rank_genes_groups (t-test)...")
+    sc.tl.rank_genes_groups(adata, "leiden", method="t-test", reference="rest")
+    ttest_full = extract_results(adata, "t-test")
+    ttest_small = extract_small(ttest_full)
+
+    write_json(ttest_full, f"{FIXTURE_DIR}/dge_reference_ttest.json")
+    write_json(ttest_small, f"{FIXTURE_DIR}/dge_reference_ttest_small.json")
+
+    # --- Wilcoxon ---
+    print("\nRunning rank_genes_groups (wilcoxon)...")
+    sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon", reference="rest")
+    wilcoxon_full = extract_results(adata, "wilcoxon")
+    wilcoxon_small = extract_small(wilcoxon_full)
+
+    write_json(wilcoxon_full, f"{FIXTURE_DIR}/dge_reference_wilcoxon.json")
+    write_json(wilcoxon_small, f"{FIXTURE_DIR}/dge_reference_wilcoxon_small.json")
+
+    # --- Validation expression data ---
+    print("\nExtracting validation expression data...")
+    validation = extract_validation_data(adata)
+    write_json(validation, f"{FIXTURE_DIR}/dge_validation_data.json")
+
+    print("\nDone. Generated fixtures:")
+    print(f"  {FIXTURE_DIR}/dge_reference_ttest.json")
+    print(f"  {FIXTURE_DIR}/dge_reference_ttest_small.json")
+    print(f"  {FIXTURE_DIR}/dge_reference_wilcoxon.json")
+    print(f"  {FIXTURE_DIR}/dge_reference_wilcoxon_small.json")
+    print(f"  {FIXTURE_DIR}/dge_validation_data.json")
+
+
+if __name__ == "__main__":
+    main()

src/charts/ChartManager.js

@@ -63,6 +63,7 @@ import ErrorComponentReactWrapper from "@/react/components/ErrorComponentReactWr
 import ViewDialogWrapper from "./dialogs/ViewDialogWrapper";
 import { deserialiseParam, getConcreteFieldNames } from "./chartConfigUtils";
 import AddChartDialogReact from "./dialogs/AddChartDialogReact";
+import DGEDialogReact from "./dialogs/DGEDialogReact";
 import MenuBarWrapper from "@/react/components/MenuBarComponent";
 
 
@@ -1582,6 +1583,17 @@ export class ChartManager {
             new BaseDialog.experiment["AnnotationDialogReact"](ds.dataStore);
         });
 
+        if (dataStore.links) {
+            const hasRac = Object.values(dataStore.links).some(
+                (linkSet) => linkSet.rows_as_columns,
+            );
+            if (hasRac) {
+                this.addMenuIcon(ds.name, "fas fa-dna", "Run DGE Analysis", () => {
+                    new DGEDialogReact(dataStore);
+                });
+            }
+        }
+
         const idiv = createEl(
             "div",
             {