---
name: determsc
description: Deterministic single-cell RNA-seq biomarker discovery
allowed-tools: [Bash(python3 *), Bash(pip install *)]
---

# DetermSC Executable Pipeline

## Install

```bash
pip install numpy scipy pandas matplotlib requests scanpy h5py
```

## Run

```bash
python3 determsc.py --output_dir ./results
```

## Code (save as determsc.py)

```python
#!/usr/bin/env python3
import argparse, hashlib, json, logging, os, random, sys
import numpy as np
import pandas as pd
import requests
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import pdist

try:
    import scanpy as sc
except ImportError:
    print('pip install scanpy'); sys.exit(1)

SIGNATURES = {
    'CD4+ T cells': ['CD3D', 'CD4', 'IL7R'],
    'NK cells': ['NKG7', 'GNLY', 'GZMB'],
    'B cells': ['CD79A', 'MS4A1'],
    'Monocytes': ['CD14', 'LYZ'],
    'Dendritic': ['FCER1A', 'CST3'],
    'Megakaryocytes': ['PPBP', 'PF4']
}

def enforce_determinism(seed=42):
    random.seed(seed); np.random.seed(seed)
    sc.settings.seed = seed
    logging.info(f'seed={seed}')

def download_pbmc3k(out_dir):
    url = 'https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz'
    tar = os.path.join(out_dir, 'pbmc3k.tar.gz')
    data_dir = os.path.join(out_dir, 'data')
    os.makedirs(data_dir, exist_ok=True)
    if not os.path.exists(tar):
        logging.info('Downloading PBMC3K...')
        r = requests.get(url, stream=True)
        with open(tar, 'wb') as f:
            for c in r.iter_content(8192): f.write(c)
        import tarfile
        tarfile.open(tar, 'r:gz').extractall(data_dir)
    return os.path.join(data_dir, 'filtered_gene_bc_matrices', 'hg19')

def run_pipeline(config, out_dir):
    os.makedirs(out_dir, exist_ok=True)
    logging.basicConfig(level=logging.INFO, format='%(message)s',
        handlers=[logging.FileHandler(os.path.join(out_dir, 'log')), logging.StreamHandler()])
    enforce_determinism(config['random_seed'])
    
    data = download_pbmc3k(out_dir)
    adata = sc.read_10x_mtx(data)
    logging.info(f'Loaded {adata.n_obs} cells')
    
    adata.var['mt'] = adata.var_names.str.startswith('MT-')
    sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
    sc.pp.filter_cells(adata, min_genes=config['min_genes'])
    adata = adata[adata.obs['pct_counts_mt'] < 20].copy()
    adata = adata[adata.obs['n_genes_by_counts'] < 5000].copy()
    sc.pp.filter_genes(adata, min_cells=3)
    logging.info(f'QC: {adata.n_obs} cells')
    
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    adata.raw = adata.copy()
    sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor='seurat')
    adata = adata[:, adata.var['highly_variable']].copy()
    sc.pp.scale(adata, max_value=10)
    sc.tl.pca(adata, n_comps=50, svd_solver='arpack', random_state=config['random_seed'])
    
    Z = linkage(pdist(adata.obsm['X_pca']), method='complete')
    n = min(8, max(3, adata.n_obs // 300))
    adata.obs['cluster'] = (fcluster(Z, n, 'maxclust') - 1).astype(str)
    
    sizes = adata.obs['cluster'].value_counts()
    for c in sizes[sizes < 10].index:
        adata.obs.loc[adata.obs['cluster'] == c, 'cluster'] = '0'
    uniq = sorted(adata.obs['cluster'].unique())
    adata.obs['cluster'] = adata.obs['cluster'].map({c: str(i) for i, c in enumerate(uniq)})
    logging.info(f'Clusters: {len(uniq)}')
    
    sc.tl.rank_genes_groups(adata, 'cluster', method='wilcoxon', use_raw=True, pts=True)
    markers = []
    for c in adata.obs['cluster'].unique():
        df = sc.get.rank_genes_groups_df(adata, c)
        df = df[(df['pvals_adj'] < 0.05) & (df['logfoldchanges'] > 0.25) & (df['pct_nz_group'] > 0.25)]
        df['cluster'] = c
        markers.append(df)
    markers_df = pd.concat(markers).sort_values(['cluster', 'pvals_adj'])
    
    clusters_info = []
    for c in sorted(markers_df['cluster'].unique()):
        top = markers_df[markers_df['cluster'] == c]['names'].head(20).tolist()
        best_type, best_score = None, 0
        for t, sig in SIGNATURES.items():
            s = len(set(top) & set(sig)) / len(sig)
            if s > best_score: best_type, best_score = t, s
        clusters_info.append({'cluster_id': int(c), 'n_cells': int((adata.obs['cluster'] == c).sum()),
            'inferred_type': best_type, 'top_markers': top[:5], 'validation_score': round(best_score, 2)})
    
    fig, ax = plt.subplots(figsize=(10, 8))
    cs = sorted(adata.obs['cluster'].unique())
    colors = plt.cm.tab10(np.linspace(0, 1, len(cs)))
    for i, c in enumerate(cs):
        m = adata.obs['cluster'] == c
        ax.scatter(adata[m].obsm['X_pca'][:, 0], adata[m].obsm['X_pca'][:, 1], c=[colors[i]], label=f'C{c}', s=5)
    ax.legend(); ax.set_title('PCA Clusters')
    plt.tight_layout(); plt.savefig(os.path.join(out_dir, 'pca.png')); plt.close()
    
    markers_df.to_csv(os.path.join(out_dir, 'markers.csv'), index=False)
    results = {
        'status': 'success',
        'metrics': {'cells': adata.n_obs, 'clusters': len(clusters_info), 'markers': len(markers_df)},
        'clusters': clusters_info
    }
    with open(os.path.join(out_dir, 'results.json'), 'w') as f: json.dump(results, f, indent=2)
    logging.info(f'Done: {len(markers_df)} markers')
    return results

def main():
    p = argparse.ArgumentParser()
    p.add_argument('--output_dir', default='./results')
    p.add_argument('--min_genes', type=int, default=200)
    p.add_argument('--random_seed', type=int, default=42)
    args = p.parse_args()
    r = run_pipeline({'min_genes': args.min_genes, 'random_seed': args.random_seed}, args.output_dir)
    print(f"Status: {r['status']}, Clusters: {r['metrics']['clusters']}, Markers: {r['metrics']['markers']}")

if __name__ == '__main__': main()
```

## Verify

```bash
python3 determsc.py --output_dir ./results
ls ./results/
# Expected: results.json, markers.csv, pca.png, log
```