In the following exercises, we will use the GenomicRanges package to explore range operations. We will use the AnnotationHub package to load DNAse hypersensitivity data from the ENCODE project. In practice, the ENCODE project published datasets like these as bed files. AnnotationHub has packaged these into GRanges objects that we can load and use directly. However, if you have a bed file of your own (peak calls, enhancer regions, etc.), you can load them into GRanges objects using rtracklayer::import.
In this exercise, we will use DNAse hypersensitivity data to practice working with a GRanges object.
AnnotationHub package to find the goldenpath/hg19/encodeDCC/wgEncodeUwDnase/wgEncodeUwDnaseK562PkRep1.narrowPeak.gz GRanges object. Load that into R as the variable dnase.library(AnnotationHub)
ah = AnnotationHub()
query(ah, "goldenpath/hg19/encodeDCC/wgEncodeUwDnase/wgEncodeUwDnaseK562PkRep1.narrowPeak.gz")
# the thing above should have only one record, so we can
# just grab it
dnase = query(ah, "goldenpath/hg19/encodeDCC/wgEncodeUwDnase/wgEncodeUwDnaseK562PkRep1.narrowPeak.gz")[[1]]dnase
class(dnase)mcols(dnase)library(GenomicFeatures)
table(seqnames(dnase))seqlevelsStyle to adjust the sequence names.seqlevels(dnase)
seqlevelsStyle(dnase)
seqlevelsStyle(dnase) = 'ensembl'
seqlevelsStyle(dnase)
seqlevels(dnase)sum(width(dnase))
sum(seqlengths(dnase))
sum(width(dnase))/sum(seqlengths(dnase))In this exercise, we are going to load the second DNAse hypersensitivity replicate to investigate overlaps with the first replicate.
goldenpath/hg19/encodeDCC/wgEncodeUwDnase/wgEncodeUwDnaseK562PkRep2.narrowPeak.gz. Load that as dnase2.dnase and dnase2? Are there are similar number?dnase2?dnase2 cover?sum(width(dnase))/sum(seqlengths(dnase))dnase that overlap with dnase2.resize (and pay attention to the fix argument) and width to do this expansion to dnase and call the new GRanges object “dnase_wide”.In this exercise, we are going to look at the overlap of DNAse sites relative to genes. To get started, install and load the TxDb.Hsapiens.UCSC.hg19.knownGene txdb object.
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
kg = TxDb.Hsapiens.UCSC.hg19.knownGeneflank method for GRanges objects. How could you use that to get the “promoter” regions of the transcripts? Let’s assume that the promoter region is 2kb upstream of the gene.flank(gx,2000)?promoters).proms = promoters(kg)summary(countOverlaps(proms))reduce.# reduce takes all overlapping regions and collapses them
# into a single region that spans all of the overlapping regions
prom_regions = reduce(proms)
# now we can check for overlaps
summary(countOverlaps(prom_regions))prop_proms = sum(width(prom_regions))/sum(seqlengths(prom_regions))
prop_dnase = sum(width(dnase))/sum(seqlengths(prom_regions))
# Iff the dnase and promoter regions are
# not related, then we would expect this number
# of DNAse overlaps with promoters.
prop_proms * prop_dnase * length(dnase) We’ll be using data from histone modification ChIP-seq experiments in human cells to illustrate the concepts of genomic ranges and features. The data consists of genomic intervals representing regions of the genome where specific histone modifications are enriched. These intervals are typically identified using ChIP-seq, a technique that maps protein-DNA interactions across the genome.
The ChIP-seq data is stored in a BED file format, which is a tab-delimited text file format commonly used to represent genomic intervals. Each line in the BED file corresponds to a genomic interval and contains information about the chromosome, start and end positions, and strand orientation of the interval. Additional columns may include metadata such as the signal strength or significance of the interval.
The AnnotationHub package in Bioconductor provides access to a wide range of genomic datasets, including ChIP-seq data. We can use this package to retrieve the ChIP-seq data for histone modifications in human cells and convert it into a GenomicRanges object for further analysis.
https://www.encodeproject.org/chip-seq/histone/
Let’s start by loading the AnnotationHub package and retrieving the ChIP-seq data for histone modifications in human cells. You can read more about the AnnotationHub package and how to use it in the Bioconductor documentation.
library(AnnotationHub)
ah <- AnnotationHub()There are multiple ways to search the AnnotationHub database. We’ve done that for you and here are the GRanges objects for each of four histone marks, and one histone mark replicate.
h3k4me1 <- ah[['AH25832']]
h3k4me3 <- ah[['AH25833']]
h3k9ac <- ah[['AH25834']]
h3k27me3 <- ah[['AH25835']]
h3k4me3_2 <- ah[['AH27284']]Each of these variables now represents the peak calls after a chip-seq experiment pulling down the histone mark of interest. In the encode project these records were bed files. The bed files have been converted to GRanges objects to allow computation within R.
Let’s say that we don’t know the behavior of the histone methylation marks with respect to CpG islands. We could ask the question, “What is the overlap of the histone peaks with CpG islands?”
We might want to actually count the number of bases of overlap between the methyl mark and CpG islands.
But some methyl marks are known to have very broad signals, meaning that there is a higher chance of overlapping CpG islands just because there are more methylated bases. We can adjust for this by “normalizing” for all possible bases covered by either set of peaks, using union. We might think of this as a sort of “enrichment score” of one set in another set.
Let’s write a small function to calculate our little enrichment score.
And give it a try:
range_enrichment_score(h3k4me1, cpg)