Oncofusion

start / graph / concept-q114-a-written-bigwig-carries-only-the-chromosomes-its-data-touched

Explainerstrong confidence 4 of 5: Replicated in peer-reviewed work by more than one group.Q114-suv39h1-at-the-locus

Two coupled traps in one conversion: pybigtools writes only the chromosomes the input actually touched, so a scaffold in your query set raises KeyError on one track and not another - and when the scaffolds were dropped to fix it, 71 of 59,892 anchors (0.12%) turned out to carry enough of the top signal to HALVE a control's numerator

also called bedGraph to bigWig chrom table · unplaced scaffolds carry the top signal · the query universe must be intersected per file

Converting one file format to another quietly dropped the odd leftover pieces of the genome, which made the program crash on some files and not others. Fixing it revealed something more interesting: those leftover pieces, barely a tenth of a percent of the sites being compared, were carrying a large share of the strongest signal.

The same thing, for a scientist

**Trap one: the written chrom table is not the chrom table you passed.** `pybigtools` 0.3.0's `BigWigWriter.write(chroms, vals)` takes a chromosome dictionary, but the resulting file's own `chroms()` contains **only the chromosomes the value iterator actually emitted** [E1]. A GEO bedGraph produced by HOMER covers 25 chromosomes; the reference table used to write it had 195. Reading a query window on `chrUn_GL000219v1` from the written file then raises ``` KeyError: 'No chromomsome with name `chrUn_GL000219v1` found.' ``` **mid-run, after two tracks had already been profiled.** The failure is not silent - it throws - but it is *asymmetric*: the same query set works on the four bigWigs downloaded from GEO (195 and 455 chroms) and fails on the ones written locally. An analysis that caught the exception and continued would have scored different tracks against **different numbers of anchors** and nothing in the output would have said so. **Trap two, which is the finding and not the plumbing.** The fix was to restrict every query set to the primary assembly, dropping **71 of 59,892** SK-N-SH CTCF anchors and **32 of 56,203** MCF-7 ones - **0.12% and 0.06%** [E2]. That also changed a control: | track | top-0.1% anchor mean, all anchors | primary assembly only | |---|---:|---:| | SU-DIPG36 SUV39H1 | **42.8** | **21.5** | | SU-DIPG13 SUV39H1 | 27.6 | 22.2 | | SU-DIPG13 H3K9me3 | 26.0 | 16.5 | | SU-DIPG36 H3K9me3 | 36.5 | 18.2 | **A tenth of a percent of the anchors, all on unlocalised scaffolds, were carrying enough of the extreme signal to halve the statistic in one track and lower it in all four.** For an H3K9me3 writer this is expected rather than shocking - SUV39H1's canonical address is pericentromeric satellite, which Bowtie2 with one copy per position cannot place, so reads pile onto whatever unlocalised scaffold is the best available home. It is still the sort of thing that silently decides whether a control passes. **The general form.** A query universe is part of the measurement. **Intersect it with each file's own chrom table, count what that removed, and report the count** - the same discipline rule 78 already demands for a symbol-keyed join that drops multi-mapping genes. Here the drop was 0.12% and it moved a control by 2x, which is the argument for counting rather than for any particular threshold.

Go deeper
Three smaller mechanics from the same conversion, worth writing down once
  • pybigtools 0.3.0 refuses a Windows drive-letter path for WRITING. pybigtools.open(r'C:\Temp\x.bw', 'w') raises ValueError: Invalid file path. Writing does not support urls. - C: is parsed as a URL scheme. A relative path works, and so does the extended-length form //?/C:/Temp/x.bw. Reading a drive-letter path is unaffected, which is why this only appears once you write.
  • missing= is deprecated in favour of fillna=, and the default for uncovered positions changed from 0 to NaN between versions with a DeprecationWarning rather than an error. All three modes were probed per file here and returned identical values, because both deposits write explicit zero-valued records - but that is a property of these deposits, not of the library.
  • One deposit can carry two different chrom tables. GSE225028's two H3K9me3 bigWigs have 455 chromosome entries and its two SUV39H1 bigWigs have 195, from the same submitter on the same day. The existing registry row on bigWig practicalities already says to verify the build per file, not per series; this is an independent rediscovery of that row in a new deposit, which is worth recording because it shows the trap is live.
Facts and figures
librarypybigtools 0.3.0, Python 3.12.13, Windows Server 2022
trap 1BigWigWriter writes only chromosomes present in the value iterator; the written file's chroms() is a subset of the dict passed in
trap 1 symptomKeyError: 'No chromomsome with name `chrUn_GL000219v1` found.' - asymmetric between downloaded and locally written files
trap 2 pathpybigtools.open('C:/...','w') -> ValueError: Invalid file path. Writing does not support urls. Workarounds: a relative path, or '//?/C:/...'. Reading a drive-letter path is unaffected.
trap 3 fillnamissing= is deprecated in favour of fillna=; the uncovered-base default changed from 0 to NaN with a DeprecationWarning, not an error
interpretationSUV39H1's canonical address is pericentromeric satellite, which a unique-alignment pipeline cannot place; extreme values concentrating on unlocalised scaffolds is the expected signature and it still moved a control by 2x
ruleintersect the query universe with each file's own chrom table, count what that removed, and report the count - the same discipline rule 78 demands for a symbol-keyed join
trap 4 two chrom tables in one deposit
{
  "GSE225028_H3K9me3_bigwigs": 455,
  "GSE225028_SUV39H1_bigwigs": 195,
  "registry_row_this_rediscovers": "bigWig practicalities: padded chrom.sizes, NaN gaps, and a summary= that is silently ignored - which already says verify the build PER FILE, not per series"
}
anchors dropped by the primary assembly restriction
{
  "SK-N-SH ENCFF685KTA": "71 of 59,892 (0.12%)",
  "MCF-7 ENCFF942TCG": "32 of 56,203 (0.06%)"
}
effect on the c1 numerator
{
  "GSM7037532": [
    42.8,
    21.5
  ],
  "GSM7037531": [
    27.6,
    22.2
  ],
  "GSM7037529": [
    26,
    16.5
  ],
  "GSM7037530": [
    36.5,
    18.2
  ]
}
How this could bite us back
  • The interesting half of this entry is a by-product of a crash. An analysis that had caught the KeyError and continued would have produced a clean-looking table with different denominators per track, and nothing downstream would have flagged it - which is the argument for letting instrument errors propagate rather than handling them.
Still open
  1. Which 71 anchors were they? The identities were not persisted, so the claim that they are pericentromeric-proximal is an inference from the scaffold names' class rather than a measurement. Rule 90 says persist the members of any set you report the size of, and this beat did not.
  2. Do this project's other bigWig analyses query scaffolds they should not? Q35, Q37, Q110, Q111 and Q112 all read ENCODE bigWigs against peak-derived anchor sets. ENCODE files carry full chrom tables so none of them would have thrown - which means any scaffold contamination in their nulls is silent rather than fatal.
  3. Does the SUV39H1 signal on unlocalised scaffolds recover a real pericentromeric distribution if the tracks are re-processed with multi-mapping reads retained? That would be a much better positive control for a SUV39H1 antibody than the KRAB-ZNF cluster, which failed here at every scale.
Every source behind this entry (2)
  1. E1established confidence 5 of 5: Multiple independent lines of evidence; treated as settled in the field.database

    pybigtools 0.3.0's BigWigWriter emits only the chromosomes present in the value iterator, so reading a query window on an absent scaffold from the written file raises KeyError; and it refuses a drive-letter path for writing while accepting it for reading.

    What this evidence cannot show

    Observed on Windows Server 2022 with pybigtools 0.3.0 and Python 3.12.13. The path behaviour is plainly platform-specific; the chrom-table behaviour should not be, but was only tested on this platform and this version.

    Evidence that argues against this (1)
    • This may be intended behaviour rather than a defect - writing chromosome entries with no data would bloat the header. The hazard is not that it is wrong but that it makes two files produced by one pipeline queryable over different universes. source
  2. E2strong confidence 4 of 5: Replicated in peer-reviewed work by more than one group.dataset

    Dropping 71 of 59,892 SK-N-SH CTCF anchors (0.12%), all on unlocalised or alt scaffolds, lowered the mean of the top 0.1% of anchor windows from 42.8 to 21.5 in one SUV39H1 track and from 27.6 to 22.2 in the other, and from 26.0 to 16.5 and 36.5 to 18.2 in the two H3K9me3 tracks.

    What this evidence cannot show

    The two runs differ only in the chromosome filter, and the before-values are quoted from the first run's console output rather than from a retained results file; the after-values are in results_bigwig_posthoc.json. The individual scaffold anchors were not listed, which is a persist-the-identities failure of exactly the kind rule 90 describes and should be fixed if this is re-run.

    Evidence that argues against this (1)
    • A 0.1% top-tail statistic over 60 windows is unstable by construction, so a shift of this size from removing 71 candidate windows is arithmetically unsurprising. What is not automatic is that the removed windows were concentrated in the extreme tail at all. source
Connected to