Integrating multi-tissue transcriptomic signatures from Multiple Sclerosis studies with subgraph isomorphism algorithms will enable the identification of conserved molecular interaction motifs related to disease progression.

Subgraph isomorphism algorithms applied to multi-tissue transcriptomic co-expression networks derived from ≥3 independent MS datasets (GSE193770, GSE108000, GSE138614) will identify conserved molecular interaction motifs — defined as recurring subgraph patterns of ≥4 nodes with edge-weight Pearson r ≥ 0.6 — that are (a) statistically enriched at chronic active rim (CA-RIM) lesion tissue relative to normal-appearing white matter (NAWR) and control tissue (FDR < 0.05, enrichment ratio ≥ 1.5×), (b) reproducible across ≥2 independent cohorts, and (c) contain ≥1 druggable node from the validated target hierarchy (DNMT1, ZNF740/BRD3, CTSS) with composite validation score ≥ 0.578. The null hypothesis is that subgraph isomorphism detects no motifs beyond those recoverable by standard pairwise DEG overlap (Jaccard similarity ≤ 0.15 improvement over baseline).

MINIMUM VIABLE TEST (MVT) — 3-phase design targeting 45-day completion:

PHASE A — Network Construction & Motif Discovery (Days 1–20): Construct tissue-stratified co-expression networks from three GEO datasets. For bulk datasets (GSE108000, GSE138614): VST-normalize counts, compute FDR-corrected partial correlations (PCIT or GeneNet R package), threshold at r ≥ 0.6, build adjacency matrices. For single-cell dataset (GSE193770): use scVI latent space correlations from the pre-built atlas (gs://aegismind-tpu-results/ms_phase2/results/) to construct cell-type-stratified networks, focusing on CD8+ T-cell cluster (Leiden cluster identity to be confirmed from atlas). Apply VF2++ subgraph isomorphism algorithm (NetworkX implementation, GPU-accelerated via cuGraph) to enumerate recurring subgraphs of size k = 4, 5, 6, 7, 8 nodes. Use gSpan or GRAMI for frequent subgraph mining with minimum support threshold = 2/3 datasets.

PHASE B — Statistical Validation & Enrichment (Days 21–35): Compute motif enrichment scores versus 10,000 permuted networks (node-label permutation preserving degree distribution). Apply Bonferroni correction. Map motif hub nodes to: (i) the 50-gene MS seed set for proximity scoring, (ii) ChEMBL druggability database, (iii) the validated target hierarchy. Perform GO/KEGG enrichment on motif node sets (clusterProfiler). Compute Jaccard similarity of motif node sets between cohort pairs.

PHASE C — Biological Validation (Days 36–45): For top 3 conserved motifs: validate hub node expression in CD8+ T cells using FACS-sorted PBMC from 10 smoldering MS patients versus 10 healthy controls (if samples available) OR re-analyze existing sorted data from GSE193770. Test whether motif hub nodes predict CA-RIM pathology score using logistic regression (AUC as primary metric). Generate network visualization and report.

# ============================================================
# EVP IMPLEMENTATION SKETCH: Subgraph Isomorphism MS Motif Discovery
# Target runtime: 45 days | GPU: 85h | CPU: 320h | RAM: 512GB
# ============================================================

# --- PHASE A: DATA LOADING & NETWORK CONSTRUCTION ---

import scanpy as sc
import scvi
import networkx as nx
import numpy as np
import pandas as pd
from scipy.stats import pearsonr, spearmanr
from gspan_mining import gSpan  # pip install gspan-mining
import cugraph  # GPU-accelerated graph ops
import anndata

# Step 1: Load pre-built scVI atlas (saves ~2 weeks retraining)
atlas = sc.read_h5ad("gs://aegismind-tpu-results/ms_phase2/results/atlas.h5ad")
# Confirm CD8+ T-cell cluster
cd8_mask = atlas.obs['leiden'].isin(identify_cd8_clusters(atlas))
cd8_adata = atlas[cd8_mask].copy()
# Expected: ~4,000–8,000 CD8+ T cells from 36,966 total

# Step 2: Load bulk datasets
bulk_datasets = {
    'GSE108000': load_geo_bulk('GSE108000', normalize='VST'),
    'GSE138614': load_geo_bulk('GSE138614', normalize='VST')
}

# Step 3: Define gene universe (DEG union + MS seed set)
deg_genes = load_deg_list('carrim_degs_1065.txt')  # from GitHub repo
ms_seed = load_seed_set('ms_seed_50genes.txt')
gene_universe = list(set(deg_genes) | set(ms_seed))  # ~1,100 genes

# Step 4: Calibration checkpoint — verify published DEG values
def calibrate_degs(adata, expected_values):
    """Verify DNMT1, CTSS, ZNF740 log2FC within ±10% of published"""
    results = run_deseq2(adata, contrast=['condition', 'CA-RIM', 'NAWM'])
    for gene, expected_fc in expected_values.items():
        observed_fc = results.loc[gene, 'log2FoldChange']
        assert abs(observed_fc - expected_fc) / abs(expected_fc) < 0.10, \
            f"CALIBRATION FAIL: {gene} FC={observed_fc:.3f}, expected={expected_fc:.3f}"
    return results

expected = {'DNMT1': 1.59, 'CTSS': 1.16, 'ZNF740': 1.15}
calibrate_degs(bulk_datasets['GSE108000'], expected)  # ABORT if fails

# --- NETWORK CONSTRUCTION ---

def build_coexpression_network(expr_matrix, genes, method='partial_corr',
                                threshold=0.6):
    """
    Build co-expression network with FDR-corrected partial correlations.
    Returns: NetworkX Graph with edge weights
    """
    expr_subset = expr_matrix[genes].dropna(axis=1)
    
    if method == 'partial_corr':
        # Use R GeneNet via rpy2 for FDR-corrected partial correlations
        import rpy2.robjects as ro
        pcor_matrix = ro.r(f'''
            library(GeneNet)
            data <- as.matrix(read.csv("expr_temp.csv"))
            pcor <- ggm.estimate.pcor(data)
            test.results <- network.test