1 Getting started with Bioconductor ecosystem

Aims

Importing genomic files (coverage tracks and features) into R
Discovering the Rle / RleList representation of genome-wide coverage
Manipulating GRanges objects (resizing, centering, filtering)
Extracting and aggregating a coverage signal over a set of features
Prototyping the first three functions of the package we will build this week: importFiles(), filterGRanges() and computeCoverage()

Tip

At any time, if you are lost or do not understand how functions in the proposed solution work, type ?<function> in the R console and a help menu will appear.

You can also check the help tab in the corresponding quadrant.

Where we are heading

Over the week, we will build a package which can plot the aggregated coverage of a genomic track over a set of GRanges of interest (at a fixed width).

Today, we prototype the core steps interactively, in a plain R script. Over the next days, we will progressively turn this prototype into a proper, documented, tested Bioconductor package. By the end of today, you will have written draft versions of three of the four functions that make up the package:

importFiles(bw_file, features_file, width) |>
    filterGRanges() |>
    computeCoverage() |>
    plotCoverage()

1.1 Set up your working environment

Let’s create a project for this course! In RStudio: File > New Project..., create a project entitled “Bioc-workshop”. Open the newly created project. Download files required for the workshop from the following GitHub repository:

https://github.com/js2264/Bioc-workshop

For today, we will work with two files stored in Share/data/:

Scc1-vs-input.bw: a coverage track (bigwig format) of a cohesin (Scc1) ChIP-seq experiment in yeast (S. cerevisiae);
Scc1-peaks.narrowPeak: the set of Scc1 binding sites (genomic features), as called from the same experiment.

These are the very same files we will later ship inside our package, so getting comfortable with them now will pay off all week.

1.2 Importing a coverage track in R

A coverage track stores, for each position of the genome, how much signal (e.g. how many sequencing reads) was observed there. bigwig is the standard binary format to store such tracks.

Question

Find the appropriate Bioconductor function to import a bigwig file in R. Import the coverage track as an RleList object (as = "RleList").

Hint

The rtracklayer package provides a generic import() function that recognises many genomic file formats. Its as = argument controls the class of the returned object.

library(rtracklayer)
bw_file <- "_____"
coverage <- import(______, as = ______)

Answer

library(rtracklayer)
bw_file <- "Share/data/Scc1-vs-input.bw"
coverage <- import(bw_file, as = "RleList")
coverage

Question

What kind of object is a RleList? How many elements does it contain, and what do they correspond to? What does the Rle (“Run-Length Encoding”) representation buy us here?

Hint

An RleList is a list of Rle vectors. Try class(coverage), length(coverage), names(coverage) and inspect one element with coverage[[1]]. A genome is mostly made of long stretches of identical values — think about how run-length encoding stores 0 0 0 0 0 3 3 3 compactly.

names(_____)
_____[[1]]

Answer

class(coverage)

## One Rle per chromosome
length(coverage)
names(coverage)

## Each Rle stores runs of identical coverage values, instead of one value
## per base pair. This is extremely memory-efficient for genomic coverage,
## which is piecewise-constant.
coverage[[1]]
lengths(coverage)  ## the length of each Rle = the length of each chromosome

1.3 Importing genomic features in R

Genomic features (peaks, genes, regulatory elements, …) are typically stored as bed-like files. A narrowPeak file is a bed variant used for ChIP-seq peaks.

Question

Import the Scc1-peaks.narrowPeak file. What class of object do you obtain?

Hint

The same rtracklayer::import() function handles bed/narrowPeak files; here you do not need the as = argument. The result is a GRanges.

features <- ______("Share/data/______")

Answer

features <- import("Share/data/Scc1-peaks.narrowPeak")
features
class(features)

Question

How many features are there? On which chromosomes do they sit? What is their median width?

Hint

GRanges behave a lot like vectors. Useful accessors: length(), seqnames(), width().

length(______)
table(______(features))
median(______(features))

Answer

length(features)
table(seqnames(features))
median(width(features))

1.4 Resizing features to a fixed width

To compare and aggregate the coverage across many features, every feature must span the same number of positions. We therefore resize all features to a fixed width, centered on the middle of each feature.

Question

Resize every feature to a window of width = 2000 bp, centered on the feature center. Which GenomicRanges function does this?

Hint

Look at ?GenomicRanges::resize. Its fix = argument decides which anchor point is kept fixed while resizing: "start", "end" or "center".

width <- 2000
resized <- resize(______, width = width, fix = ______)

Answer

width <- 2000
resized <- resize(features, width = width, fix = "center")
resized
table(width(resized))  ## all features now span exactly 2000 bp

1.5 Wrap it up: your first function, `importFiles()`

You now have all the pieces to import a coverage track and a set of (fixed-width) features together. Let’s package these steps into a reusable function — the first function of our package.

Question

Write a function importFiles(bw_file, features_file, width) that:

imports the coverage track as an RleList;
imports the features as a GRanges;
resizes the features to width, centered on their center;
returns both in a named list (coverage and features).

Hint

Start from the skeleton below and fill in the # TODOs with the code you just wrote interactively. The function should simply glue those steps together.

importFiles <- function(bw_file, features_file, width) {
    coverage <- # TODO: import bw_file as an RleList
    features <- # TODO: import features_file as a GRanges
    features <- # TODO: resize features to `width`, centered
    list(coverage = ______, features = ______)
}

Answer

importFiles <- function(bw_file, features_file, width) {
    coverage <- rtracklayer::import(bw_file, as = "RleList")
    features <- rtracklayer::import(features_file)
    features <- GenomicRanges::resize(features, width = width, fix = "center")
    list(coverage = coverage, features = features)
}

l <- importFiles(
    "Share/data/Scc1-vs-input.bw",
    "Share/data/Scc1-peaks.narrowPeak",
    width = 2000
)
l

1.6 Extracting the coverage over features

We can now extract the coverage track values sitting under each feature. Subsetting an RleList by a GRanges returns, for each feature, the Rle of coverage values spanning it.

Question

Extract the coverage over the resized features. What object do you get, and how long is each of its elements?

Hint

RleList objects can be subset by a GRanges, i.e. coverage[features]. The result is another RleList, with one element per feature.

cov_by_feature <- l$______[l$______]

Answer

cov_by_feature <- l$coverage[l$features]
cov_by_feature
## Each element spans `width` positions (here 2000)
head(lengths(cov_by_feature))

1.7 Filtering features: your second function, `filterGRanges()`

Some features may sit too close to a chromosome end: once resized, their window runs past the chromosome boundary. Others may sit on a chromosome that is absent from the coverage track. In both cases, extracting the coverage would fail or return truncated profiles.

Question

Before extracting coverage, we should remove such problematic features. Write a function filterGRanges(l) that takes the list returned by importFiles() and keeps only the features that (a) sit on a chromosome present in the coverage track, and (b) are fully contained within their chromosome.

Hint

The length of each chromosome is lengths(l$coverage) (a named vector). A feature is safe to keep if its start() is >= 1 and its end() is <= the length of its chromosome. Use %in% names(l$coverage) to check the chromosome is covered.

filterGRanges <- function(l) {
    features <- l$features
    chroms <- as.character(seqnames(features))
    chrom_lengths <- lengths(l$coverage)
    keep <- # TODO: on a covered chromosome AND within bounds
    l$features <- features[keep]
    l
}

Answer

filterGRanges <- function(l) {
    features <- l$features
    chroms <- as.character(seqnames(features))
    chrom_lengths <- lengths(l$coverage)
    on_covered_chrom <- chroms %in% names(l$coverage)
    within_bounds <- start(features) >= 1 &
        end(features) <= chrom_lengths[chroms]
    l$features <- features[on_covered_chrom & within_bounds]
    l
}

l <- filterGRanges(l)
length(l$features)  ## possibly fewer features than before

1.8 Aggregating the coverage: your third function, `computeCoverage()`

We finally turn the per-feature coverage into a single, aggregated profile: for each position along the window, the mean coverage across all features, together with a confidence interval.

Question

Reshape the extracted coverage into a features x positions numeric matrix. Why a matrix, and why byrow = TRUE?

Hint

Unlisting the RleList gives one long numeric vector (feature 1’s positions, then feature 2’s, …). Folding it back into a matrix with ncol = width and byrow = TRUE puts one feature per row and one position per column — exactly what we need to average column by column.

mat <- matrix(
    as.numeric(unlist(l$coverage[l$features])),
    ncol = ______, byrow = ______
)

Answer

mat <- matrix(
    as.numeric(unlist(l$coverage[l$features])),
    ncol = width, byrow = TRUE
)
dim(mat)  ## nrow = number of features, ncol = width

Question

Now wrap everything into a function computeCoverage(l) that returns a data.frame with one row per position and the columns distance (from -width/2 to width/2 - 1), mean, ci_low and ci_high (bounds of the 95% confidence interval).

Hint

Per-position mean is colMeans(mat). For a 95% CI, use 1.96 * sd / sqrt(n), where sd is apply(mat, 2, sd) and n = nrow(mat). The distance column is seq(-width/2, width/2 - 1, by = 1).

computeCoverage <- function(l) {
    features <- l$features
    width <- unique(width(features))
    mat <- # TODO: build the features x positions matrix
    means <- # TODO: per-position mean
    ci <- # TODO: 1.96 * per-position sd / sqrt(n)
    data.frame(
        distance = seq(-width/2, width/2 - 1, by = 1),
        mean = means,
        ci_low = means - ci,
        ci_high = means + ci
    )
}

Answer

computeCoverage <- function(l) {
    features <- l$features
    width <- unique(width(features))

    ## features x positions matrix
    mat <- matrix(
        as.numeric(unlist(l$coverage[features])),
        ncol = width, byrow = TRUE
    )

    ## per-position mean and 95% CI
    n <- nrow(mat)
    means <- colMeans(mat, na.rm = TRUE)
    ci <- 1.96 * apply(mat, 2, sd, na.rm = TRUE) / sqrt(n)

    ## return a clean data.frame
    data.frame(
        distance = seq(-width/2, width/2 - 1, by = 1),
        mean = means,
        ci_low = means - ci,
        ci_high = means + ci
    )
}

df <- computeCoverage(l)
head(df)

1.9 A first look at the result

We now have a tidy data.frame describing the aggregated coverage. Let’s plot it — this preview is exactly what tomorrow’s fourth function, plotCoverage(), will formalize.

Question

Plot mean (with its confidence interval ribbon) as a function of distance. Do you see an enrichment of Scc1 signal around the center of its binding sites?

Hint

Use ggplot2, mapping distance to x and mean to y. Add a geom_ribbon(aes(ymin = ci_low, ymax = ci_high)) for the confidence interval and a geom_line() for the mean.

library(ggplot2)
ggplot(df, aes(x = ______, y = ______)) +
    geom_ribbon(aes(ymin = ______, ymax = ______), alpha = 0.2) +
    geom_line()

Answer

library(ggplot2)
ggplot(df, aes(x = distance, y = mean)) +
    geom_ribbon(aes(ymin = ci_low, ymax = ci_high), fill = "steelblue", alpha = 0.2) +
    geom_line(color = "steelblue") +
    labs(x = "Distance to center (bp)", y = "Mean coverage") +
    theme_bw()

Recap

In a single script, you have prototyped three of the four functions of our package (importFiles(), filterGRanges(), computeCoverage()) and sketched the fourth (plotCoverage()). Tomorrow, we will turn this prototype into a proper Bioconductor package: create the package skeleton, move these functions into R/ files, and document them.

The reference implementation of all four functions is available in Share/functions/, should you want to compare with your own version.

1.1 Set up your working environment

1.2 Importing a coverage track in R

1.3 Importing genomic features in R

1.4 Resizing features to a fixed width

1.5 Wrap it up: your first function, importFiles()

1.6 Extracting the coverage over features

1.7 Filtering features: your second function, filterGRanges()

1.8 Aggregating the coverage: your third function, computeCoverage()

1.9 A first look at the result

1.5 Wrap it up: your first function, `importFiles()`

1.7 Filtering features: your second function, `filterGRanges()`

1.8 Aggregating the coverage: your third function, `computeCoverage()`