# Genomic Ranges The `GenomicRanges` package provides data structures for representing and manipulating genomic intervals. It is the cornerstone of the BiocPy ecosystem, allowing you to assign genomic coordinates (chromosome, start, end, strand) to data. ## The Data Model At its core, a `GenomicRanges` object is a vector-like container where each element represents a genomic interval. It is built upon simpler primitives. ### 1. IRanges: The Foundation Before understanding `GenomicRanges`, we must look at `IRanges`. An `IRanges` object represents integer intervals defined by a **start** position and a **width**. * **1-based indexing:** BiocPy follows the Bioconductor/R convention for intervals (1-based, closed). * **Closed intervals:** `[start, end]` are both included. * **Relationship:** `end = start + width - 1`. ```{python} from iranges import IRanges # Construct IRanges ir = IRanges(start=[1, 10, 20], width=[5, 5, 5]) # Accessors print(f"Starts: {ir.start}") print(f"Ends: {ir.end}") print(f"Widths: {ir.width}") ``` ### 2. GenomicRanges: Adding Context An `IRanges` tells you *where* (coordinate-wise), but not *on what* (chromosome/sequence). `GenomicRanges` extends `IRanges` by adding: * **Seqnames:** The chromosome or scaffold name (e.g., "chr1"). * **Strand:** The direction of the feature (`+`, `-`, or `*`). * **Metadata (mcols):** Arbitrary annotations like scores, gene IDs, or names. ```{mermaid} classDiagram class GenomicRanges { +seqnames: vector +strand: vector +mcols: DataFrame } class IRanges { +start: vector +width: vector } GenomicRanges *-- IRanges : contains ``` If you represent regions across a genome, you use `GenomicRanges`. If you represent regions within a single unstranded sequence (like a protein or specific transcript), you might use `IRanges` directly. ## Core Concepts Summary * **Seqnames:** "chr1", "chrX", etc. * **Ranges:** The underlying `IRanges` (start, width). * **Strand:** * `+` (Forward/Watson) * `-` (Reverse/Crick) * `*` (Unstranded/Unknown) **Metadata:** A `BiocFrame` holding any other data columns aligned to the ranges. --- # Robust Genomic Integration One of the key advantages of BiocPy is the ability to leverage the vast annotation resources of Bioconductor without managing flat files yourself. The `txdb` package allows you to query standard Bioconductor `TxDb` SQLite databases and populate `GenomicRanges` objects directly. This bypasses the need to parse error-prone GTF/GFF files and ensures you are working with curated, versioned gene models. See [Accessing Annotations](../biochubs/annotations.qmd) for more details on `TxDb` usage. By loading these annotations into `GenomicRanges`, you can immediately leverage interval arithmetic for analysis, as shown in the workflow below. --- # Workflow: ChIP-seq Peak Annotation To demonstrate the power of `GenomicRanges`, let's walk through a common bioinformatics task: **annotating ChIP-seq peaks**. **Goal:** Given a set of binding peaks and gene annotations, identify which peaks fall within "promoter" regions (upstream of genes). ## 1. Setup Mock Data First, we generate mock data representing our "Genes" and our "Peaks". ```{python} import numpy as np from genomicranges import GenomicRanges from iranges import IRanges from biocframe import BiocFrame # --- 1. Create Gene Annotations --- # Imagine 5 genes on chr1 and chr2 genes = GenomicRanges( seqnames=["chr1", "chr1", "chr2", "chr2", "chr2"], ranges=IRanges([1000, 5000, 1000, 4000, 8000], width=[500]*5), # 500bp genes strand=["+", "-", "+", "-", "+"], mcols=BiocFrame({ "gene_id": ["Gene_A", "Gene_B", "Gene_C", "Gene_D", "Gene_E"], "type": ["protein_coding"] * 5 }) ) # --- 2. Create ChIP-seq Peaks --- # Imagine 6 peaks found by an uncited peak caller peaks = GenomicRanges( seqnames=["chr1", "chr1", "chr1", "chr2", "chr2", "chr3"], ranges=IRanges([900, 1200, 4800, 950, 4100, 10000], width=[50]*6), # 50bp peaks strand=["*"] * 6, # Peaks are usually unstranded mcols=BiocFrame({ "peak_id": [f"peak_{i}" for i in range(1, 7)], "score": [10.5, 8.2, 15.0, 5.5, 20.1, 3.0] }) ) print("Genes:\n", genes) print("\nPeaks:\n", peaks) ``` ## 2. Defining Promoters We define a promoter as the region **200bp upstream** of a gene's transcription start site (TSS). * For `+` strand genes: Upstream is `start - 200`. * For `-` strand genes: Upstream is `end + 200`. The `flank()` method handles this strand-aware logic automatically. ```{mermaid} graph TD A[Gene (+)] -- flank(start=True) --> B[Promoter] B --> C[Upstream Region] D[Gene (-)] -- flank(start=True) --> E[Promoter] E --> F[Downstream in coordinates\n(Upstream biologically)] ``` ```{python} # Define promoters: 200bp upstream of the TSS (start=True) promoters = genes.flank(width=200, start=True) # Add the gene_id to the promoter metadata for tracking promoters.mcols["gene_id"] = genes.mcols["gene_id"] print("Promoters:\n", promoters) ``` Notice how `Gene_B` (on `chr1`, `-` strand, original range `5000-5499`) has a promoter at `5500-5699`—biologically upstream on the negative strand. ## 3. Finding Overlaps Now we identify which peaks overlap with these defined promoters. ```{python} # Find overlaps: query=peaks, subject=promoters hits = peaks.find_overlaps(promoters) print(hits) ``` `hits` is a collection of indices indicating matches. We can use these indices to extract the matching records. ## 4. Annotating the Peaks We can now merge the peak information with the gene information. ```{python} # Indices of the overlapping peaks and promoters peak_indices = hits.get_column("query_hits") promoter_indices = hits.get_column("self_hits") # Extract the actual objects overlapping_peaks = peaks[peak_indices] overlapping_promoters = promoters[promoter_indices] # Create a summary result result = BiocFrame({ "peak_id": overlapping_peaks.mcols["peak_id"], "score": overlapping_peaks.mcols["score"], "overlaps_gene": overlapping_promoters.mcols["gene_id"], "promoter_loc": list(str(x) for x in overlapping_promoters) }) print("Annotated Peaks:\n", result) ``` ## 5. Statistical Summary How many peaks fell into promoters vs intergenic regions? ```{python} total_peaks = len(peaks) annotated_count = len(set(peak_indices)) # Use set to count unique peaks intergenic_count = total_peaks - annotated_count import pandas as pd df = pd.DataFrame({ "Category": ["Promoter", "Intergenic"], "Count": [annotated_count, intergenic_count] }) print(df) ``` --- # Advanced Operations ## Nearest Neighbors Sometimes a peak doesn't overlap a promoter directly, but you want to know the *closest* gene. ```{python} # Find the nearest gene for every peak # Returns indices in 'genes' that are closest to each 'peak' nearest_inds = genes.nearest(peaks) print("Nearest Gene Indices:", nearest_inds) # Add nearest gene info to peaks peaks.mcols["nearest_gene"] = [genes.mcols["gene_id"][int(x)] if x is not None else None for x in nearest_inds] print("\nPeaks with Nearest Genes:\n", peaks) ``` ## Set Operations Standard set operations allow you to manipulate genomic territory. * `reduce()`: Collapse overlapping regions into single contiguous blocks (e.g., getting total coverage footprint). * `setdiff()`: Subtract regions (e.g., remove blacklisted regions). ```{python} # Example: What is the total genomic footprint of our genes? footprint = genes.reduce() print("Gene Footprint:\n", footprint) # Example: Remove the promoter regions from the peaks (finding only non-promoter binding) non_promoter_peaks = peaks.setdiff(promoters) print("\nNon-promoter Binding Sites:\n", non_promoter_peaks) ``` ## Summary `GenomicRanges` provides a rich grammar for genomic arithmetic. By combining `flank`, `find_overlaps`, `nearest`, and set operations, you can build complex analysis pipelines purely in Python.