How do I actually build a GRanges object from scratch?

If you've spent any time in bioinformatics, you know that genomic coordinates are a nightmare to track in a standard data frame. You've got chromosomes, start positions, end positions, and the strand—and if you get one wrong, your whole analysis is garbage. That's why we use GRanges from the GenomicRanges package. It's not just a list; it's a specialized object that understands the "geometry" of a genome.

I usually start by defining my coordinates explicitly. Let's say we're looking at a few specific regions on Chromosome 17, perhaps around the TP53 gene. Here is how I'd set that up:

library(GenomicRanges)

# Let's define some regions of interest
my_ranges <- GRanges(
    seqnames = "chr17", 
    ranges = IRanges(start = c(7661779, 7686769), end = c(7686448, 7686915)), 
    strand = c("+", "-"), 
    name = c("region_A", "region_B")
)

my_ranges

Notice that I used IRanges inside the GRanges call. Think of IRanges as the "integer range" (the start and end) and GRanges as the "genomic range" (the integer range + the chromosome + the strand). I've also added a name column. You can add almost any metadata you want here, and it stays glued to the coordinates.

How do I find overlaps between two different genomic datasets?

This is where GRanges really earns its keep. Imagine you have one set of ranges representing your "peaks" from a ChIP-seq experiment and another set representing known gene promoters. You want to know which peaks hit which promoters. In a standard data frame, you'd be writing nested loops or complex joins. In GenomicRanges, you just use findOverlaps.

# Dataset 1: My experimental peaks
peaks <- GRanges(seqnames = "chr17", ranges = IRanges(start = c(7660000, 7670000), end = c(7665000, 7675000)))

# Dataset 2: Known promoter regions
promoters <- GRanges(seqnames = "chr17", ranges = IRanges(start = c(7661000, 7680000), end = c(7662000, 7681000)))

# Find the overlaps
hits <- findOverlaps(peaks, promoters)

# 'hits' is a Hits object, which is basically a mapping of indices.
# To see which peak index matches which promoter index:
hits

The Hits object is incredibly memory efficient, which is vital when you're dealing with millions of reads. If you want to actually extract the overlapping ranges, you can use the hits to subset your objects, like peaks[queryHits(hits)]. I've found that most people struggle with the Hits object at first, but once you realize it's just a map of "Index A matches Index B," it becomes second nature.

What does "reducing" actually do to my ranges?

You'll see the reduce() function everywhere in genomic pipelines. In plain English: reduce takes a set of overlapping or abutting ranges and merges them into a single, larger range. It's basically "collapsing" the intervals.

I often use this when I have a bunch of fragmented reads that I know belong to the same biological feature. I don't care about the individual gaps; I just want the total footprint of the activity.

# Create some overlapping ranges
messy_ranges <- GRanges(
    seqnames = "chr17", 
    ranges = IRanges(start = c(100, 150, 500), end = c(200, 250, 600))
)

# Merge them
clean_ranges <- reduce(messy_ranges)

# messy_ranges had (100-200) and (150-250). 
# clean_ranges will have (100-250) and (500-600).
clean_ranges

Just a heads-up: reduce() will merge ranges on the same chromosome and strand. If you have ranges on different strands, they won't be merged even if they overlap perfectly. That's usually what you want, because the strand carries critical biological meaning.

📋 Practical Task

Exercise: Merging and Filtering ChIP-seq Peak Intervals

You have been given a set of raw ChIP-seq peak coordinates on Chromosome 1. Some of these peaks overlap because the peak-calling software was too sensitive. Your goal is to merge these overlapping regions into "consolidated peaks" and then filter out any consolidated peak that is shorter than 200 base pairs (as these are likely noise).

Instructions:

1. Create a GRanges object called raw_peaks with the following coordinates on "chr1":
   • Start: 1000, End: 1200
   • Start: 1100, End: 1300
   • Start: 5000, End: 5100
   • Start: 8000, End: 8500
   • Start: 8400, End: 8600
2. Use the reduce() function to merge overlapping intervals into a new object called merged_peaks.
3. Calculate the width of these merged peaks using the width() function.
4. Filter merged_peaks to keep only those with a width ≥ 200.
5. Print the final filtered GRanges object.