Add attribution analyses, reproducibility metadata, and CI

Author: shraddhapiparia

.github/workflows/ci.yml

@@ -0,0 +1,29 @@
+name: CI
+
+on:
+  push:
+  pull_request:
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    defaults:
+      run:
+        shell: bash -l {0}
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: mamba-org/setup-micromamba@v2
+        with:
+          micromamba-version: latest
+          environment-file: environment.yml
+          init-shell: bash
+          cache-environment: true
+          create-args: cpuonly   # CPU-only PyTorch; avoids large CUDA download on runners
+
+      - name: Smoke test — Phase 1 → Phase 2 on synthetic data
+        run: bash test_run.sh
+
+      - name: Unit tests
+        run: pytest -q

.gitignore

@@ -5,8 +5,7 @@ geno_ld/
 metadata/
 
 # Un-ignore the synthetic smoke-test generator (script only, not generated files).
-# data/* is used (not data/) so git can traverse into data/synthetic/ for exceptions.
-data/
+data/*
 !data/synthetic/
 data/synthetic/region_blocks/
 data/synthetic/block_plan/
@@ -42,3 +41,4 @@ build/
 # IDE
 .vscode/
 .idea/
+.claude/

CLAUDE.md

@@ -1,5 +1,9 @@
 # Block-Based Genotype Embedding Analysis
 
+[![CI](https://github.com/shraddhapiparia/blockbased-genotype-embedding-analysis/actions/workflows/ci.yml/badge.svg)](https://github.com/shraddhapiparia/blockbased-genotype-embedding-analysis/actions/workflows/ci.yml)
+
+Status: ongoing research. README reflects current findings; documentation, CI, and tests are being added incrementally.
+
 Unsupervised learning of subject-level genomic representations from LD-aware blocks,
 with application to asthma-relevant loci and downstream phenotype association.
 
@@ -32,14 +36,19 @@ The result is an embedding that is biologically interpretable (attention scores)
   axis of subject-level variation.
 - **HLA block embeddings outperform ancestry PCs.** HLA block_PC1 explains subject
   cluster structure beyond what genotype PC1–10 can account for, confirming genuine
-  biological signal rather than ancestry confounding.
+  biological signal rather than ancestry confounding. HLA sb15 η² = 0.767 vs η² = 0.051
+  for genotype PC3 (script `02_subject_cluster_analysis.py`, step C;
+  `results/output_regions2/ORD/clustering/hla_block_cluster_association.tsv` and
+  `results/output_regions2/ORD/clustering/genotype_pc_cluster_association.tsv`).
 - **PDE4D emerges after masking HLA.** A leave-HLA-out re-clustering experiment
   (script `03_leave_hla_out_analysis.py`) reveals PDE4D as the next most structurally
   informative block — consistent with its established role in asthma and β-agonist
   pharmacogenomics.
 - **Phenotype signal is real but subtle; IgE is the strongest.** Continuous phenotypes
-  (blood eosinophil count, IgE, lung function spirometry test, exacerbation) show the most consistent association
-  with block-level PC features across subjects.
+  (blood eosinophil count, IgE, lung function spirometry test, exacerbation) show the
+  most consistent association with block-level PC features across subjects
+  (script `01_block_embedding_phenotype_analysis.py`;
+  `results/output_regions2/ORD/all_blocks_pheno_analysis/phenotype_block_associations.tsv`).
 - **Biology recovered without phenotype labels.** The model was trained unsupervised on
   genotype data only. The emergence of HLA class II and PDE4D in post-hoc analysis
   validates that the learned geometry reflects known asthma biology.
@@ -53,10 +62,12 @@ components reflect ancestry rather than disease-relevant biology. Phase 1 VAE em
 preserve local LD-block haplotype structure that SNP-level PCA discards. Phase 2 adds
 cross-block context: the Transformer learns which blocks co-vary meaningfully across
 subjects, reorganizing rather than destroying the Phase 1 geometry (Phase 1 vs Phase 2
-pairwise-distance correlation ≈ 0.68). The result is a subject-level space where HLA
-class II dominates the primary axis, PDE4D emerges as the next structurally informative
-block after HLA removal, and IgE shows stronger phenotype association — biological signal
-that ancestry-adjusted PCA does not recover.
+median pairwise-distance Spearman r ≈ 0.68 across blocks; script
+`06_phase1_phase2_block_comparison.py`,
+`results/phase_comparison/pdm_correlations.csv`). The result is a subject-level space
+where HLA class II dominates the primary axis, PDE4D emerges as the next structurally
+informative block after HLA removal, and IgE shows stronger phenotype association —
+biological signal that ancestry-adjusted PCA does not recover.
 
 ---
 
@@ -76,11 +87,15 @@ cross-block Transformer (Phase 2) to produce subject embeddings is shown above.
 
 ![Subject PCA colored by cluster](docs/images/subject_pca_clusters.png)
 
-PCA of Phase 2 subject embeddings reveals three reproducible strata (`k=3`; ARI = 0.999).
-The weak silhouette score (`0.139`) indicates the learned space is structured as a
-continuous gradient rather than sharply separated clinical subtypes.
+PCA of Phase 2 subject embeddings reveals three reproducible strata (`k=3`; mean
+pairwise ARI = 0.999 across 50 random seeds). The weak silhouette score (`0.139`)
+indicates the learned space is structured as a continuous gradient rather than sharply
+separated clinical subtypes.
 
-*Source: `scripts/analysis/02_subject_cluster_analysis.py`.
+*Source: ARI from `scripts/analysis/04_cluster_stability_analysis.py`
+(`results/output_regions2/ORD/cluster_stability/ari_kmeans_seeds.tsv`); silhouette from
+`scripts/analysis/02_subject_cluster_analysis.py` step A
+(`results/output_regions2/ORD/clustering/clustering_metrics.csv`).
 Filename: `docs/images/subject_pca_clusters.png`*
 
 ---
@@ -90,8 +105,12 @@ Filename: `docs/images/subject_pca_clusters.png`*
 ![HLA class II dominates the learned embedding space](docs/images/hla_embedding_dominance.png)
 
 HLA class II subblocks strongly organize the Phase 2 embedding space. HLA sb15 explains
-far more cluster variance than ancestry PCs (η² = 0.767 vs 0.051 for genotype PC3) and
-correlates strongly with the main embedding axis (EmbedPC1 r = −0.88).
+far more cluster variance than ancestry PCs (η² = 0.767 vs 0.051 for genotype PC3;
+script `02_subject_cluster_analysis.py` step C,
+`results/output_regions2/ORD/clustering/hla_block_cluster_association.tsv`). HLA sb15
+block_PC1 also shows strong negative Spearman correlation with the main embedding axis;
+the specific value (r ≈ −0.88) was observed in exploratory analysis and is not yet
+tracked in a dedicated output table.
 *Source: `scripts/analysis/02_subject_cluster_analysis.py`.*
 
 ---
@@ -108,27 +127,62 @@ biologically coherent without supervised training.
 
 ---
 
+## Limitations
+
+- **Internal validation only.** All results are from a single cohort (COS/TRIO).
+  Clustering structure, attention rankings, and phenotype associations have not been
+  validated in an independent dataset.
+- **No causal inference.** Attention weights and LOBO/perturbation attribution scores
+  identify blocks that are statistically influential for the learned embedding geometry.
+  They are heuristic measures, not evidence of mechanistic importance or causality.
+- **Post-hoc phenotype associations.** Associations between embedding features and
+  clinical phenotypes (IgE, eosinophils, FEV1, exacerbation) were tested after
+  unsupervised training. They should be treated as hypothesis-generating and require
+  external validation.
+- **Ancestry and generalizability.** The embedding was trained within a single ancestry
+  stratum. Generalization to other ancestries or multi-ancestry cohorts has not been
+  assessed.
+- **Data availability.** Input genotype and phenotype data are access-restricted and not
+  version-controlled. Full reproduction from this repository requires obtaining restricted
+  data access separately. The synthetic smoke test (`bash test_run.sh`) validates pipeline
+  wiring only.
+- **Embedding geometry is not clinical classification.** The learned subject-level
+  embedding space reflects genomic variation structure, not disease phenotype. Cluster
+  membership does not imply clinical disease subtypes or therapeutic relevance.
+
+---
+
 ## Repository Structure
 
 ```
 scripts/
   core/       Core pipeline — Phase 1 VAE, Phase 2 Transformer, block analysis, plotting
-  analysis/   Numbered post-hoc scripts (01–07): phenotype association, clustering,
-              HLA validation, confounder analysis, 17q21 validation
+  analysis/   Numbered post-hoc scripts:
+                01–07  phenotype association, clustering, HLA validation,
+                       confounder analysis, 17q21 validation
+                08     clinical PC / embedding alignment (Pearson + Spearman, OLS, Ridge)
+                09     Phase 2 block attribution — leave-one-block-out (LOBO)
+                10     Phase 1 SNP attribution within selected blocks
   archive/    Superseded wrappers, exploratory one-offs, debug scripts
-configs/      YAML configs for Phase 1, Phase 2, and no-HLA variant
+configs/      YAML configs for Phase 1, Phase 2, no-HLA variant, and synthetic smoke test
 data/         Genotype block files and block manifest (access-restricted, not tracked)
+              data/synthetic/  — generated synthetic data for smoke testing (gitignored)
 results/      Pipeline outputs (access-restricted, not tracked)
 metadata/     Phenotype table, eigenvec file (access-restricted, not tracked)
 docs/         Method notes and figures
+tests/        Lightweight unit tests (pytest); no real data required
 environment.yml
+pytest.ini
 WORKFLOW.md   Step-by-step execution guide with CLI examples
 run_pipeline.sh  Single entry point — runs full pipeline or --dry-run input check
+test_run.sh   Smoke test — Phase 1 → Phase 2 on synthetic data (no restricted data needed)
 CLAUDE.md     AI assistance constraints and workflow summary
 ```
 
 See [WORKFLOW.md](WORKFLOW.md) for full CLI instructions, expected inputs/outputs per
-step, and execution order.
+step, and execution order. See [docs/data_contract.md](docs/data_contract.md) for the
+Phase 1 → Phase 2 data format contract, subject/block ordering conventions, and common
+failure modes.
 
 ---
 
@@ -150,6 +204,12 @@ python scripts/core/attention_phase2.py --config configs/config_phase2.yaml
 
 # Post-hoc analysis (example)
 python scripts/analysis/03_leave_hla_out_analysis.py
+
+# Attribution pipeline (scripts 08–10; requires Phase 2 outputs)
+python scripts/analysis/08_clinical_pc_embedding_alignment.py
+python scripts/analysis/09_phase2_block_attribution.py
+python scripts/analysis/10_phase1_snp_attribution_within_blocks.py \
+    --selected-blocks region_9p24_IL33
 ```
 
 Use `--dry-run` on `run_pipeline.sh` or on either phase script to validate inputs without

README.md

@@ -88,6 +88,34 @@ Scripts are numbered in intended execution order.
 
 ---
 
+## Attribution scripts (08–10)
+
+Run after Phase 2 and the core analysis scripts. Scripts 08–10 form a
+hierarchical attribution pipeline: embedding alignment → block attribution →
+SNP attribution.
+
+| Script | Purpose | Key inputs | Key outputs |
+|---|---|---|---|
+| `08_clinical_pc_embedding_alignment.py` | Align learned Phase 2 embedding PCs with clinical phenotype PCs; Pearson + Spearman heatmaps; OLS and Ridge baselines including ancestry-PC-only model | Phase 2 `individual_embeddings.csv`, `ldpruned_997subs.eigenvec`, phenotype CSV | `correlation_heatmap_pearson.png`, `correlation_heatmap_spearman.png`, `embedding_pc_correlations.tsv`, R² comparison tables |
+| `09_phase2_block_attribution.py` | Leave-one-block-out (LOBO) attribution: for each of the 174 Phase 1 blocks, mean-mask it in the Phase 2 input and measure change in embedding PC1/PC2 and IgE ridge score | Phase 1 `all_blocks.npy`, Phase 2 checkpoint + `individual_embeddings.npy`, `block_order.csv`, `subjects.csv` | `phase2_PC1_leave_one_block_out.csv`, `phase2_PC2_leave_one_block_out.csv`, `phase2_log10Ige_leave_one_block_out.csv`, top-20 bar plots |
+| `10_phase1_snp_attribution_within_blocks.py` | For user-selected blocks from script 09: mean-mask each SNP in the Phase 1 input, re-encode through the frozen VAE, re-run Phase 2, project onto fixed PCA axes, and rank SNPs by mean absolute delta | Phase 1 block checkpoints (`.pt`), per-block `.raw` genotype files, Phase 2 checkpoint + embeddings, script 09 LOBO CSVs (optional, enables block-weighted priority) | Per-block and combined `*_snp_attribution.csv`, bar plots |
+
+```bash
+# Run 08 (requires Phase 2 outputs at default paths)
+python scripts/analysis/08_clinical_pc_embedding_alignment.py
+
+# Run 09
+python scripts/analysis/09_phase2_block_attribution.py
+
+# Run 10 (requires --selected-blocks from 09 results)
+python scripts/analysis/10_phase1_snp_attribution_within_blocks.py \
+    --selected-blocks region_9p24_IL33,region_6p21_HLA_classII_sb15 \
+    --lobo-pc1-csv results/analysis/phase2_block_attribution/phase2_PC1_leave_one_block_out.csv \
+    --lobo-pc2-csv results/analysis/phase2_block_attribution/phase2_PC2_leave_one_block_out.csv
+```
+
+---
+
 ## Archive
 
 `scripts/archive/` contains scripts that are **not part of the active workflow**: