Reprocess snm3C-seq with new samples

[adamklie]

Context

For igvf-data/igvf_sc-islet_snm3c, reprocess with new samples.

Tasks

• Transfer new data from TSCC to cluster
• Run through pipeline (bin/1_process) and sample QC (bin/2_sample_qc)
• Deal with batch correction
• Continue downstream analysis steps

[adamklie]

Dependencies

• Feeds into Integrative analysis: snm3C-seq × Dominguez2020 bulk methylation #12 (snm3C × Dominguez2020 integration), Retrain CpGNet methylation models on new pseudobulk data #14 (CpGNet retraining)

[adamklie]

Notes from project reorganization (2026-03-23)

Detailed processing and integration documentation filed in igvf-data/igvf_sc-islet_snm3c/README.md:

Data integration approach

• Modality preprocessing: snmCG gene body methylation (55,486 features), snRNA, snATAC 5kb bins, snATAC gene activity, snmCG 5kb hypomethylation
• Joint embedding via CCA across modalities
• Consensus cell type annotation by majority label across 3 integration strategies

scHiCluster pipeline status

Imputation, gene scoring, compartments (single-cell), pseudobulking, loop calling, and embedding are done. Domains failed (R issue). Differential analysis still TODO.

Outstanding

• Differential methylation: pseudobulk DMG (t-test + ALLCools), single-cell DMG, DMRs (MethSCAN, methylpy)
• Figure 5 plans: differential loops (5a), compartments (5b), boundaries (5c), gene scores (5d)
• mCG pseudobulking for celltype+treatment: ALLC + bigwig generation done

[adamklie]

Status Update (2026-04-09)

Unblocked but deferred — this has many steps and we're prioritizing multiome analyses (#1, #4, #8) first. Will come back to this.

[adamklie]

Status Update (2026-04-16)

Scope has grown — now blocking #2 on HiC normalization

Sara (Gaulton lab) is re-engaging on peak-to-gene linking (#2) and specifically needs normalized HiC contact maps that should come out of this reprocessing. Her current results use default gamma/scale params on unnormalized maps and are degrading output quality. This promotes #11 from "deferred, many steps" to partial blocker for #2.

Audit in flight

Kicked off a background audit of igvf-data/igvf_sc-islet_snm3c/ to report:

• Current processing state + sample counts
• Current HiC contact-map state (what's there, from which pipeline, normalized or not)
• What it takes for full reprocessing (per the original issue scope)
• Whether there's a faster HiC-only renormalization path to unblock Sara without doing the whole reprocessing

Will fold findings into a follow-up comment.

Priority

Bumping up based on the #2 dependency. Not the top priority this week (that's #1 motif discovery completion + #8 ATAC integration + #6 GCP layout), but jumping ahead of the integrative analyses (#9, #10, #12, #14).

[adamklie]

Audit results — snm3C-seq state (2026-04-16)

Background audit returned.

What's there

• Raw: 114 sample directories in fastq/ (1,546 FASTQ files). Three differentiation batches × timepoints (6/24/48/72 h) × six conditions (control, 3-cyt, IFNγ, Ex-4_HG, palmitate, dex).
• Pipeline scaffolding: bin/1_get_data, bin/2_process_data, bin/3_methscan, bin/4_allcools, bin/5_hicluster, bin/6_pseudobulk, bin/7_cpgnet, bin/8_upload, bin/9_submission. Some steps present but not run against new samples.
• Current v1 YAP-based outputs (ca. Jan 2025): 3 pseudobulk HiC files (SC.alpha, SC.beta, SC.EC — 5.8–6.9 GB each in .hic format, generated via juicer_tools pre). Plus imputed contacts (25K/100K resolution), loop calls, compartments, gene scores, CpG ratios.

HiC — current normalization state

• .hic files produced with juicer_tools pre and NO explicit normalization flag — no -k option was passed. This matches what Sara is reporting: the files don't have stored normalization vectors, which is why she's been using ad-hoc default gamma/scale params.
• No mcool / cooler-format alternatives present.
• No renormalization tooling is present in the bin/ dirs either.

Gap analysis — full reprocessing (original issue scope)

Step	Status
1_process (alignment / demux)	❌ absent — v1 relied on YAP upstream; IGVF-standard snm3C alignment not implemented
2_sample_qc	❌ not implemented
3_methscan / 4_allcools	✅ code exists; not run on new samples
5_hicluster (contact maps)	✅ pseudobulk .hic exist for v1
Batch correction	❌ no code
Downstream pseudobulk / gene scoring	✅ partial for v1

Gap analysis — HiC-only quick-win for #2 unblock

• Option A: juicer_tools post -k KR on existing 3 .hic files → embeds KR normalization vectors in-place. ~5 min/file, ~15 min total. No reprocessing needed.
• Option B: convert to cooler (.hic → pairs → cooler), then cooler balance (ICE normalization). More portable but more setup.
• Caveat: no documentation of which normalization method is "correct" for snm3C pseudobulk. Needs a short conversation with Sara to pin down her preference before either approach.

Next concrete action

• Full reprocessing (medium-term): transfer remaining new samples from TSCC (script bin/1_get_data/1_transfer_fastqs.sh present but incomplete), decide on IGVF-standard snm3C alignment pipeline to adopt, fill in bin/2_sample_qc, run through 3–7 → final HiC.
• HiC-only quick-win (if we want to unblock Sara fast): ping Sara for preferred normalization method (KR / VC / SCALE / ICE), then run the chosen approach on the 3 existing .hic files and share. ~1 afternoon end-to-end.

Not pursuing the quick-win today — the method choice needs a conversation with Sara first.

[adamklie]

Audit + v0 archive committed (2026-04-28; commit `f227efc` in nested `igvf-data/igvf_sc-islet_snm3c/` repo).

Sample state — clearer picture

All 54 metadata samples are present in `processed/` already (no transfers needed). Mapping written: `results/1_get_data/sample_processed_mapping.tsv`.

differentiation_batch	downstream_status	n
JE002	`incorporated`	42
DM060	`sample_qc_only`	12

The 12 DM060 samples got sample_qc but were never merged downstream. Stage 4 merge / 3 methscan endo / 5 hicluster cell_table all show 0 DM060 cells. Original task framing of "transfer additional samples" was inaccurate — the samples are present; the gap is downstream incorporation.

Edge cases captured in the mapping file:

• 1 sample (`IGVF_44.12_2`) only has its `_deep` variant; no shallow. Use deep.
• 12 DM060 samples appear in both `processed/2024_06_06/` and `processed/2025_07_23/` (re-sequenced). Mapping uses 2025_07_23 as canonical.
• 5 orphan dirs in `processed/` not in the canonical 54-sample manifest (likely failed QC or replaced).

v0 archive (42-sample JE002 baseline)

Frozen ahead of the planned v1 re-run with all 54 incorporated:

• `archive/v0/bin/` — 17 MB code-only snapshot (98 files)
• `archive/v0/results/` — moved 42-sample outputs (3_methscan, 4_allcools/merge, 5_hicluster, 6_upload). `4_allcools/sample_qc/` stays at top level (54-sample state).
• `archive/v0/MANIFEST.md` describes baseline + what's NOT here.

Light cleanup (also committed)

Per audit findings:

• Stage 3/4/5 READMEs were identical copies of MethScan reference docs — rewrote each stage-specific
• Bespoke files moved to `scratch/2026_04_28/snm3c_v0_cleanup_archive/` (the `X_*`, exploratory DMR notebooks, dead-ref scripts)
• `bin/7_cpgnet/4_variant_effect.ipynb` (0 bytes) deleted
• Stage 5 hicluster duplicate numbering resolved (`4_make_cell_table.sh` → `4a_`, `4_embed.sh` → `4b_`)
• Stage 8 path bugs fixed in both upload notebooks (`results/6_upload/` → `results/8_upload/`)
• New top-level `bin/README.md`

Deferred for v1 (per "light" cleanup scope): `${DATASET_ROOT}` parameterization, 9 → 7 stage redesign.

Track A (the actual re-run) — still pending

Re-run from stage 4 merge through methscan / hicluster / pseudobulk / upload with all 54 samples incorporated. Multi-day cluster compute. Outputs land at `results/N_stage/` (overwriting v0 layout) and become v1.

Sara's HiC normalization concern (no `-k KR` in `bin/5_hicluster/scripts/hic.sh`) remains to be addressed in the v1 re-run. Audit confirms no normalization flag in `juicer_tools pre` invocation.

[adamklie]

v1 re-run plan written (commit `6861388`): `bin/v1_RERUN_PLAN.md`, 479 lines.

Stage-by-stage execution doc for incorporating the 12 DM060 samples downstream. Per-stage inputs, run commands, expected outputs, compute estimates, verification steps.

Compute estimates (from v0 timings):

• Stage 3 methscan: ~1 day pipeline + 1 day interactive annotation
• Stage 4 allcools merge + 3-res embed: ~1 day (256G / 2d slurm reservation)
• Stage 5 hicluster: ~3-7 days (impute is the dominant cost)
• Stage 6 pseudobulk: ~half a day
• Stage 8 upload: hours
• End-to-end: ~2 weeks

Open questions to resolve before kickoff:

1. HiC normalization (`-k KR` vs alternative) — Sara's call (carried from 04/16)
2. Cell-type annotation: re-cluster from scratch vs project v0 labels onto DM060 cells
3. Batch covariate handling (DM060 vs JE002) for any DE/DMR
4. DM060 condition coverage incomplete — only `control` + `dex`; cross-batch contrasts only valid for these
5. CGN-extraction prerequisite — `processed/2025_07_23/.../allc/` only has `.allc.tsv.gz`, not `.allc.CGN-Both.allc.tsv.gz`. Stage 2 may need a CGN extraction sub-step before stage 3 can run for DM060
6. v1 output paths: `results//v1/` subtree (recommended) vs direct overwrite

Smoke test proposed: `1_prepare.sh` on a single DM060 sample (~384 cells, IGVF_6D_12_1) before kicking off heavy stage-5 imputation.

Cluster execution will not happen without explicit user sign-off.

[TheShivaniLakkaraju]

here are my objects incl. dex and processed with Adam's methscan pipeline (just including 3C filters) on Adam's GCS bucket: they need celltype annotations. included my processing notebooks and scripts so you can verify too.

Parent: gs://2025_01_06/shivani_snm3c_methscan/

Outputs — filtered methscan matrices (41 GB, incl. dexamethasone), verified size-matched

Matrix: ATAC (6.4 GB)
gs:// path: gs://2025_01_06/shivani_snm3c_methscan/filtered_matrices/ATAC_matrix/
HTTPS base: https://storage.googleapis.com/2025_01_06/shivani_snm3c_methscan/filtered_matrices/ATAC_matrix/
────────────────────────────────────────
Matrix: geneslop2k (14.0 GB, incl. annotated_with_AK.h5ad)
gs:// path: .../filtered_matrices/geneslop2k_matrix/
HTTPS base: https://storage.googleapis.com/2025_01_06/shivani_snm3c_methscan/filtered_matrices/geneslop2k_matrix/
────────────────────────────────────────
Matrix: promoter (2.5 GB)
gs:// path: .../filtered_matrices/promoter_matrix/
HTTPS base: https://storage.googleapis.com/2025_01_06/shivani_snm3c_methscan/filtered_matrices/promoter_matrix/
────────────────────────────────────────
Matrix: VMR (17.2 GB, incl. annotated.h5ad + per-condition *_celltype.txt incl. dex_celltype.txt)
gs:// path: .../filtered_matrices/VMR_matrix/
HTTPS base: https://storage.googleapis.com/2025_01_06/shivani_snm3c_methscan/filtered_matrices/VMR_matrix/

Inputs — processing notebooks (both dirs)
https://storage.googleapis.com/storage/v1/b/2025_01_06/o?prefix=shivani_snm3c_methscan/processing_notebooks/

• gs://2025_01_06/shivani_snm3c_methscan/processing_notebooks/2_process_data/ (8 files)
• gs://2025_01_06/shivani_snm3c_methscan/processing_notebooks/3_methscan/ (52 files)