Harden end-to-end HiChIP workflow, QC, and provenance - #1
Merged
Conversation
Runs the pipeline end to end on GSE101498 (Mumbach et al., Nat Genet 2017): primary human CD4+ T-helper subsets, H3K27ac HiChIP, three biological replicates each, plus CTCF HiChIP in naive cells as an architectural-anchor arm. 11 libraries, 2.15 billion read pairs. MACS2 could not run. The bioconda macs2 2.2.x binaries link __log_finite, a libm symbol glibc dropped at 2.27, so on any current distribution it dies inside `import MACS2.IO.ScoreTrack` before parsing an argument. The anchor step is the step every loop depends on, so nothing downstream of it was reachable. Moved to MACS3: same model, same lab, flag-compatible. H3K4me1 was classified narrow. It is broad — the signal is a wide bimodal shoulder around the nucleosome-depleted region, not a sharp peak. For HiChIP this matters more than for ChIP-seq: anchors are what loops hang off, so a wrong anchor shape does not merely blur a peak, it merges distinct enhancer-promoter contacts into one uninterpretable loop. Stripe calling added (rules/06b_stripes.smk, Stripenn). A stripe is a continuous line of contact enrichment anchored at one end — the signature of one anchor reeling in a domain by loop extrusion — and it is not recoverable from a loop list, because a loop caller tests discrete pixel pairs and fragments a stripe into individually unconvincing marginal pixels. The pipeline called loops, insulation and compartments but had no stripe caller. Prompted by Banecki et al. (Commun Biol 2025), who make stripes a primary readout of HiChIP. Technical replicates are merged into their biological replicate at fetch time. They share a library, so carrying them separately would inflate n and understate the between-donor variance the differential test has to beat. Differential comparisons defined (were empty): Th17 vs Naive, Treg vs Naive and Th17 vs Treg on H3K27ac. CTCF is in no comparison — it exists in naive cells only, and contrasting loops called on a CTCF anchor set against loops called on an H3K27ac anchor set would compare anchor definitions, not biology. Viz loci switched from generic cancer genes to the lineage-defining ones (RORC, IL17A, FOXP3, IL2RA, CTLA4, plus MYC as shared control), with coordinates generated from the GENCODE GTF rather than typed by hand. download_resources.sh no longer skips index building in silence when the indexer is absent; a skipped index failed hours later inside the first alignment job. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
…loop counts FitHiChIP is not on bioconda. envs/fithichip.yaml declared `- fithichip`, which resolves from no channel, and the run rule's fallback told the user to `mamba install -c bioconda fithichip`. So the loop caller could not be installed by any path the repo described — and every loop, every differential loop and every APA in this pipeline depends on it. It is distributed as a GitHub release of shell/R/python scripts; the pinned 11.0 release is now fetched by a rule and its dependencies come from the upstream Dockerfile. Pinned to a tag rather than master because both things this workflow touches have changed across versions: the config keys it writes, and the output directory nesting it parses. The generated config did not match 11.0. It carried Draw, TimeProf and HiCProBasedir, which 11.0 does not read, and omitted MergeInt, which it does. Without MergeInt, adjacent significant bins are emitted as separate loops instead of being merged into one contact — inflating the loop count and letting the same interaction be counted several times in the differential test downstream. Stripenn could not import cv2: the opencv-python wheel it pulls needs a newer libstdc++ than the host provides and carries no RPATH back into the environment, so the loader fell through to /lib/x86_64-linux-gnu and died on a missing CXXABI version. The environment was right and the search order was wrong — the same shape as a host ~/.Rprofile or ~/.local numpy shadowing a conda env. Fixed by putting the env's own lib first on LD_LIBRARY_PATH. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
~/.Rprofile on this machine runs .libPaths("~/Rlibs"), which prepends a library
built against a different R ahead of the conda environment's own. FitHiChIP's R
steps then cannot load edgeR, ggplot2, data.table or dplyr, even though all four
are installed in the environment it is running in — requireNamespace returns
FALSE for packages whose directories are plainly there.
The profile is sourced on every R startup, after the environment is read, so
clearing R_LIBS_USER is not enough; R_PROFILE_USER has to be neutralised too.
This would have surfaced only after the alignment of 2.15 billion read pairs, as
an R error inside the loop caller.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
MACS was being run with -f BAMPE on the HiChIP 1D BAM. That asks it to read the two ends of a pair as the two ends of one sonicated fragment — but in HiChIP the two ends are ligation partners, routinely megabases apart or on different chromosomes, because that separation IS the contact being measured. The implied "fragment length" is the loop span. MACS3 will not do it: it builds an empty paired-end track and dies dividing by the fragment count. Anchors are now called from the read ends individually, each extended to a nucleosome (--nomodel --extsize 147), which is what FitHiChIP's own peak inference does. That also means the anchors this stage produces and the anchors FitHiChIP would infer for itself are defined the same way, rather than by two different models. Input is BED rather than BAM because MACS3's BAM reader returns zero tags from this file and then divides by that zero; bedtools bamtobed of the same BAM yields 5,957,598 tags and 11,056 anchors on the shallowest library. FitHiChIP feeds MACS a BED here too. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
The sorted pairsam is consumed by pairtools_dedup and never read again, but it is the largest thing the workflow writes: roughly 0.15 GB per million pairs, so around 500 GB across a 2.15-billion-pair cohort, sitting on disk for the whole run alongside a deduped copy of itself. The trimmed FASTQs are read only by the aligner. Both are now temp(), so Snakemake reclaims them as soon as their consumer finishes. Without this the run outgrows the array before it reaches loop calling — on this cohort it was heading for ~1.2 TB of intermediates against 975 GB free. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
count_per_loop snapped each anchor to its bin start and read a single res-sized bin (e1 = s1 + res). That was right when FitHiChIP emitted one bin per anchor, but MergeInt=1 merges adjacent significant bins, so anchors are commonly 2-4 bins wide. Each loop was then counted on a fraction of its own footprint, and the fraction shrinks as the anchor widens -- i.e. it tracks ChIP enrichment, the very thing that differs between the groups being compared, so it enters DESeq2 as a group-correlated bias rather than as noise. On 300 synthetic multi-bin anchors against a real 5 kb matrix, single-bin counting recovered 22% of the pairs the full anchor rectangle contains. Counts are now summed over the full rectangle of bins the two anchors span, and taken from one streaming pass over the pixel table instead of a cooler .fetch() per loop: the union set is O(1e5) loops and each fetch is an indexed HDF5 range read, which costs hours per sample. Verified pixel-identical to per-loop rectangle fetch on Th17_H3K27ac_rep1 (300/300 loops). Also drop the blanket `except Exception: counts.append(0)`, which turned a systematic failure (mismatched chromosome naming, wrong bin grid) into a silent all-zero count matrix, and stop FitHiChIP's text header row from being emitted into the union BEDPE as a loop. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
…cordance Three separate faults in the only replicate-concordance metric in the pipeline. hicrepSCC returns a numpy array indexed by chromosome, not a dict, so the `scc.items()` call raised AttributeError -- the step could not complete at all. The array is pre-filled with the sentinel -2.0 for chromosomes it cannot score (`np.full(len(chrNames), -2.0)`), which np.nanmean does not skip, so unscored chromosomes would have dragged SCC toward -2 had it run. Chromosomes are now excluded through hicrepSCC's own excludeChr argument, and the sentinel is masked. bDownSample is passed explicitly instead of receiving the bin size positionally. Replicates were grouped on subject_id + mark, which in a cohort of three T-cell subsets taken from the same three donors puts Naive, Th17 and Treg from one donor into a single "replicate" group: it reported the correlation BETWEEN CELL TYPES as replicate concordance -- scoring the exact biology the differential test is looking for against a threshold that assumes it is noise. Grouping is now cell_type + mark. The absolute 0.85 pass threshold is a Hi-C convention and does not transfer. HiCRep downsamples the deeper library to the shallower one, so a pair is only as good as its worst member. Measured here at 25 kb: two libraries from DIFFERENT cell types (90.8M and 47.6M contacts) score 0.743, while two genuine replicates score 0.221 because one holds 3.0M. The cutoff would have failed a real replicate pair and passed a cross-cell-type pair. SCC is now a sanity floor, and any pair whose shallower member falls under min_contacts_for_scc is recorded as depth_confounded and cannot decide PASS/FAIL. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
APA divided the centre pixel by the mean of all four 3x3 corners. Contact frequency inside the window is a function of distance from the diagonal: a pixel (i, j) sits at genomic separation D + (j - i) * bin_size, so the four corners span D +/- 2 * win * bin_size -- at window_size=20 and 10 kb bins, a 400 kb swing. The bottom-left corner is therefore ~400 kb closer to the diagonal than the loop itself and carries far more contacts for reasons that have nothing to do with the loop. Averaging it with the others inflated the denominator and deflated every APA score, and that deflated score was the one gating pass/fail. Background is now the two corners on the j - i = 0 anti-diagonal, which sit at the same genomic separation as the centre pixel, and pass/fail is decided on the random-shift control -- which moves both anchors by one offset and so is distance-matched by construction. It was already being computed and then ignored. The random shift could also be arbitrarily small, sliding the "control" window back onto the real loop so the denominator absorbed the signal it exists to measure; it is now required to clear the window. apa.score_min and hicrep.threshold_pass are read from config rather than hardcoded. loop_qc_summary carried its own copy of the HiCRep threshold (0.85), which disagreed with config and silently won in the report. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Every plot the pipeline produced was per-sample: one APA per library, one browser
panel per locus, one volcano per comparison. Nothing was assembled across the
cohort, so the questions a reader actually asks of a HiChIP experiment had no
figure that answered them -- is the loop count tracking biology or sequencing
depth, do the T-cell subsets separate, are the stripes a property of the anchor
mark rather than of the cells.
Five figures, vector PDF plus 400 dpi PNG:
1 library QC depth against the usability floor, duplication, cis fraction,
and P(s) decay; libraries under the floor are hatched and
dashed rather than dropped, because a shallow library that
silently disappears is how a depth artefact becomes a result
2 reproducibility compartment E1 correlation across all libraries (robust to
depth) beside HiCRep SCC with depth-confounded pairs flagged
3 loops loop counts, loop yield AGAINST depth, span distribution, and
APA aggregated per cell type
4 differential volcano per comparison plus the span distribution of gained
and lost loops
5 stripes split by anchor mark, never averaged across them: on CTCF a
stripe is an extrusion anchor, on H3K27ac it is not, and
pooling the two would manufacture a difference
Every panel that can be empty renders "no data" instead of raising, so one
loopless or failed library cannot destroy the figure set at the end of a
multi-day run. apa_plot now persists its aggregate matrix (.npy) so figure 3 can
render real APA panels rather than re-deriving them.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
virtual_4c read one row of the contact matrix by calling clr.matrix(...).fetch(chrom), which returns a DENSE n x n array for the whole chromosome. chr1 at 5 kb is 49,792 bins, so reading a single viewpoint row asked the allocator for 19.8 GB -- per job, with viz jobs running concurrently. It now fetches the viewpoint row against the chromosome, which is the same numbers out of the same matrix: measured peak RSS falls from 19.8 GB of allocation to 0.19 GB, and the profile is identical (49,792 bins). The rule was also reachable from no target, so it had never run. Virtual 4C is the view that shows an enhancer looking at its target gene -- the RORC / IL17A / FOXP3 question this cohort exists to answer -- so it is now built by `rule all`. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Twelve tests over the parts that were wrong and are easy to get wrong again: anchor width (a 15 kb anchor covers three 5 kb bins, not one), chromosome bin offsets, the streaming loop counter against a direct rectangle fetch on a cooler whose every pixel is known, FitHiChIP's non-comment header row, which APA corners are comparable to the centre pixel, and HiCRep's -2.0 sentinel surviving np.nanmean. Two are regression tests that fail against the old code. The scripts now guard `main(snakemake)` on `snakemake` being in globals, which Snakemake injects and a test import does not, so they can be imported without executing. count_per_loop's counting is lifted into count_loops(clr, loops, res). loop_qc_summary was gating on apa_score (the corner ratio) while apa_plot decided its own `pass` on apa_vs_random_shift (the distance-matched control), so the two could reach opposite verdicts about the same library. Both now read the control. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
TUTORIAL told users the `mark` column must match a key in `config.yaml:macs2.marks` -- a key that no longer exists -- and BEST_PRACTICES listed H3K4me1 as a narrow mark, which is the exact error that was fixed in the config: it is broad per ENCODE. A doc that names a dead config key or recommends the bug is worse than no doc. Also: MACS2 -> MACS3 throughout (MACS2 cannot load on current glibc), the stale "HiCRep SCC >= 0.85" acceptance criterion replaced with what the metric actually supports on a HiChIP matrix, and stripes / cohort figures added to the stage list and the results tree. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
It appeared in four places across README, TUTORIAL and BEST_PRACTICES as a non-negotiable, and this cohort disproves it: HiCRep downsamples the deeper library to the shallower, so SCC grades the worst member's depth. At 25 kb, two libraries from different cell types (90.8M and 47.6M contacts) score 0.743 while two genuine replicates score 0.221 because one holds 3.0M. Following the old rule fails the real replicate pair and passes the cross-cell-type pair. Replaced with what the metric supports: group on condition + protein rather than donor, treat SCC as a sanity floor, and disregard it entirely below the contact floor. Two new checklist items make depth the first confounder to rule out, since loop count, stripe count and SCC all rise with unique cis pairs. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Smoke-tested the figure stage against the partial results rather than waiting for the end of the run to find out. It rendered all five figures without raising -- the graceful-empty handling works on maximally incomplete input (8 of 11 libraries with no stats, no loops, no APA, no HiCRep) -- but figure 1 was not publishable: the 90-degree sample labels under panels a and c ran into the title of the panel below, the usability-floor annotation sat on top of the tick labels, and panels b/c spent half their axes on samples that had no data. Sample ids are shortened for display (Naive_H3K27ac_rep1 -> "Naive K27ac r1"), panels drop rows with no measurement, and the depth-floor label moved outside the axes. Colour was also the only thing separating cell types and nothing said so, so the key is now on the figure -- where it survives an empty panel -- rather than inside panel d. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
…HiC-Pro
Ran FitHiChIP against a real finished library rather than waiting for the cohort
to finish, and it quit before calling a single loop:
ERROR ===>>>> HiC-pro is not installed in the system - FitHiChIP quits !!!
FitHiChIP 11.0 accepts validPairs only if it can build the contact matrix and bin
interval files from them, and it does that by shelling out to HiC-Pro
(FitHiChIP_HiCPro.sh:399-406). HiC-Pro is not installed, is not on bioconda, and is
precisely the legacy stack this pipeline replaced with bwa-mem2 + pairtools +
cooler. So loop calling would have failed for all 11 samples -- and it would have
failed only after every alignment completed, since nothing upstream depends on it.
cooler already builds the matrix, so it is handed over directly. `COOL=` takes the
5 kb single-resolution {sample}.base.cool, which FitHiChIP reads with
`cooler dump -t pixels --join`. It must be the plain .cool file and not an
.mcool::/resolutions/5000 URI: FitHiChIP validates the path with `[ ! -f ... ]`,
which a URI fails. Verified past the gate on Th17_H3K27ac_rep1 -- it now builds the
matrix chromosome by chromosome.
This also deletes rule pairs_to_validpairs, which existed solely to feed that dead
path and wrote a 2.4 GB intermediate per sample (~26 GB across the cohort).
Separately confirmed the R isolation added earlier is load-bearing, not defensive:
with the host ~/.Rprofile active, .libPaths() is prepended with ~/Rlibs and edgeR,
ggplot2, data.table and dplyr all resolve FALSE in the FitHiChIP env; with
R_PROFILE_USER=/dev/null they all resolve TRUE.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Watching the real FitHiChIP run showed it builds TWO interaction-type directories:
the requested FitHiChIP_Peak2ALL_..., and an ALL2ALL one produced while computing
the background model. The rule collected its result with
find "$(dirname {output.loops})" -name "*interactions_FitHiC_Q0.01.bed" | head -n 1
which searches the whole sample tree and takes whatever comes first. `find` returns
directory order, not sorted order, so the rule could pick up the ALL2ALL call set
instead of the Peak-to-ALL one that was configured -- silently, with no error, and
not necessarily making the same choice for every sample.
The output directory name is now derived from the same numeric codes that get
written into the .config (IntType 3 -> Peak2ALL, BiasType 1 -> Coverage_Bias,
UseP2PBackgrnd 1 -> P2PBckgr_1), so the path read from cannot drift from the config
written. Confirmed identical to the directory the live run produced. Any fallback
glob is confined to that directory, and failure now prints which BEDs did exist.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Ran FitHiChIP to completion on Th17_H3K27ac_rep1 (exit 0, 30,903 loops at q<0.01, median span 105 kb, 98.7% of loops with both anchors on a peak, no trans) and then looked at what it actually wrote. Two corrections follow. 1. MergeInt=1 was having no effect on our output. FitHiChIP writes TWO results: the raw significant bin pairs, and -- under Merge_Nearby_Interactions/ -- the same set with adjacent significant bin pairs collapsed into one contact. The rule was reading the raw file. Setting MergeInt=1 only makes FitHiChIP WRITE the merged one; it still has to be read. So the double-counting that commit 9c2b794 claimed to fix was still entirely present: 30,903 raw calls against 17,875 merged, i.e. 42% of the raw set is the same physical loop reported once per adjacent bin pair it spans, and every one of those duplicates would have entered DESeq2 as an independent observation. 2. Anchors are NOT multi-bin. Commit 16118c1 justified counting over the whole anchor rectangle by asserting that MergeInt=1 makes anchors "routinely 2-4 bins wide", and quoted 22% of signal recovered. That is wrong: FitHiChIP's merge picks a representative bin pair rather than widening the anchor, and every anchor in the real output -- raw and merged alike -- is exactly one 5 kb bin. The 22% came from synthetic multi-bin anchors that FitHiChIP never emits. The rectangle logic is kept, because the union BEDPE is not required to come from FitHiChIP and a wider-anchor caller would otherwise be counted on an arbitrary fraction of its footprint. But it is defensive here, not load-bearing, and the docstring now says so. The load-bearing half of that commit -- one streaming pass instead of a cooler .fetch() per loop -- stands. The result path is derived from config (fithichip.merge_nearby), verified identical to the file the live run produced. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
results/pairs was 86 GB for five libraries and heading for ~375 GB across eleven,
against 764 GB free with 383 GB of transient fasterq-dump scratch also on the disk.
Two of the three files pairtools_dedup writes per library are dead weight:
unmapped.pairs.gz 12.6 GB/library, consumed by NO rule at all. It is a pure
byproduct of --output-unmapped. What QC actually wants is
the unmapped COUNT, and that is already in the stats file.
~139 GB across the cohort.
dedup.pairsam.gz 18.6 GB/library, and its only consumer is pairs_to_1d_bam.
Once the 1D BAM exists nothing reads it again. ~205 GB.
Both are now temp(). dedup.pairs.gz (the cooler input) and the stats file are
untouched.
Marked temp() BEFORE deleting the existing copies, which is the part that matters:
delete a non-temp output and Snakemake sees a missing output and re-runs the rule
that made it, cascading back through alignment. Verified with a dry run either side
of the deletion -- pairtools_dedup stays at 6 jobs (the libraries that still need it)
rather than jumping to 11, and the five finished libraries keep their mcool, peaks
and 1D BAM. Reclaimed 82 GB immediately; free space 761 GB -> 843 GB.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Snakemake prepends its own preamble (the `snakemake` object, logging setup, imports)
to the top of every file used by a `script:` directive. That pushes any
`from __future__` import below other statements, which Python rejects outright:
SyntaxError: from __future__ imports must occur at the beginning of the file
So EVERY script-based rule failed on contact with real data -- virtual_4c,
cooltools_insulation, cooltools_eigs_cis and hicrep_replicate_qc all died in the
first minutes of the full run, and count_per_loop, apa_plot, loop_qc_summary,
differential_loops, export_oracle_cos, bedpe_annotate, stripe_summary and figures
were all queued to die the same way. The unit tests did not catch it because they
import the modules directly, which is exactly the path Snakemake does not take.
The import was only there to allow `str | Path` and `dict[str, float]` annotations
on older Pythons. Every environment in this workflow runs 3.10.20 or newer (3.10,
3.11, 3.12, 3.14), where both are native syntax, so the import bought nothing.
Stripped from all 14 scripts; all 14 still parse and the tests still pass.
Also fix fastqc_raw, which failed on every sample: it renames FastQC's output to a
canonical name, but these FASTQs are already named {sample}_R1.fastq.gz, so source
and destination are the same path and `mv x x` exits 1 with "are the same file".
The rename now only happens when the names actually differ.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Two regression tests that would each have caught a failure that only appeared when
the full DAG met real data:
test_no_future_imports_in_workflow_scripts
Snakemake prepends its preamble to a `script:` file, so a __future__ import
ends up below other statements and Python refuses to compile it. Every script
rule in the workflow died on this.
test_workflow_scripts_guard_their_main_call
A script that calls main(snakemake) at import time cannot be imported by a
test at all, so none of its logic is reachable from the suite. That is why the
existing tests were blind to the above -- they import the modules directly,
which is not how Snakemake runs them.
The second test immediately found seven more scripts with an unguarded call
(bedpe_annotate, compartments_to_bigwig, cooltools_eigs_cis, differential_loops,
export_oracle_cos, pygenometracks_loops, virtual_4c); all are now guarded, which
costs nothing at run time -- Snakemake still injects `snakemake` into globals -- and
makes them testable.
14 tests pass.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
FitHiChIP wrote no loops for any sample and left an empty log behind every
failure. Two separate faults:
- The generated config resolved COOL to an absolute path but left ChrSizeFile,
PeakFile and OutDir relative. FitHiChIP chdir's into its own release
directory before invoking its R and python steps, so those paths resolved
against the wrong root: it read no peaks and no chromosome sizes and quit
before writing anything. Every path is now resolved.
- The rule redirected `2> {log}`, but FitHiChIP reports its fatal errors on
stdout. The diagnostic was being thrown away, which is why the failure
presented only as "produced no BED". Now `> {log} 2>&1`.
cooltools and HiCRep died on hg38's ~160 unplaced scaffolds, which they do not
skip:
cooltools insulation IndexError: index 0 is out of bounds for axis 0 with size 0
hicrepSCC AssertionError: Contact matrix 1 of chromosome GL000208.1 is empty
A new main_chrom_view rule emits the assembled-chromosome view that
expected-cis, insulation and eigs-cis now all take via --view, and from which
HiCRep takes its chrNames. This is not only a workaround: insulation,
compartments and stratum-adjusted correlation are all defined on a chromosome
with a real distance-decay profile, and mean nothing on a 60 kb contig even
where they happen to compute. HiCRep additionally drops chrX, since donor sex
is not recorded in this cohort and X copy number therefore varies between the
libraries being correlated.
`cooltools eigs-cis` also took a -p/--nproc it does not accept (expected-cis and
insulation do); it aborted with "Error: No such option '-p'".
Verified: Treg_H3K27ac_rep1 now runs FitHiChIP to exit 0 and yields 354 loops at
q<=0.01. 14/14 tests pass.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
The stripe summary reported 0 stripes for all 11 libraries and read as a clean
result. It was four faults stacked, and the first one hid the rest:
- `|| true` on the stripenn call, plus a fallback that materialised an empty
table when no output appeared, made a CRASH indistinguishable from a genuine
"no stripes found" -- for every library, silently. Removed; the rule now
fails when stripenn fails.
- stripenn asks on stdin whether it may clear a pre-existing output directory,
and Snakemake always pre-creates the parent directory of a rule's output, so
that prompt always fires. With no terminal it reads EOF and click aborts.
Every zero came from this. Answered via a redirect -- NOT `yes Y |`, which
trades the bug for another: yes takes SIGPIPE when stripenn stops reading and
under `set -o pipefail` that failed the rule after a run that had succeeded,
whose output Snakemake then deleted as incomplete.
- --minL/--maxW are counts of BINS, not base pairs (minL gates a detected line
length in bins; maxW is compared against bin indices). Passing bp straight
through demanded a 200 Mb minimum stripe length at 10 kb -- longer than any
human chromosome, so nothing could ever be called even had it run. Config now
states them in bp and the rule converts.
- hg38's ~160 unplaced scaffolds killed it outright: stripenn takes a quantile
of each chromosome's non-zero pixels, which is empty on a scaffold with no
contacts (IndexError). Its own filter only drops JH5/GL4/RANDOM, which does
not match GL000008.2 / KI270*. Restricted to the same assembled-chromosome
view cooltools and HiCRep already take.
And it must run on RAW counts. Handing it the ICE-balanced matrix returns zero
stripes genome-wide: stripenn is an image method, cooler's balancing leaves NaN
in every unmappable bin (24.2% of the chr1 matrix; 3,224 rows entirely NaN), and
NaN survives the clip so Canny finds no edges. Measured on Naive_CTCF_rep1/chr1,
identical thresholds: balanced 0 candidates / 0 stripes, raw 55 / 38. stripenn's
own matrix handle is `unbalLib`; it controls coverage bias with percentile
saturation and a background null model, not with balancing. --norm is now set to
None explicitly, so a cooler that later carries a KR column cannot silently
change the answer.
HiCRep, separately: every sample in a replicate group is handed the same input
list, computes the same pairs, and averaged all of them -- so the group got one
shared verdict, dominated by its good pairs. Naive_H3K27ac_rep1 holds 3.3M reads
and correlates with its replicates at 0.23, but its group also contains
rep2-vs-rep3 at 0.867, and rep1 was written out as PASS. A per-sample QC file now
scores only the pairs that sample is in, and grades on the best-agreeing
replicate rather than the mean, so one bad sibling cannot fail a good library.
loop_qc_summary reads the same statistic, so the two files cannot disagree.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
…hing
All three differential contrasts failed. pyDESeq2 0.4 does
self.varm["non_zero"] = ~(self.X == 0).all(axis=0)
self.non_zero_idx = np.arange(self.n_vars)[self.varm["non_zero"]]
assuming varm returns the 1-D boolean it was handed. From anndata 0.11 varm is
strictly 2-D and coerces that assignment to (n_vars, 1), so the second line
indexes a 2-D mask:
IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed
The env pinned pyDESeq2 but left anndata to pip, which resolved 0.12.19 against
pyDESeq2 0.4.12. Nothing in pyDESeq2's metadata rules that out. anndata is now
pinned to 0.10.9, the series 0.4 was built against. The count matrices themselves
were fine (6 samples x 135,474 loops, int64, unique keys, no all-zero rows) --
this was purely a dependency clash.
The ORACLE COS export died on Naive_H3K27ac_rep1 with
ValueError: all the input arrays must have same number of dimensions
_loops_to_edges guarded the empty-BEDPE case but not the case where a NON-empty
loop set contributes no edges -- every loop trans, or both anchors landing in the
same bin (i == j), which is routine at coarse resolutions where a sub-100 kb loop
is shorter than one bin. `np.asarray([])` is 1-D, and concatenating it with the
2-D adjacency attributes raises. Naive_H3K27ac_rep1 calls 17 loops, all of which
collapse into single bins, so it was the library that found it. Both empty paths
now return (2, 0) and (0, 3).
Regression test pins the shape contract.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
Both bugs produced figures that looked publishable and were wrong. APA was aggregating the main diagonal, not the loops. A pixel at offset (dy, dx) from a loop of span D sits at separation D + (dx - dy)*bin, and dx - dy spans +/- 2*window -- so the diagonal falls inside the window for ANY loop with D <= 2*window*bin. At the shipped +/-20 bins / 10 kb that is every loop under 400 kb, and the median loop in this cohort is 105-175 kb. Every APA panel, per sample and cohort-wide, was a picture of the diagonal with no corner peak, and it still scored 2-20x "enrichment" because the centre was being compared against corners that were equally contaminated. A passing score and a meaningless image. Now a 21x21 window (+/-10 bins, the standard APA geometry) with a 210 kb floor. The floor is derived and enforced in code as (2*window + 1) * bin_size rather than trusted from YAML, so editing a config value cannot silently reintroduce the diagonal, and a test pins the geometry. How many loops the floor discards is logged and written into the JSON: an APA of 3.5 over 14,000 loops and one over 60 are not the same claim. That fix then exposed a worse one. Naive_H3K27ac_rep1 keeps exactly ONE loop after the distance filter, and one loop's centre pixel over a random-shift control that landed on empty matrix scored 4,874,766 -- reported as a PASS. A ratio built from a single observation cannot fail. APA is now three-state like HiCRep: under min_loops_for_apa (20) the sample is NOT_ASSESSED, recorded and never passed. loop_qc_summary did `(apa_score or 0) >= threshold`, which turned a null score into a silent 0 and reported it as a FAIL -- "these loops are not real" when the truth is "there are too few loops to say". It now carries NOT_ASSESSED through, and a null APA no longer counts toward the overall pass. Figure 1d plotted `dist` (a count of BINS) on an axis labelled "genomic separation (bp)". cooltools emits dist_bp alongside it. The shape of P(s) is unchanged by the rescaling, which is precisely why it was invisible -- the curve looked correct while claiming contacts decay out to 1e4 bp when the data spans ~250 Mb. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
…against depth Figure 2 claimed that biological replicates do not agree on A/B compartments. Naive_H3K27ac rep2 vs rep3 -- a pair HiCRep scores at 0.867 -- came out at 0.20, and several genuine replicate pairs were NEGATIVE. That is not a result, it is a sign error: the sign of an eigenvector is arbitrary, and cooltools solves E1 one chromosome at a time, so without a phasing track "A" is positive on some chromosomes and negative on others, independently in every sample. Correlating raw E1 across libraries averages those random flips toward zero. The compartments were right the whole time; only their orientation was unfixed. eigs-cis now takes a GC content phasing track (the A compartment is the GC-rich, gene-rich one), built at 100 kb by a new gc_phasing_track rule over the same assembled-chromosome view the other cooltools rules use. Measured on the pair above: r = 0.20 -> 0.979, per-chromosome median r = 0.98, and 0 of 23 chromosomes negative, where before many were. Figure 3d grouped APA by cell type alone, which averaged the two Naive CTCF libraries into the same panel as the Naive H3K27ac ones. A TF anchor set and a histone anchor set have different APA profiles by construction -- CTCF scores 20-39x here against 2.8-4x for H3K27ac -- and the mean of the two describes neither. Split by (cell type, mark); the CTCF panel now shows the tight focal peak the H3K27ac panels do not, which is the point. Figure 5b plotted stripes per anchor type as a dot plot, which reads as "CTCF anchors yield fewer stripes". That comparison cannot be made in this cohort: the two CTCF libraries are also the two lowest-complexity ones (55-58% duplicates), and stripe count tracks depth exactly as loop count does. Replaced with stripe yield against depth, where CTCF sits where its depth puts it rather than below where its anchor type would, and the confound is stated on the panel. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
The Pearson r (E1) colourbar label is drawn to the right of the bar and ran into panel b's replicate-pair labels. Co-Authored-By: Claude Opus 4.8 <[email protected]> Claude-Session: https://claude.ai/code/session_01RvYGT9EwbMeuedgJ56FQPE
challagandla
changed the title
challagandla
changed the title
challagandla
marked this pull request as ready for review
July 16, 2026 11:36
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What changed
This PR carries the real-cohort HiChIP workflow through a broader safety and reproducibility audit. It adds fail-fast sample/config validation, restriction-fragment and balancing QC, consensus anchors and loop universes, paired differential contracts, balance-aware Mustache handling, transparent MultiQC/publication reporting, provenance manifests, and hardened ORACLE COS exports.
The audit also corrected silent-failure paths discovered during the earlier GSE101498 run: FitHiChIP path handling, APA geometry and assessment states, per-sample HiCRep verdicts, compartment phasing, stripe runtime behavior, exact differential hypotheses, missing-value masks, and misleading QC labels.
Why
The previous workflow could produce plausible-looking output after several scientifically important failures. The revised contracts make missing references, duplicate SRA runs, invalid replicate designs, non-converged balancing, underpowered comparisons, malformed reporting inputs, and incomplete provenance explicit instead of silently accepting or overstating them.
Validation
git diff --check, secret scanning, and public-files-only checks passed.Interpretation boundary
The earlier committed workflow was run on 11 GSE101498 libraries and produced the historical figures and outputs described in this PR. The newest rerun against the hardened contracts stopped at 21 of 488 scheduled jobs, so those historical biological outputs are stale relative to the final code. The bundled paired contrasts use two complete donor pairs and are deliberately labelled
PILOT_UNDERPOWERED; they are not publication-ready evidence.No raw data, references, generated results, credentials, or internal audit artifacts are included.