If you have a set of genomic ranges, the GenomicDistributions R package can help you visualize properties of your region set. GenomicDistributions produces these nine types of plot:
GenomicDistributions can work with any reference genome, as long as you have some annotation data for it (like chromosome sizes and locations of genes). To make things easier for the common use cases, I’ve included in the package basic metadata for the most commonly used features from the reference genomes I use most (hg19, hg38, and mm10). If you need to produce similar plots with different features, partitions, or reference assemblies, that’s also possible, and not much more difficult; GenomicDistributions is very modular and will work with other bioconductor packages to process that data, but it requires one or two additional steps to curate your reference data.
In this vignette, we’ll go through examples of each of the plots using my common built-in features and partitions. If you want more control, there’s another advanced vignette that will introduce you how to define your own features, partitions, and chromosome sizes for custom analysis.
Before we start, I want to explain the design philosophy for functions in this package. Many R plotting packages combine calculations and plotting into one function. This may seem convenient, but if you want to plot the calculation results in a different way or combine them with something else, you can’t because you only have access to the final plot. GenomicDistributions divides these tasks so you can use the intermediate data to design your own custom plot, or use the calculated results directly for other analysis.
In GenomicDistributions, each plot type has two functions: a calculate function and a plot function. The calculate functions take your GRanges object and return a table of summary results. You can use these summary statistics how you like – aggregate them across multiple region sets, insert them into other plots you have, and so forth; or, you can simply plug that result directly into the corresponding plot function, which returns a ggplot2 object. Separating the calculation and plotting functions like this gives you more control over your results.
GenomicDistributions like this:
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("GenomicDistributions")
Start by loading up the package and getting your query set of regions as a GenomicRanges object. I’ve included an example bed file to demonstrate how these plots look. You can load it up like this:
library("GenomicDistributions") queryFile = system.file("extdata", "vistaEnhancers.bed.gz", package="GenomicDistributions") query = rtracklayer::import(queryFile)
Chromosome distribution plots help you visualize how your regions are distributed across chromosomes. To produce these, you’ll need to specify the chromosome lengths for your reference assembly. There are a few ways to do this.
For the common reference assemblies that I use (hg19, hg38, mm9, and mm10), I’ve included the metadata in the package. If you’re working with one of these genomes, making a plot of the distribution across chromosomes takes just a couple of lines of code:
# First, calculate the distribution: x = calcChromBinsRef(query, "hg19") # Then, plot the result: plotChromBins(x)
What if we want to do the same thing but on 2 query sets at the same time? No problem:
# Let's fudge a second region set by shifting the first one over query2 = GenomicRanges::shift(query, 1e6) queryList = GRangesList(vistaEnhancers=query, shifted=query2) x2 = calcChromBinsRef(queryList, "hg19") plotChromBins(x2)
## Warning: Removed 1 rows containing missing values (geom_bar).
These functions just do a naive binning across the genome. If you want to tweak the way the bins are handled, or use a different reference assembly, that’s also possible and is only slightly more complicated. There are other functions you can use for that, which are outlined in another vignette.
Feature distance distribution plots will show you how your regions are distributed with respect to the nearest feature of interest. To illustrate, we’ll use Transcription Start Sites (TSS) as our example feature of interest (but really, you can use any region set).
For TSS plots, since this is such a common use case, we can use a handy built-in function that does everything for us. It’s just one line of code to check distances from query to your TSSs (for common genomes), and then a second line of code to plot those distances:
# Calculate the distances: TSSdist = calcFeatureDistRefTSS(query, "hg19") # Then plot the result: plotFeatureDist(TSSdist, featureName="TSS")
This plot uses log-scale increasing bins to show how your regions are distributed. Now, let’s make a similar plot with multiple region sets input:
TSSdist2 = calcFeatureDistRefTSS(queryList, "hg19") plotFeatureDist(TSSdist2, featureName="TSS")
You can also plot a tiled version that aligns them all vertically:
plotFeatureDist(TSSdist2, featureName="TSS", tile=TRUE)
If you want to check distances to other features, that’s no problem;
calcFeatureDistRefTSS() is really just a wrapper for the workhorse function,
calcFeatureDist(). To show how this works, get some features you want to check the distance to. Here, let’s just shift our query set by a normally distributed random number:
featureExample = GenomicRanges::shift(query, round(rnorm(length(query), 0,1000)))
Now, with these features, we just use the
calcFeatureDist function to calculate the distances. This function uses the fast rolling joins from
data.table under the hood, so it completes very quickly. The result of this gets piped right into the plotting function as before:
fdd = calcFeatureDist(query, featureExample) plotFeatureDist(fdd)
Genomic partition distribution plots show you how your regions are distributed across genome annotation classes. This is most commonly used to show the distribution over element types, such as promoters, exons, introns, or intergenic regions. GenomicDistributions provides 3 types of partition distribution plots: percentages, expected, and cumulative.
The most basic partition plot just provides a barplot of percentage region overlaps. You can produce one or two-set plots like so:
gp = calcPartitionsRef(query, "hg19") plotPartitions(gp)
gp2 = calcPartitionsRef(queryList, "hg19") plotPartitions(gp2)