SEPIA — current state & handoff

SEPIA — current state & handoff

Living snapshot, last updated 2026-07-24. Its job: let anyone (or a fresh, post-compaction session) pick up cold without losing decisions or context. Companions: changelog.md = dated log of what was done/found/decided; plan.md = the forward roadmap; this file = the current-state + open-questions snapshot. Update this whenever a decision is made or a track finishes.

Orientation — where everything lives

• GitHub repo: lasisilab/sepia (renamed from paint-workflowr). Branch: main only. This is the published anchor; Lily Heald's PAINT thesis is the prior work (her own repo is separate + untouched).
• Laptop checkout: ~/GitHub/sepia — note the on-disk folder may still literally be named paint-workflowr (a cosmetic rename is pending; git + remote are already sepia).
• Cluster checkout (repo = workdir): /nfs/turbo/lsa-tlasisi1/sepia on UMich Great Lakes. It can now push (SSH remote + an SSH key SSO-authorized for lasisilab).
• Shared lab genomes (outside the project, symlinked in as resources/genomes; reusable by other lab studies): /nfs/turbo/lsa-tlasisi1/genomes/ → reference/ (hs37d5, hg19), SGDP/ (completed 159-sample reference), archaic/{vcf,mask,bam}/ + MANIFEST.tsv.
• Lily's original data (READ-ONLY seed; do not modify/delete — Tina will clean it later): /nfs/turbo/lsa-tlasisi1/lheald_thesis/aDNA_data/.
• Git tracking model: big data is gitignored (resources/, scratch/, *.bam/.bcf/.vcf.gz/.fa/.bed/.bim/.fam…); small results sync via output/.

Decisions locked in (with the why)

1. Project = SEPIA (Systematic Evaluation of Pigmentation In Archaic hominins); a QC-first re-analysis continuing PAINT.
2. Build strategy — report/integrate on hg38; genotypes physically stay hg19; bridge by a verified rsID↔︎hg19↔︎hg38 map. A genotype is build-independent, so this is lossless relabeling — no re-alignment. The genome-wide sanity PCA stays hg19 (a control). Genome-wide hg38 (B2b) is optional/deferred. (This reconciled an earlier hg38-monolithic drift and a mid-session hg19 over-correction.)
   • Replication path is hg19-native (confirmed 07-24): step one (undo A2, redo Lily's pigmentation figures) = extract the 222 panel at the correct hg19 coords (pos_hg19 in the map, all 222 present) from the hg19 data — cheapest correct fix, no hg38 needed. hg38 is only the reporting/integration layer for the later GWAS-superset + HIrisPlex-S work; the map carries both builds so nothing is duplicated.
3. Site build: migrated workflowr → Quarto (reconciled onto main); auto-deploy paused (workflow_dispatch) until the R-analysis pages are rendered (they need data / will be regenerated by SEPIA).
4. Trust: rebuild from raw/published sources with QC gates; do not trust Lily's derived data.
5. Sample subset (papers/SAMPLE_INCLUSION.md): high-cov (Altai/Den5, Vindija 33.19, Denisova 3, Chagyrskaya 8, Denisova 25) → published snpAD VCFs, no BAM; Mez1 → its VCF; low-cov Les Cottés, Denisova 11, Goyet → BAMs (ANGSD GLs + mapDamage); excluded: El Sidrón (exome), Vindija 87 (= Vi 33.19), Mez2, Spy, HST, Scladina.
6. Modern reference: complete the 166-SGDP merge → 159 usable (7 were truncated local copies, dropped). Sanity 1 PASSED — recovers textbook global structure.
7. Shared genomes = full canonical genomes (reusable), + cluster→GitHub push via SSH + SAML-SSO.
8. Merge architecture lesson: NFS punishes random seeks — stream-thin (bcftools view -T) in parallel, don't merge -R over scattered regions.

Current state

• DONE: SEPIA rename+rebrand; Quarto reconciled (deploy gated); cluster reorg + shared genomes; SGDP 159-sample reference + Sanity-1 PCA (output/pca_wg/); SGDP metadata (data/sgdp_metadata.tsv, Mallick 2016/Reich); archaic high-cov VCFs (5 genomes, ~258 GB) + Den25 masks; cluster push fixed.
• DONE → data/panel_hg19_hg38_map.tsv (+ .SOURCE.md) — the B2a rsID↔︎hg19↔︎hg38 map, 222/222 resolved (Ensembl REST + dbSNP), per-build ref/alt; 11 palindromic flagged (Q4), 9 merged rsIDs. FOUND a uniform +1 off-by-one in the panel's hg38 coords (snps.txt − true GRCh38 = +1 across all 222; a 0/1-based error). True tracers: SLC24A5 rs1426654 = 15:48134287 (panel stored 48134288), HERC2 rs12913832 = 15:28120472. hg19 was correct → genome-wide/SGDP work unaffected. The map re-pulls authoritative coords, so it corrects this. Doc follow-up: fix the tracer definitions (glossary/Q2) 48134288→48134287; leave BUG_EVIDENCE/A2 (they quote the panel's buggy value on purpose).
• DONE (B2a — the pigmentation bug is fixed) — code/pig_panel_extract.slurm (job 54641658): 222 panel @ correct hg19 coords → JOINT PCA of 159 SGDP + 5 high-cov archaics, 160 informative SNPs. Eigenvalues 31/15/8.4/6.6/5.7… (PC1≈35%) — real structure, not rank-1. Archaics cluster near the African end. Output output/pca_pig/ (committed 8a1d5e2). Figure plan/figures/pig_pca_before_after.png → Fig 2 + A2 toggle on the plan page. (Lesson: archaic extraction uses -R seek, not -T stream — the 4 GB per-chr all-sites VCFs.)
• DONE (site live + current) — Quarto deploy un-stuck: project.render restricted to CI-buildable pages (index/about/license + frozen sgdp_reference_pca); R analysis pages dropped from navbar until a data-bearing local render populates _freeze/. Live: https://lasisilab.github.io/sepia/ (roadmap /plan/sepia-plan.html). Redeploy is manual (gh workflow run publish-site.yml) — push alone does NOT update the site.
• DONE (figure) → plan/figures/sgdp_pca_before_after.png — SGDP WG PCA, Lily's 15 vs SEPIA's 159, coloured by region (Fig 0). "For Lily" section now teaches A2/A1/A3 + fixes.
• DONE (B1 — gene-network integration) → data/pigment_snp_map.tsv — the extended 1,155-SNP map (built by code/build_pigment_snp_map.py): union of GWAS Catalog + HIrisPlex-S + Martin2017 + Kim2024 + discordance, re-resolved to verified hg38+hg19 (tracers exact). Strict superset of the 222 panel (222 in + 933 new). Confidence: 453 high / 60 mod / 642 low. Documented + LIVE at snp_panel.qmd → /snp_panel.html. Inventory: papers/GENE_NETWORK_DATASETS.md.
• DONE (B3 — wider pigmentation PCA) — code/wider_pca.slurm (job 54657946): extended set, all ~923 SNPs + high-confidence-only (~350), joint 159 SGDP + 5 archaics. Fig 4 on the plan page. output/pca_wider/.
• DONE (archaic PROJECTION, WG + pigmentation = Sanity 2) — installed newer plink2 (tools/plink2, v2.0.0-a.7, gitignored; the ARC module is 2019 + can't emit loadings). Modern-PCA-then-project-archaics for both spaces: code/project_archaics.slurm → output/projections/, validated r = 1. Workbook archaic_projection.qmd (live, navbar "Projection").
• DONE (B6/B7 — gene-region PCA) — pigmentation-gene map data/pigment_gene_map.tsv (1,004 genes → hg19+hg38 regions, Ensembl) + code/gene_region_pca.slurm (all variants in gene bodies ±10 kb, LD-pruned) → output/pca_generegion/.
• ★ KEY FINDING — whole-genome and pigmentation-gene-region PCAs are near-identical (archaics = divergent outgroup at the non-African edge, high-cov PC1 ≈ 0.15); only the ascertained pigmentation-SNP panel places archaics with Africans. ⇒ genes = ancestry, ascertained SNPs = selection: "archaics with Africans" is the footprint of selection at derived light-skin loci, not a gene-ancestry signal. Three-way figure plan/figures/three_way_projection.png; full logic in the workbook.
• IN FLIGHT: nothing running — all cluster jobs COMPLETED (SGDP merge, B2a, B3, projection, gene-region, acquisition).
• DATA NOW AVAILABLE: acquisition COMPLETED — low-cov BAMs landed (12 files, ~29 GB, genomes/archaic/bam/) → ANGSD GLs + mapDamage (Q7/Q8) unblocked.
• KNOWN GAPS (small): (a) Chagyrskaya = 0 records in the panel/gene extraction from the high-cov snpAD VCF dir (.noRB files) → B2a/B3 have 5 not 6 archaics (the projection uses ancient_sgdp_wg, which DOES include Chagyrskaya); (b) LP6005441-DNA_A09 missing from sgdp_metadata.tsv (shows "unknown"); (c) acquire_archaic.slurm built malformed https://ERR… URLs for some BAM accessions (still got 12).

Open questions / awaiting Tina

• Full 166 SGDP? The 7 dropped are truncated local copies; a clean re-acquire needs the CGC all-sites data (Seven Bridges login) — the public Reich route is variant-only and would bias them. Default: stay at 159 (ample) unless you want the CGC pull.
• Quarto site go-live timing — the R-analysis pages need their data (laptop-only) rendered, or SEPIA regenerates them; also git rm docs/ (vestigial workflowr) + wire the introgression_app into Quarto before deploy.
• Gene-network integration scope — DECIDED (TL, 07-24): fold in everything SNP-level (done → the 1,155-SNP map), do the wider PCA first (B3, running), tag high/low confidence (done). Defer the gene-level tables (Baxter/Bajpai/Raghunath/D'Arcy → B6/B7). Next open call: HIrisPlex-S phenotype prediction (B4) needs the numeric weights pulled from Walsh 2017 / the Erasmus tool.
• Local folder rename (~/GitHub/paint-workflowr → sepia) — cosmetic, do at a break (git works regardless).

Immediate next steps

1. B4 — HIrisPlex-S phenotype prediction (the "what did they look like" question the projection can't answer): pull the numeric weights (Walsh 2017 / Erasmus tool), score the archaics at the HIrisPlex-S markers. (This is the natural next scientific step.)
2. ANGSD GLs + mapDamage (Q7/Q8) on the low-cov BAMs (now present in genomes/archaic/bam/) — genotype likelihoods + aDNA-damage authenticity.
3. f/D-statistics to quantify the genome-wide archaic position / Neanderthal-admixture lean (PCA is too coarse for it).
4. Small fixes: Chagyrskaya .noRB panel extraction; add LP6005441-DNA_A09 to sgdp_metadata.tsv; fix the acquisition BAM-URL resolution.
5. R analysis pages → a data-bearing local quarto render to populate _freeze/, then restore introduction/depth/pca to the navbar.

Handoff note (for a fresh Claude / Claude Science)

Read this file + changelog.md top-to-bottom, then plan.md. The site's step-by-step workbooks are the clearest narrative: sgdp_reference_pca.qmd (Sanity 1), snp_panel.qmd (the SNP set), archaic_projection.qmd (projection + the key gene-vs-SNP finding). The projection recipe + its hard-won gotchas are documented in code/project_archaics.slurm (and reused in code/gene_region_pca.slurm). Cluster access needs Tina to hold the SSH tunnel open (Kerberos+Duo) — Claude cannot log in itself; ssh -O check greatlakes first, always. Newer plink2 lives at tools/plink2 (gitignored). Coordinate maps (data/pigment_snp_map.tsv, data/pigment_gene_map.tsv) are the reusable substrates.

Cluster-access rules (still binding)

ssh -O check greatlakes before batched commands; one batched command at a time; never hammer; heavy work via sbatch, not the login node. Full rules in CLUSTER_ACCESS.md.