I was chatting with a colleague recently who was trying to pull some RNA-seq analysis into R. They kept complaining that they couldn't find the DESeq2 package on CRAN. When I looked at their screen, they were just hammering away at install.packages("DESeq2") and getting a "package not available" error. It reminded me of the first time I hit this wall. You assume that if it's a famous R package, it's on CRAN. But biological data is a different beast entirely.

The CRAN wall

Let's try to reproduce exactly what my colleague did. If I just try to install a heavy-hitter biological package the standard way, here is what happens:

install.packages("DESeq2")
# Error in install.packages("DESeq2") : unable to find available package DESeq2

Right. It's not there. Now, I could spend an hour searching Google, but if you're doing bioinformatics in R, you have to realize that the community essentially built their own parallel universe called Bioconductor. It's not just a collection of packages; it's a whole ecosystem with its own repository and its own versioning system that stays in sync with specific versions of R.

Bridging the gap with BiocManager

To get into this ecosystem, we can't use the default R installation tool. We need a manager that knows how to talk to the Bioconductor servers. I'll start by installing BiocManager—which, ironically, *is* on CRAN—and then use it to grab the biological tools.

install.packages("BiocManager")
library(BiocManager)

# Now I'll try that DESeq2 install again, but through the manager
BiocManager::install("DESeq2")

You'll notice during the installation that BiocManager does something install.packages doesn't: it checks your R version and ensures you're getting the version of the package that is compatible with your specific R release. This is a lifesaver because biological packages often depend on very specific C++ libraries that break if the versions are mismatched.

Not your average data frame

Once I have a Bioconductor package loaded, I usually run into the same second hurdle: the data structures. I'm going to load GenomicRanges, another staple of the ecosystem, and try to create a simple genomic interval. I'll treat it like a data frame because that's how we do everything in base R, right?

BiocManager::install("GenomicRanges")
library(GenomicRanges)

# Let's define a region on chromosome 1
my_region <- GRanges(seqnames = "chr1", 
                     ranges = IRanges(start = 100, end = 200))

print(my_region)
# GRanges object with 1 range
# genomic ranges: 1
# seqnames ranges start end width strand
# chr1    1       100 200     101  *

Wait a minute. If I try to access the start position using my_region$start, it returns NULL. That's annoying. Why? Because Bioconductor doesn't rely heavily on the standard S3 objects (like data frames and lists) that we've used in previous lessons. They use S4 objects.

S4 is a more formal system. It's stricter. It prevents you from accidentally shoving a character string into a column that must be an integer. To get data out of these objects, you often need specific "accessor" functions rather than the $ operator.

# This doesn't work
# my_region$start 

# This is the Bioconductor way
start(my_region)
# [1] 100

end(my_region)
# [1] 200

It feels a bit clunky at first, but when you're dealing with a genome that has billions of base pairs, that strictness is what keeps your analysis from collapsing into a heap of type-conversion errors.

📋 Practical Task

Mapping Gene Coordinates with GRanges

In this exercise, you will practice navigating the Bioconductor ecosystem by creating a genomic range object and extracting specific metadata using S4 accessor functions.

• First, ensure you have BiocManager installed and use it to install the GenomicRanges package.
• Create a GRanges object named gene_coords that contains two ranges:
  • Range 1: chr1, start 1000, end 2000
  • Range 2: chr2, start 5000, end 6000
• Add a metadata column to this object called gene_name with the values "Gene_A" and "Gene_B".
• Using the correct Bioconductor accessor functions, extract the widths of these two ranges and the seqnames (chromosome names) and print them to the console.