Geniml: Genomic Interval Machine Learning
Overview
Geniml is a Python library that bridges genomic interval biology and machine learning. It provides region2vec for learning dense vector representations of genomic regions from BED files, BEDSpace for nearest-neighbor search in embedding space, dataset classes for ML-ready genomic interval loading, and evaluation utilities for embedding quality. Geniml is designed for researchers who want to apply modern ML techniques to chromatin accessibility, histone modification, or other region-based genomic data.
When to Use
- Learn dense embeddings of genomic regions from a collection of BED files to enable ML-based analysis (region2vec)
- Cluster chromatin accessibility peaks or histone modification sites by embedding similarity
- Search for genomic regions similar to a query region using approximate nearest-neighbor search (BEDSpace)
- Build training datasets for ML models from BED-format genomic intervals with a PyTorch-compatible interface
- Compare embedding quality across training runs or datasets using quantitative metrics
- Integrate genomic region representations into custom neural network architectures
- For basic BED file parsing and set operations without ML, use
gtarsorpysam-genomic-filesinstead
Prerequisites
- Python packages:
geniml,torch,numpy,pandas,anndata - Data requirements: BED files (minimum 3 columns: chr, start, end); optionally a pre-built universe file
- Environment: Python 3.8+; GPU optional but recommended for region2vec training on large datasets
pip install geniml torch numpy pandas anndata
Quick Start
from geniml.region2vec import Region2VecExModel
from geniml.io import RegionSet
# Load a collection of BED files and train region2vec embeddings
region_sets = [RegionSet("sample1.bed"), RegionSet("sample2.bed"), RegionSet("sample3.bed")]
model = Region2VecExModel("path/to/universe.bed")
model.train(region_sets, epochs=10, batch_size=32)
# Get embedding for a specific region
embedding = model.encode("chr1", 1000000, 1500000)
print(f"Embedding shape: {embedding.shape}")
# Embedding shape: (100,)
Core API
Module 1: RegionSet — Genomic Interval I/O
Load BED files into geniml's primary data structure for downstream operations.
from geniml.io import RegionSet
# Load a BED file
rs = RegionSet("peaks.bed")
print(f"Loaded {len(rs)} regions")
print(f"First region: {rs[0]}") # Region object: chr, start, end
print(f"Chromosomes: {set(r.chr for r in rs)}")
# Convert to list of Region objects
regions = list(rs)
for r in regions[:3]:
print(f" {r.chr}:{r.start}-{r.end}")
from geniml.io import RegionSet
# Create RegionSet from a list of (chr, start, end) tuples
regions_data = [
("chr1", 100000, 101000),
("chr1", 200000, 201500),
("chr2", 50000, 51200),
]
rs = RegionSet(regions_data)
print(f"RegionSet with {len(rs)} regions from list")
# Access by index
r = rs[0]
print(f"chr={r.chr}, start={r.start}, end={r.end}, width={r.end - r.start}")
Module 2: Universe Building
A universe defines the set of consensus regions used as the vocabulary for region2vec. Build it from a collection of BED files.
from geniml.universe import UniverseBuilder
from geniml.io import RegionSet
# Collect BED files representing diverse samples
bed_files = ["sample1.bed", "sample2.bed", "sample3.bed", "sample4.bed"]
region_sets = [RegionSet(f) for f in bed_files]
# Build universe (consensus non-overlapping regions)
builder = UniverseBuilder()
universe = builder.build(region_sets)
# Save universe to BED file
universe.to_bed("universe.bed")
print(f"Universe size: {len(universe)} consensus regions")
from geniml.universe import UniverseBuilder
from geniml.io import RegionSet
# Build universe with custom parameters
builder = UniverseBuilder(
fraction=0.5, # Region must appear in >= 50% of samples to be included
merge_dist=0, # Merge adjacent regions within this distance (bp)
)
bed_files = [f"sample_{i}.bed" for i in range(1, 11)]
region_sets = [RegionSet(f) for f in bed_files]
universe = builder.build(region_sets)
print(f"Filtered universe: {len(universe)} regions (fraction >= 0.5)")
Module 3: Region2Vec — Training Embeddings
Train word2vec-style embeddings on genomic regions, treating each BED file as a "document" and each region as a "word."
from geniml.region2vec import Region2VecExModel
from geniml.io import RegionSet
# Initialize model with a pre-built universe
model = Region2VecExModel(universe="universe.bed", embedding_dim=100)
# Load training data (collection of BED files = corpus)
bed_files = [f"atac_{i}.bed" for i in range(1, 51)]
region_sets = [RegionSet(f) for f in bed_files]
# Train
model.train(
region_sets,
epochs=20,
batch_size=64,
window_size=5,
min_count=1,
)
print("Training complete")
# Save trained model
model.save("region2vec_model/")
print("Model saved to region2vec_model/")
from geniml.region2vec import Region2VecExModel
# Load a pre-trained model
model = Region2VecExModel.load("region2vec_model/")
# Encode a single genomic region
embedding = model.encode("chr1", 1_000_000, 1_500_000)
print(f"Single region embedding shape: {embedding.shape}")
# Single region embedding shape: (100,)
# Encode an entire BED file → matrix of region embeddings
from geniml.io import RegionSet
rs = RegionSet("query_peaks.bed")
embeddings = model.encode_region_set(rs)
print(f"BED file embeddings shape: {embeddings.shape}")
# BED file embeddings shape: (N_regions, 100)
Module 4: BEDSpace — Embedding Nearest-Neighbor Search
Index a corpus of BED file embeddings for fast similarity search.
from geniml.bedspace import BEDSpace
from geniml.region2vec import Region2VecExModel
from geniml.io import RegionSet
# Build a BEDSpace index from a set of BED files
model = Region2VecExModel.load("region2vec_model/")
bed_files = [f"dataset_{i}.bed" for i in range(1, 101)]
region_sets = [RegionSet(f) for f in bed_files]
bedspace = BEDSpace(model)
bedspace.fit(region_sets, labels=[f"dataset_{i}" for i in range(1, 101)])
bedspace.save("bedspace_index/")
print(f"BEDSpace index built with {len(region_sets)} datasets")
from geniml.bedspace import BEDSpace
from geniml.io import RegionSet
# Load index and query
bedspace = BEDSpace.load("bedspace_index/")
query = RegionSet("query_sample.bed")
# Find top-k most similar datasets
results = bedspace.query(query, k=5)
for rank, (label, score) in enumerate(results, 1):
print(f" Rank {rank}: {label} similarity={score:.4f}")
Module 5: Genomic Interval Datasets for ML
PyTorch-compatible Dataset classes for training ML models on genomic intervals.
from geniml.datasets import TokenizedBEDDataset
from torch.utils.data import DataLoader
# Create a tokenized dataset from multiple BED files + labels
bed_files = ["condition_A_rep1.bed", "condition_A_rep2.bed",
"condition_B_rep1.bed", "condition_B_rep2.bed"]
labels = [0, 0, 1, 1]
dataset = TokenizedBEDDataset(
bed_files=bed_files,
universe="universe.bed",
labels=labels,
)
print(f"Dataset size: {len(dataset)} samples")
# Use with PyTorch DataLoader
loader = DataLoader(dataset, batch_size=8, shuffle=True)
for batch_tokens, batch_labels in loader:
print(f"Batch tokens shape: {batch_tokens.shape}")
print(f"Batch labels: {batch_labels}")
break
from geniml.datasets import RegionEmbeddingDataset
from geniml.region2vec import Region2VecExModel
import torch
# Pre-embed BED files and create a dataset for a classifier
model = Region2VecExModel.load("region2vec_model/")
bed_files = ["pos_1.bed", "pos_2.bed", "neg_1.bed", "neg_2.bed"]
labels = [1, 1, 0, 0]
emb_dataset = RegionEmbeddingDataset(
bed_files=bed_files,
model=model,
labels=labels,
aggregation="mean", # aggregate region embeddings per sample
)
X, y = emb_dataset[0]
print(f"Sample embedding shape: {X.shape}, label: {y}")
# Sample embedding shape: (100,), label: 1
Module 6: Embedding Evaluation
Assess embedding quality using neighborhood overlap and intrinsic metrics.
from geniml.eval import EmbeddingEvaluator
from geniml.region2vec import Region2VecExModel
from geniml.io import RegionSet
model = Region2VecExModel.load("region2vec_model/")
bed_files = [f"sample_{i}.bed" for i in range(1, 21)]
region_sets = [RegionSet(f) for f in bed_files]
evaluator = EmbeddingEvaluator(model)
metrics = evaluator.evaluate(region_sets)
print(f"Neighborhood overlap score: {metrics['neighborhood_overlap']:.4f}")
print(f"Silhouette score: {metrics.get('silhouette', 'N/A')}")
# Higher neighborhood overlap → more biologically coherent embeddings
Key Concepts
region2vec Training Analogy
Region2vec treats each BED file as a "sentence" (a document of co-occurring genomic regions) and each genomic region token (from the universe) as a "word." Word2vec's skip-gram objective is applied: regions that frequently co-occur in BED files are learned to have similar embeddings. This means two regions that tend to be open/active in the same set of samples will have nearby embeddings, enabling meaningful similarity search without explicit labels.
# Conceptual analogy: universe regions = vocabulary tokens
# BED file = sentence = [region_token_1, region_token_2, ...]
# Training: predict co-occurring regions within a sliding window
from geniml.io import RegionSet
rs = RegionSet("sample.bed")
print(f"This BED file contains {len(rs)} 'words' (tokens) from the universe vocabulary")
Universe as Vocabulary
The universe is a non-overlapping, genome-wide set of consensus regions that serves as the token vocabulary for region2vec. A well-chosen universe should cover the genomic regions present in your dataset while being compact enough for efficient training. Regions in a BED file that do not overlap the universe are ignored during training.
Common Workflows
Workflow 1: Train and Search Embeddings from ATAC-seq Peaks
Goal: Embed ATAC-seq peak sets, then find samples most similar to a query.
from geniml.universe import UniverseBuilder
from geniml.region2vec import Region2VecExModel
from geniml.bedspace import BEDSpace
from geniml.io import RegionSet
from pathlib import Path
# Step 1: Collect BED files
bed_dir = Path("atac_peaks/")
bed_files = sorted(bed_dir.glob("*.bed"))
region_sets = [RegionSet(str(f)) for f in bed_files]
labels = [f.stem for f in bed_files]
print(f"Loaded {len(region_sets)} ATAC-seq peak sets")
# Step 2: Build universe
builder = UniverseBuilder(fraction=0.3)
universe = builder.build(region_sets)
universe.to_bed("atac_universe.bed")
print(f"Universe: {len(universe)} regions")
# Step 3: Train region2vec
model = Region2VecExModel(universe="atac_universe.bed", embedding_dim=100)
model.train(region_sets, epochs=15, batch_size=64, window_size=5)
model.save("atac_region2vec/")
# Step 4: Build BEDSpace index
bedspace = BEDSpace(model)
bedspace.fit(region_sets, labels=labels)
bedspace.save("atac_bedspace/")
# Step 5: Query with a new sample
query = RegionSet("new_sample.bed")
results = bedspace.query(query, k=5)
print("Top 5 most similar samples:")
for rank, (label, score) in enumerate(results, 1):
print(f" {rank}. {label} (similarity={score:.4f})")
Workflow 2: ML Classification on Genomic Intervals
Goal: Train a logistic regression classifier on region2vec embeddings to classify cell types.
from geniml.region2vec import Region2VecExModel
from geniml.datasets import RegionEmbeddingDataset
from geniml.io import RegionSet
from torch.utils.data import DataLoader
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
import numpy as np
model = Region2VecExModel.load("atac_region2vec/")
# Prepare labeled dataset (e.g., T-cell vs B-cell ATAC peaks)
bed_files = ["tcell_1.bed", "tcell_2.bed", "tcell_3.bed",
"bcell_1.bed", "bcell_2.bed", "bcell_3.bed"]
labels = [0, 0, 0, 1, 1, 1]
dataset = RegionEmbeddingDataset(
bed_files=bed_files, model=model, labels=labels, aggregation="mean"
)
# Collect embeddings
X = np.array([dataset[i][0].numpy() for i in range(len(dataset))])
y = np.array([dataset[i][1] for i in range(len(dataset))])
print(f"Feature matrix: {X.shape}, labels: {y}")
# Cross-validate a simple classifier
clf = LogisticRegression(max_iter=1000)
scores = cross_val_score(clf, X, y, cv=3, scoring="accuracy")
print(f"Cross-val accuracy: {scores.mean():.3f} ± {scores.std():.3f}")
Key Parameters
| Parameter | Module | Default | Range / Options | Effect |
|---|---|---|---|---|
embedding_dim |
Region2Vec | 100 |
50–512 |
Dimensionality of region embeddings; higher = more expressive but slower |
window_size |
Region2Vec training | 5 |
2–20 |
Context window for skip-gram co-occurrence; larger = broader co-occurrence |
epochs |
Region2Vec training | 10 |
5–100 |
Training iterations; increase for large datasets |
batch_size |
Region2Vec training | 32 |
16–512 |
Mini-batch size; larger batches improve stability |
fraction |
UniverseBuilder | 0.5 |
0.0–1.0 |
Minimum fraction of samples a region must appear in to enter the universe |
aggregation |
RegionEmbeddingDataset | "mean" |
"mean", "sum", "max" |
How per-region embeddings are pooled into a sample-level vector |
k |
BEDSpace query | 10 |
1–corpus size |
Number of nearest neighbors to return |
Common Recipes
Recipe: Load a Pre-Trained Hugging Face region2vec Model
When to use: Apply a community-trained model without local training.
from geniml.region2vec import Region2VecExModel
# Load a pre-trained model from HuggingFace Hub
model = Region2VecExModel("databio/r2v-ChIP-atlas-hg38-v2")
print(f"Model embedding dim: {model.embedding_dim}")
from geniml.io import RegionSet
rs = RegionSet("my_peaks.bed")
embeddings = model.encode_region_set(rs)
print(f"Embeddings shape: {embeddings.shape}")
Recipe: UMAP Visualization of Sample Embeddings
When to use: Visually inspect clustering of samples in embedding space.
from geniml.region2vec import Region2VecExModel
from geniml.io import RegionSet
import numpy as np
model = Region2VecExModel.load("atac_region2vec/")
bed_files = [f"sample_{i}.bed" for i in range(1, 21)]
sample_embeddings = []
for f in bed_files:
rs = RegionSet(f)
emb = model.encode_region_set(rs)
sample_embeddings.append(emb.mean(axis=0)) # mean-pool regions
X = np.array(sample_embeddings)
print(f"Sample embedding matrix: {X.shape}")
try:
import umap
import matplotlib.pyplot as plt
reducer = umap.UMAP(n_components=2, random_state=42)
X_2d = reducer.fit_transform(X)
plt.figure(figsize=(6, 5))
plt.scatter(X_2d[:, 0], X_2d[:, 1], s=40)
for i, f in enumerate(bed_files):
plt.annotate(f, (X_2d[i, 0], X_2d[i, 1]), fontsize=7)
plt.title("Sample embeddings (UMAP)")
plt.savefig("sample_umap.png", dpi=150, bbox_inches="tight")
print("Saved sample_umap.png")
except ImportError:
print("Install umap-learn for UMAP visualization")
Troubleshooting
| Problem | Cause | Solution |
|---|---|---|
FileNotFoundError on universe |
Universe BED path incorrect or not built yet | Run UniverseBuilder.build() and save with .to_bed() before training |
| All embeddings identical or near-zero | Universe has no overlap with training BED files | Verify BED coordinate system matches universe (hg38 vs hg19, chr prefix) |
RuntimeError: CUDA out of memory |
Batch size too large for GPU | Reduce batch_size or use CPU training |
| Very low neighborhood overlap score | Too few training samples or too many epochs (overfitting) | Use ≥20 BED files; tune epochs; try smaller embedding_dim |
| BEDSpace query returns no results | Query regions not present in index universe | Check query BED overlaps universe; use fraction parameter to broaden universe |
| Slow training | Large universe + many BED files | Reduce universe size with higher fraction threshold; use GPU |
Related Skills
- gtars — fast BED file I/O and interval operations for preprocessing before geniml training
- scanpy-scrna-seq — downstream clustering and UMAP visualization of embeddings stored in AnnData
References
- Geniml documentation — official docs and tutorials
- Geniml GitHub repository — source code and examples
- PyPI: geniml — installation and version info
- HuggingFace: databio region2vec models — pre-trained region2vec model collection