Edit on GitHub

starplast.chromatin

Chromatin signal per gene, summarised from coverage tracks over promoters.

GEO serves these two experiments as bigWig only -- no peak calls, no per-gene table -- so the number has to be computed here. What is computed is deliberately the dullest possible summary: the mean coverage over a fixed window around each gene's transcription start, normalised to the genome mean. No peak calling, no background model, no thresholds. Every one of those is a modelling choice that would be invented here rather than taken from the authors, and the project's rule is that a number which looks like a measurement must be one.

Why the promoter and not the gene body

Accessibility and factor occupancy are promoter phenomena; over a gene body they mostly report how long the gene is. The window is PROMOTER bases either side of the start, taken from the strand, so a reverse-strand gene's promoter is at its higher coordinate. Getting that backwards is silent -- the column still looks like data -- which is why the loader is tested on a reverse-strand gene.

What only the unperturbed arm is used

Both deposits are knockdowns, and both include an untreated arm. Only the untreated arm is read: it is the measurement of what chromatin looks like in a normal parasite, which is what the slot asks. The perturbed arms measure what a specific depletion does, which is a different question and would need its own slot rather than being averaged into this one.

PROMOTER = 1000
LOCATION_TABLE = 'toxodb_gene_location.tsv'
LOCATION = re.compile('^(?P<chrom>[^:]+):(?P<start>[\\d,]+)\\.\\.(?P<end>[\\d,]+)\\((?P<strand>[+-])\\)')
def gene_windows(base: str, resolve=None) -> pandas.DataFrame:

Promoter windows per gene: sequence, start, end, taken from the strand.

def track_means(path: str, windows: pandas.DataFrame) -> pandas.Series:

Mean coverage in each window, as a fraction of the track's genome-wide mean.

Normalised to the track rather than left raw, because two tracks differ in sequencing depth by whatever the submitters happened to load. A ratio to the track's own mean is comparable between them and is what makes averaging replicates meaningful.

def readable(path: str) -> bool:

Whether this file is a bigWig at all, by its magic number.

GSM8524430_UT_2.bw in GSE277553 begins with eight 0xFF bytes and is not a bigWig. It is the same 6,943,536 bytes whether taken from the series tar or fetched from GEO as a sample file, so the corruption is in the deposit and not in the download -- which is the only reason it is correct to skip it rather than to re-fetch. It is skipped LOUDLY: a replicate silently dropped is a mean over fewer samples than the note beside the column claims.

def chromatin_signals( base: str, resolve=None, log=<built-in function print>) -> pandas.DataFrame:

ATAC accessibility and HDAC3 occupancy over promoters, from the untreated arms.