Quickstart

A 10-minute tour of genomeblocks. We’ll load ATAC peaks, annotate them with a GTF, build a chromatin-contact graph from HiChIP loops, compute O/E weights from a Hi-C matrix, and render a multi-track browser view.

Table of contents

  1. 1. Load peaks as Loci
  2. 2. Interval algebra
  3. 3. Annotate with a GTF
  4. 4. Annotate and group CREs
  5. 5. Build a chromatin-contact graph
  6. 6. Compare two Loci sets as a heatmap
  7. 7. Render a browser view
  8. Next steps

1. Load peaks as Loci

from genomeblocks import Loci

peaks = Loci.make("atac_peaks.narrowPeak")
print(peaks)                     # Loci(n=73412)
print(peaks[0])                  # Locus(chrom='chr1', start=9800, end=10100, strand='.')
print(peaks[0].uid)              # 'chr1:9800-10100(.)'

Loci is a list subclass, so indexing, iteration, and len() work as expected. It also exposes an on-demand cgranges index for O(log n) overlap queries.


2. Interval algebra

# Extend each peak by ±200 bp, then collapse overlapping peaks.
cre = peaks.slop(200).sort().merge()

# Boolean set operations
superenhancers = Loci.make("H3K27ac_SE.bed")
se_cre = cre & superenhancers    # intersection
only_cre = cre - superenhancers  # difference
combined = cre | superenhancers  # union
xor = cre ^ superenhancers       # symmetric difference

See the Loci guide for the full API.


3. Annotate with a GTF

from genomeblocks import Genes

genes = Genes.make("gencode.v38.annotation.gtf", promoter_r=1000)

# Per-CRE genomic context (Promoter-TSS / 5UTR / 3UTR / Exonic / Intronic / Intergenic)
annot_df = genes.annotations(cre)
print(annot_df.groupby("annotation").size())

# Nearest gene via pyranges
near_df = genes.nearest_genes(cre)

Genes also parses UCSC RefSeq tables (Genes.make_ucsc(...)) and supports alt-promoter-aware transcript-level operations. See Genes guide.


4. Annotate and group CREs

Label CREs by gene context, and split them with set algebra. Grouping for plots is just a plain dict[str, Loci] — no separate annotation object:

prom = Loci(list(genes.annot['prom']))

annot  = genes.annotations(cre)                  # region class per CRE uid
counts = annot["annotation"].value_counts()

groups = {"promoters": cre & prom,               # CREs overlapping a promoter
          "enhancers": cre - prom}               # the rest

groups feeds straight into plot_heatmap(..., groups=groups) (see the AR & FOXA1 walkthrough).


5. Build a chromatin-contact graph

from genomeblocks import Architecture

arch = (Architecture.make(cre, "HiChIP_loops.bedpe", r=2500)
                    .add_mcool(cre, "cohesin.mcool", resolution=5000)
                    .normalize(cre, source="w", name="n"))

print(arch)
# Architecture(name='Skeleton', loci=42311, links=185093, vertex_props=[uid], edge_props=[w, n, d])

Annotate the graph and extract hub genes in one line:

arch.annotate(cre, genes)
result = arch.prime_hubs(key="n")
print(sorted(result["prime_genes"]))

Full API: Architecture guide.


6. Compare two Loci sets as a heatmap

from genomeblocks import compare_heatmap

bigwigs = ["ATAC.bw", "H3K4me1.bw", "H3K27ac.bw"]

fig, union, S, groups = compare_heatmap(
    a=cre,
    b=superenhancers,
    bigwigs=bigwigs,
    a_name="CRE-only",
    b_name="SE-only",
    common_name="shared",
    n_bins=200,
    flank=3000,
    normalize=True,
    cmap=["Blues", "Reds", "Greens"],
    vmax=[10, 5, 8],
)
fig.savefig("compare.pdf")

See Signal guide for the full pipeline.


7. Render a browser view

from genomeblocks import browser
from genomeblocks.bedpe import read_bedpe

fig, _ = browser(
    region=("chr6", 122_600_000, 122_800_000),
    tracks={
        "HiChIP loops": read_bedpe("loops.bedpe"),
        "ATAC signal":  "ATAC.bw",
        "ATAC peaks":   peaks,
        "Genes":        genes,
    },
    figsize=(12, None),      # auto-height
    bw_n_bins=120,
)
fig.savefig("browser.svg")

A full real-data example is in examples/ar_foxa1_lncap/ and the browser walkthrough; both render into an IGV-like, fully-vectorial SVG.


Next steps


Copyright © 2024–2026 Umut Berkay Altintas. MIT Licensed.

This site uses Just the Docs, a documentation theme for Jekyll.