BEDPE

Paired intervals as one table: read loop files into Pairs, filter them, intersect their anchors with Loci (bedtools pairtobed), and count raw Hi-C / HiChIP contacts from .pairs, allValidPairs and Juicer files into windows.

Table of contents

  1. What it is
  2. Reading and building
    1. Pairs.make
    2. Pairs.from_frame and as_pairs
  3. Rows, columns and anchors
  4. Filtering
  5. Intersecting anchors with loci: overlapping / anchors_overlap
  6. Export
  7. Pairs files: streaming raw contacts
    1. Contacts per window: count_pairs
    2. Window-to-window matrix: count_pairs_2d
  8. Hand-off to Architecture
  9. End-to-end: enhancer-promoter loops
  10. See also

What it is

gb.Pairs holds a BEDPE as two row-aligned Loci plus extra columns. Pair i is anchor P.a[i] and anchor P.b[i]; both tables share one Genome.

Attribute Content
P.a a Loci: chrom1, start1, end1, strand1
P.b a Loci: chrom2, start2, end2, strand2
P.cols name, score when the file has them, then any further columns (col11, …)
LOOPS.BEDPE chr8 127.73M chr8 127.80M l1 5 chr8 127.73M chr8 128.21M l2 3 chr8 127.74M chr2 41.05M l3 1 … any frame with these columns: Pairs.from_frame · as_pairs Pairs.make tables backend Pairs two Loci aligned by row + columns P.a Loci codes starts ends ± P.b Loci codes starts ends ± 0 127,73… 127,73… + 0 127,80… 127,80… − 0 127,73… 127,73… + 0 128,21… 128,21… + 0 127,74… 127,74… − 1 41,05… 41,05… + P.cols: name · score · … is_cis = a.codes == b.codes · distance = |mid2 − mid1|, ∞ across chromosomes Architecture.make anchors (mid ± r) → overlap_pairs → edges browser · View · igv_html arcs between anchor midpoints Pairs.overlapping(loci) pairtobed: an anchor touches the set every consumer calls as_pairs()
One BEDPE parser, two aligned Loci. Pairs.make reads the file through the tables backend into anchor a and anchor b — two Loci whose row i is pair i — plus the name / score columns, and is_cis / distance are computed from the two coordinate arrays. Everything that takes loops (Architecture.make, the browser, overlapping) calls as_pairs on its argument, so a path, a frame or a Pairs all work there. 
import genomeblocks as gb

P = gb.Pairs.make("loops.bedpe")
P
# -> Pairs(n=4, cis=3, trans=1, cols=[name, score])
P.shape, P.columns
# -> ((4, 10), ['chrom1', 'start1', 'end1', 'chrom2', 'start2', 'end2', 'name', 'score', 'strand1', 'strand2'])
P.to_pandas()
# ->   chrom1  start1  end1 chrom2  start2   end2 name  score strand1 strand2
# -> 0   chr1     900  1100   chr1    4900   5100   l1    5.0       +       -
# -> 1   chr1    1900  2100   chr1   10900  11100   l2    3.0       +       +
# -> 2   chr1     900  1100   chr2     500    600   l3    1.0       -       +
# -> 3   chr2     500   600   chr2    5000   5100   l4    2.0       +       +

Like every table it has len, shape, columns, head(), tail(), describe(), P['score'], an HTML repr, and the Arrow / dataframe protocols (pl.DataFrame(P), duckdb.sql("select * from P")). This module is the one BEDPE reader in genomeblocks: Architecture.make and the browser call it for you.

Reading and building

Pairs.make 

P = gb.Pairs.make("loops.bedpe", min_score=3.0, max_distance=2e6)
P
# -> Pairs(n=2, cis=2, trans=0, cols=[name, score])
arg default meaning
min_score None drop pairs with score below this (needs a score column)
max_distance None drop cis pairs whose anchor midpoints are farther apart — and every trans pair, whose distance is inf
genome None a Genome to share with other tables
backend None the table parser: 'polars' or 'pandas'

The six anchor columns are required (fewer raises); name, score, strand1, strand2 are read when present; columns beyond 10 are kept as col11, col12, … .gz files and leading # header lines are fine. genomeblocks.bedpe.read_bedpe(path, **kw) is the same call.

Pairs.from_frame and as_pairs

Any pandas / polars / pyarrow frame with chrom1 start1 end1 chrom2 start2 end2 columns (by name, case-insensitive, else the first six); strand1 / strand2 are read by name and every other column is kept.

import pandas as pd
df = pd.DataFrame({"chrom1": ["chr1"], "start1": [100], "end1": [200],
                   "chrom2": ["chr1"], "start2": [900], "end2": [1000], "fdr": [0.01]})
F = gb.Pairs.from_frame(df)
F, F.columns
# -> (Pairs(n=1, cis=1, trans=0, cols=[fdr]), ['chrom1', 'start1', 'end1', 'chrom2', 'start2', 'end2', 'strand1', 'strand2', 'fdr'])

genomeblocks.bedpe.as_pairs(x) takes a Pairs (returned as is), a BEDPE or .parquet path, or a frame — it is what Architecture.make and the browser call on their bedpe argument, so all of those forms work there too.

A .bedpe handed to Loci.make raises “use Pairs.make()”.

Rows, columns and anchors 

P = gb.Pairs.make("loops.bedpe")
P[0]                                   # the two anchors of pair 0, as Locus views
# -> (Locus[0](chr1:900-1100(+)), Locus[0](chr1:4900-5100(-)))
P["score"], P["name"]                  # a column
# -> (array([5., 3., 1., 2.]), array(['l1', 'l2', 'l3', 'l4'], dtype=object))
P[P["score"] >= 3]                     # a mask, row numbers or a slice: a new Pairs
# -> Pairs(n=2, cis=2, trans=0, cols=[name, score])
for a, b in P.head(2):                 # iteration yields (anchor a, anchor b)
    print(a.uid, "->", b.uid)
# -> chr1:900-1100(+) -> chr1:4900-5100(-)
# -> chr1:1900-2100(+) -> chr1:10900-11100(+)

Per-pair arrays:

Property Meaning
P.is_cis both anchors on one chromosome
P.distance \|mid2 - mid1\| for cis pairs, inf for trans
P.mids1, P.mids2 anchor midpoints (P.a.centers, P.b.centers) 
P.is_cis, P.distance
# -> (array([ True,  True, False,  True]), array([4000., 9000.,   inf, 4500.]))
P.describe()
# ->                      value
# -> pairs                     4
# -> cis                       3
# -> trans                     1
# -> chromosomes               2
# -> cis distance min       4000
# -> cis distance median  4500.0
# -> cis distance max       9000
# -> score min               1.0
# -> score median            2.5
# -> score max               5.0

The anchors are full Loci, so everything on the Loci page applies to them. The set of distinct anchors is one line:

anchors = (P.a + P.b).merge()
anchors
# -> Loci(n=6, chroms=2, sorted)

Filtering 

P.filter(min_score=3)                  # needs a score column
# -> Pairs(n=2, cis=2, trans=0, cols=[name, score])
P.filter(max_distance=5000)            # trans pairs (distance inf) go too
# -> Pairs(n=2, cis=2, trans=0, cols=[name, score])
P[P.is_cis & (P.distance < 5000)]      # any mask works
# -> Pairs(n=2, cis=2, trans=0, cols=[name, score])

filter returns the same object when nothing is dropped, a new Pairs otherwise.

Intersecting anchors with loci: overlapping / anchors_overlap

P.overlapping(loci) keeps the pairs with at least one anchor overlapping loci — bedtools pairtobed; both=True requires both anchors. r widens the anchors by r bp on each side before the test. loci is anything as_loci takes.

cre  = gb.Loci.make("peaks.bed")
prom = gb.as_loci("chr1:800-1,200")
P.overlapping(prom)                    # either anchor
# -> Pairs(n=2, cis=1, trans=1, cols=[name, score])
P.overlapping(prom, both=True)         # both anchors
# -> Pairs(n=0, cis=0, trans=0, cols=[name, score])
P.overlapping(cre, both=True)          # CRE-CRE loops only
# -> Pairs(n=4, cis=3, trans=1, cols=[name, score])
P.overlapping(gb.as_loci("chr1:1,200-1,300"), r=200)     # anchors widened by ±200 bp
# -> Pairs(n=2, cis=1, trans=1, cols=[name, score])

anchors_overlap(loci, r=0) is the primitive underneath: two boolean masks, anchor a overlaps and anchor b overlaps, aligned to the pairs. Combine them when the two anchors must satisfy different conditions:

genes = gb.Genes.make("gencode.v44.annotation.gtf")
promoters = genes.annot["prom"]
enhancers = cre - promoters

pa, pb = P.anchors_overlap(promoters)          # promoter at a / at b
ea, eb = P.anchors_overlap(enhancers)          # enhancer at a / at b
ep = P[(pa & eb) | (pb & ea)]                  # enhancer-promoter loops
ep
# -> Pairs(n=2, cis=1, trans=1, cols=[name, score])

Both take backend= for the interval engine; every engine gives the same masks.

Export 

P.to_bedpe("filtered_loops.bedpe")     # ten tab-separated columns ('.' / 0 fill in name / score)
P.save("loops.parquet")                # parquet; Pairs.load / as_pairs read it back
gb.Pairs.load("loops.parquet")
# -> Pairs(n=4, cis=3, trans=1, cols=[name, score])
P.to_pandas(), P.to_polars(), P.to_arrow(), P.to_numpy()   # frames, or a structured array

Pairs files: streaming raw contacts

The functions below work on contact files — one read pair per line, as written by HiC-Pro, pairtools / 4DN and Juicer — not on BEDPE loops. They stream the file once in chunks and never build a full matrix.

read_pairs_chunks(path, format='auto', columns=None, chunksize=2_000_000) yields pandas frames of chrom1, pos1, chrom2, pos2 (categorical chromosomes). The format is detected from the first data line:

format Layout Columns used (0-based)
'pairs' pairtools / 4DN .pairs: readID chrom1 pos1 chrom2 pos2 strand1 strand2 ... 1, 2, 3, 4
'allvalidpairs' HiC-Pro: readID chr1 pos1 strand1 chr2 pos2 strand2 ... 1, 2, 4, 5
'juicer' Juicer medium: readname str1 chr1 pos1 frag1 str2 chr2 pos2 frag2 ..., strands as + / - or as SAM flags 0 / 16 2, 3, 6, 7

# header lines are skipped, .gz is read. For any other layout pass columns=(chrom1, pos1, chrom2, pos2); an unknown format name raises an error listing the known ones.

from genomeblocks import bedpe
chunk = next(bedpe.read_pairs_chunks("sample.allValidPairs"))
chunk.head(2)
# ->   chrom1  pos1 chrom2  pos2
# -> 0   chr1  8482   chr2  3933
# -> 1   chr2  6000   chr1  1727
bedpe.PAIRS_FORMAT_COLUMNS
# -> {'allvalidpairs': (1, 2, 4, 5), 'pairs': (1, 2, 3, 4), 'juicer': (2, 3, 6, 7)}

Contacts per window: count_pairs

count_pairs(loci, pairs_file) counts, for every window of loci, the contacts landing in it, split by the chromosome of the partner end. Both ends of every pair are counted: each end inside a window adds one to that window under the other end’s chromosome, so a cis pair with both ends in windows counts once per end. The result is a DataFrame aligned to the rows of loci — chrom, start, end, uid, then one column per partner chromosome that has at least one contact.

PAIRS CHUNK (pandas) chrom1 pos1 chrom2 pos2 chr1 12,040 chr1 61,500 chr1 60,380 chr2 5,120 chr2 4,990 chr2 18,200 chr1 33,700 chrX 9,410 chrom columns are categorical: names → codes once per category chr1 → 0 chr2 → 1 chrX → 2 WINDOW KEYS · code · 2⁴⁰ + start (sorted) w0 w1 w2 · w4 w0 w1 w0 chr1 chr2 chrX gap 60,380 on chr1 → window 1 i = searchsorted(keys, code·2⁴⁰ + pos) − 1 keep i if same code and pos < end, else −1 COUNTS count_pairs: windows × partner chr1 chr2 chrX count_pairs_2d: windows × windows sparse CSR bincount / unique per chunk
Each read end lands in at most one window. Windows are keyed by (chromosome code, start) once; every end of every pair is placed with one searchsorted and dropped when it falls past its window's end, then the counts are a bincount per partner chromosome. That is why windows on one chromosome must not overlap, why both ends of a cis pair add to the table, and why the file streams in chunks without a matrix. 
W = gb.Loci.tile_genome({"chr1": 10_000, "chr2": 10_000}, 2500)
df = bedpe.count_pairs(W, "sample.allValidPairs")
# -> [INFO] processed 400 pairs total
df
# ->   chrom  start    end                 uid  chr1  chr2
# -> 0  chr1      0   2500      chr1:0-2500(.)    65    31
# -> 1  chr1   2500   5000   chr1:2500-5000(.)    73    23
# -> 2  chr1   5000   7500   chr1:5000-7500(.)    53    28
# -> 3  chr1   7500  10000  chr1:7500-10000(.)    69    29
# -> 4  chr2      0   2500      chr2:0-2500(.)    38    74
# -> 5  chr2   2500   5000   chr2:2500-5000(.)    33    91
# -> 6  chr2   5000   7500   chr2:5000-7500(.)    19    82
# -> 7  chr2   7500  10000  chr2:7500-10000(.)    21    71
df[["chr1", "chr2"]].to_numpy().sum()          # 400 pairs, both ends in a window: 800
# -> 800

target_chrom="chr8" keeps only the contacts whose partner is on that chromosome and returns a single count column:

W.count_pairs("sample.allValidPairs", target_chrom="chr2", verbose=False)
# ->   chrom  start    end                 uid  count
# -> 0  chr1      0   2500      chr1:0-2500(.)     31
# -> 1  chr1   2500   5000   chr1:2500-5000(.)     23
# -> 2  chr1   5000   7500   chr1:5000-7500(.)     28
# -> 3  chr1   7500  10000  chr1:7500-10000(.)     29
# -> 4  chr2      0   2500      chr2:0-2500(.)     74
# -> 5  chr2   2500   5000   chr2:2500-5000(.)     91
# -> 6  chr2   5000   7500   chr2:5000-7500(.)     82
# -> 7  chr2   7500  10000  chr2:7500-10000(.)     71

Loci.count_pairs(file, **kw) is the method form. format, columns and chunksize are passed to the reader; verbose=False silences the progress line.

Windows on one chromosome must not overlap — each position is placed in at most one window by binary search. Tiles from Loci.tile_genome or merged peaks (L.merge()) qualify; raw overlapping peaks do not.

Window-to-window matrix: count_pairs_2d

count_pairs_2d(loci_a, pairs_file, loci_b=None) returns a scipy.sparse.csr_matrix of shape (rows of loci_a) x (rows of loci_b). With loci_b=None the matrix is over loci_a on both axes and symmetric: both orientations of every pair are counted, so a pair with both ends in window i adds 2 to cell (i, i), and a pair between i and j adds 1 to (i, j) and 1 to (j, i).

LOCI_A · 4 WINDOWS ON chr1 · 3 READ PAIRS w0 w1 w2 w3 p0: w0 – w2 p1: w1 – w1 p2: w3 – w0 each end → its window by searchsorted; a pair is counted once as (i, j) and once as (j, i) scipy.sparse CSR · shape (len(loci_a), len(loci_b)) M[i, j] = reads with one end in window i, the other in j count_pairs_2d(loci_a, pairs) · 4 × 4 0 0 1 1 0 2 0 0 1 0 0 0 1 0 0 0 2 w0 w0 w1 w1 w2 w2 w3 w3 (1, 1) = 2 both ends of p1 in w1: +1 per end (0, 2) = (2, 0) = 1 both orientations of p0 symmetric: M == M.T · M.sum() = 2 × pairs loci_b GIVEN · a × b, NOT SYMMETRIC 0 2 2 0 1 0 1 0 w0 w1 w2 w3 b0 b1 b0 = [0, 200) · b1 = [200, 400) w0: p0 → b1, p2 → b1 = 2 w1: p1, both ends → b0 = 2 w2, w3: one end each → b0
count_pairs_2d counts both orientations of every pair. Three read pairs over four windows: p0 joins w0 and w2 and adds 1 to (0, 2) and 1 to (2, 0); p1 has both ends inside w1 and adds 1 per end, so (1, 1) = 2; p2 joins w3 and w0. The matrix over one set is symmetric and sums to twice the number of pairs. With loci_b given, rows are windows of loci_a and columns windows of loci_b, and nothing is symmetric: the same three pairs give w0 → b1 = 2, w1 → b0 = 2, w2 → b0 = 1 and w3 → b0 = 1. 
M = bedpe.count_pairs_2d(W, "sample.allValidPairs", verbose=False)
M.shape, M.nnz, int(M.sum()), (M != M.T).nnz
# -> ((8, 8), 64, 800, 0)
M.toarray()[:4, :4]                                    # the chr1 x chr1 block
# -> array([[18, 22,  8, 17],
# ->        [22, 18, 18, 15],
# ->        [ 8, 18, 10, 17],
# ->        [17, 15, 17, 20]])

Two helpers read the matrix back through the windows:

bedpe.pair_2d_to_frame(M, W, W).head(3)                # non-zero cells as a long frame
# ->   chrom1  start1  end1 chrom2  start2  end2  count
# -> 0   chr1       0  2500   chr1       0  2500     18
# -> 1   chr1       0  2500   chr1    2500  5000     22
# -> 2   chr1       0  2500   chr1    5000  7500      8
block, wa, wb = bedpe.pair_2d_block(M, W, W, "chr1", "chr2")   # one dense chromosome pair
block.shape, wa, wb
# -> ((4, 4), Loci(n=4, chroms=1), Loci(n=4, chroms=1))

With a second set, the matrix is loci_a x loci_b and no longer symmetric:

W.count_pairs_2d("sample.allValidPairs", loci_b=W.by_chrom("chr2"), verbose=False).shape
# -> (8, 4)

Hand-off to Architecture

A Pairs is what Architecture.make builds its edges from: every anchor is mapped to the CREs within r bp of its midpoint and each (CRE, CRE) combination becomes an edge. Pass the Pairs, a BEDPE path or a frame — as_pairs runs on it either way — so you can filter first:

cre = gb.Loci.make("peaks.bed")
loops = gb.Pairs.make("loops.bedpe").filter(min_score=2)
A = gb.Architecture.make(cre, loops, r=100)
# -> [INFO] 3 loops | 3 mapped (100.0%) | loci=6, links=3 (0 trans)
A
# -> Architecture(name='Skeleton', loci=6, links=3 [3 cis · 0 trans], edge_props=[w], vertex_props=[])

The browser draws a Pairs or a BEDPE path as loop arcs (gb.browser(..., tracks={"loops": loops})) — see Browser.

End-to-end: enhancer-promoter loops 

import genomeblocks as gb

genes = gb.Genes.make("gencode.v44.annotation.gtf", promoter_r=1000)
cre   = gb.Loci.make("atac.narrowPeak")
loops = gb.Pairs.make("loops.bedpe", min_score=1)

promoters = genes.annot["prom"]
enhancers = cre - promoters

pa, pb = loops.anchors_overlap(promoters, r=500)
ea, eb = loops.anchors_overlap(enhancers, r=500)
ep = loops[(pa & eb) | (pb & ea)]

ep.cols["gene_a"], _ = genes.nearest_tss(ep.a)      # the gene at each anchor, as a column on the pairs
ep.cols["gene_b"], _ = genes.nearest_tss(ep.b)
ep.to_bedpe("ep_loops.bedpe")
ep.to_pandas()

ep.a and ep.b are Loci, so nearest_tss labels each anchor, and the result goes into ep.cols — the extra-column dict, which to_pandas() and to_bedpe() write out (a Pairs has no P["name"] = ... shorthand; P["name"] reads a column).

See also

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

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