```{r setup2, include=FALSE, cache=FALSE, eval=TRUE}
library(knitr)
opts_chunk$set(root.dir=getwd(),
fig.path='Figures',
dev='png',
fig.show='hold',
cache=TRUE)
library(genomation)
library(genomationData)
library(GenomicRanges)
set.seed=10
```
# `genomation`: a toolkit for annotation and visualization of genomic data
## Introduction
-----
`genomation` is a toolkit for annotation
and in bulk visualization of genomic features (scored or unscored)
over predefined regions.
The genomic features which the package can handle can
be anything with a minimal information of chromosome, start and end.
The features could have any length and most of the time they are
associated with a score. Typical examples of such data sets include aligned
reads from high-throughput sequencing (HTS) experiments, percent methylation
values for CpGs (or other cytosines), locations of transcription factor binding,
and so on. On the other hand, throughout the vignette we use the phrase
"genomic annotation" to refer to the regions of the genome associated with a
potential function which does not necessarily have a score
(examples: CpG islands, genes, enhancers, promoter, exons, etc. ).
These genomic annotations are usually the regions of interest, and distribution
of genomic features over/around the annotations are can make the way for
biological interpretation of the data.
The pipeline for computational knowledge extraction consists of three steps:
data filtering, integration of data from multiple sources or generation of
predictive models and biological interpretation of produced models, which leads
to novel hypotheses that can be tested in the wetlab.`genomation` aims
to facilitate the integration of multiple sources of genomic features with
genomic annotation or already published experimental results.
## Access the data
-----
High-throughput data which will be used to show the functonality of the
`genomation` are located in two places. The annotation and cap analysis
of gene expression (CAGE) data comes prepared with the `genomation` package, while
the raw HTS data can be found in the sister package `genomationData`.
To install the `genomation` and the complementary data package the from
the Bioconductor repository, copy and paste the following lines into your R interpreter:
```{r eval=FALSE}
biocLite('genomationData')
biocLite('genomation')
```
The r Biocpkg("genomationData") vignette contains a verbose description of
included files.
To list the available data, type:
```{r listgenomationData, eval=TRUE}
list.files(system.file('extdata',package='genomationData'))
```
To see the descriptions of the files:
```{r genomationDataInfo, eval=TRUE, tidy=TRUE}
sampleInfo = read.table(system.file('extdata/SamplesInfo.txt',
package='genomationData'),header=TRUE, sep='\t')
sampleInfo[1:5,1:5]
```
Basic annotation data and processed experimental data can be found within the
`genomation` package.
The data can be accesed throught the data command or located
in the extdata folder.
```{r genomationInfo, eval=FALSE, tidy=TRUE}
library(genomation)
data(cage)
data(cpgi)
list.files(system.file('extdata', package='genomation'))
```
## Data input
-----
One of larger hindrances in computational genomics stems from the myriad of
formats that are used to store the data. Although some formats have been
selected as de facto standards for specific kind of biological data (e.g. BAM, VCF),
almost all publications come with supplementary tables that do not have the
same structure, but hold similar information. The tables usually have a tabular
format, contain the locationof elements in genomic coordinates and various
metadata colums. `genomation` contais functions to read genomic
features and genomic annotation provided they are in a tabular format.
These functions will read the data from flat files into GRanges or GRangesList
objects.
readGeneric is the workhorse of the `genomation` package. It is a function
developed specifically for input of genomic data in tabular formats, and their
conversion to a GRanges object.
By default, the function persumes that the file is a standard .bed file
containing columns chr, start, end.
```{r readGeneric1, eval=TRUE, tidy=TRUE}
library(genomation)
tab.file1=system.file("extdata/tab1.bed", package = "genomation")
readGeneric(tab.file1, header=TRUE)
```
If the file contains meta data columns (as in extended bed format),
it is possible to read all or some of the additional columns.
To select columns which you want to read in, use the meta.col argument
```{r readGeneric2, eval=TRUE, tidy=TRUE}
readGeneric(tab.file1, header=TRUE, keep.all.metadata=TRUE)
readGeneric(tab.file1, header=TRUE, meta.col = list(CpGnum=4))
```
If the file contains header, the function will automatically recognize the
columns using the header names.
```{r readGeneric3, eval=TRUE, tidy=TRUE}
readGeneric(tab.file1, header=TRUE, keep.all.metadata=TRUE)
```
If the files have permutted columns, such that the first three do not represent
chromosome, start and end, you can select an arbitrary set of columns using the
chr, start and end arguments.
```{r readGeneric4, eval=TRUE, tidy=TRUE}
tab.file2=system.file("extdata/tab2.bed", package = "genomation")
readGeneric(tab.file2, chr=3, start=2, end=1)
```
readGeneric function can easily be extended to read almost any kind of
biological data. As an example we have provided convenience functions to read
the Encode narrowPeak and broadPeak formats, and gtf formatted files.
```{r readGff1, eval=TRUE, tidy=TRUE}
gff.file=system.file("extdata/chr21.refseq.hg19.gtf", package = "genomation")
gff = gffToGRanges(gff.file)
head(gff)
```
In order to split the last column of the gff file, use the
split.group=TRUE argument.
```{r readGff2, eval=FALSE, tidy=TRUE}
gff = gffToGRanges(gff.file, split.group=TRUE)
head(gff)
```
There are specific functions to read genomic annotation from flat bed files.
readFeatureFlank is a convenience function used to get the ranges flanking
the set of interest. As an example, it could be used to get the CpG island shores,
which have been shown to harbour condition specific differential methylation.
```{r readFeatureFlank, eval=TRUE, tidy=TRUE}
# reading genes stored as a BED files
cpgi.file = system.file("extdata/chr21.CpGi.hg19.bed", package = "genomation")
cpgi.flanks = readFeatureFlank(cpgi.file)
head(cpgi.flanks$flanks)
```
## Extraction and visualization of genomic data
-----
A standard step in a computational genomics experiment is visualization of
average enrichment over a certain predefined set of ranges, such as mean coverage
of a certain histone modification around a transcription factor binding site,
or visualization of histone positions around transcription start sites.
### Extraction of data over genomic winows
*ScoreMatrix* and *ScoreMatrixBin* are functions used to extract
data over predefined windows.
*ScoreMatrix* is used when all of the windows
have the same width, such as a designated area around the transcription start
site, while the *ScoreMatrixBin* is designed for use with windows of
unequal width (e.g. enrichment of methylation over exons).
Both functions have 2 main arguments: *target* and
*windows*. *target* is the data that we want to extract, while the
*windows* represents the regions over which we want to see the enrichment.
The target data can be in 3 forms: a GRanges, a RLeList or a path to an indexed
.bam file. The *windows* must be GRanges object.
As an example we will extract the density of cage tags around the promoters on
the human chromosome 21.
```{r ScoreMatrixExample, eval=TRUE, tidy=TRUE}
data(cage)
data(promoters)
sm = ScoreMatrix(target=cage, windows=promoters)
sm
```
*ScoreMatrixBin* function has an additional *bin.num* argument which
specifies the number of bins that will represent each windot window
(ie. it converts windows of unequal width into ones of equal width.).
By default, the binning function is set to *mean*.
```{r ScoreMatrixBinExample, eval=FALSE, tidy=TRUE}
data(cage)
gff.file = system.file('extdata/chr21.refseq.hg19.gtf', package="genomation")
exons = gffToGRanges(gff.file, filter='exon')
sm = ScoreMatrixBin(target=cage, windows=exons, bin.num=50)
sm
```
Running *ScoreMatrixBin* with *bin.num=50* on
a set of exons warned us that some of the exons are shorter than 50 bp and
were thus removed from the set before binning.
The rownames of the resulting *ScoreMatrix* object correspond
to the ranges that were used to construct the windows (e.g. row name 10 means
that the 10th element in the target GRanges object was used to extract the data).
If a certain rowname is not present in the ScoreMatrix object, that means that
the corresponding range was filtered out (e.g. the range could have been on a
chromosome that was not present in the target).
To simultaneously work on multiple files you can use the *ScoreMatrixList*
function. The function also has 2 obligatory arguments *targets* and
*windows*. While the *windows* is the same as in *ScoreMatrix*,
the targets argument contains results from multiple experiments.
It can be in one of the three formats: a list of RleLists, a list of GRanges
(or a GRangesList object), or a character vector designating a set of
.bam or .bigWig files.
```{r ScoreMatrixList, eval=TRUE, tidy=TRUE}
data(promoters)
data(cpgi)
data(cage)
cage$tpm = NULL
targets = list(cage=cage, cpgi=cpgi)
sm = ScoreMatrixList(targets=targets, windows=promoters, bin.num=50)
sm
```
If all of the windows have the same width *ScoreMatrixList* will use the
*ScoreMatrix* function. That can be overridden by explicitly specifying the
*bin.num* argument, as we did in the example.
### Visualization of multiple genomic experiments
There are 2 basic modes of visualization of enrichment over windows: either
as a heatmap, or as a histogram. *heatMatrix*, *plotMeta* and
*multiHeatMatrix* are functions for visualization of *ScoreMatrix*
and *ScoreMatrixList* objects.
We will plot the distribution of CAGE tags around promoters on human chr21.
```{r heatMatrix1, eval=FALSE, tidy=TRUE, fig.keep='all'}
data(cage)
data(promoters)
sm = ScoreMatrix(target=cage, windows=promoters)
heatMatrix(sm, xcoords=c(-1000, 1000))
plotMeta(sm, xcoords=c(-1000, 1000))
```
The *heatMatrix* function can also take a list of numeric vectors
designating row names, or a factor variable that represent our
annotation over the windows.
```{r heatMatrix2, eval=FALSE, tidy=TRUE, fig.width=4, fig.height=3}
library(GenomicRanges)
data(cage)
data(promoters)
data(cpgi)
sm = ScoreMatrix(target=cage, windows=promoters, strand.aware=TRUE)
cpg.ind = which(countOverlaps(promoters, cpgi)>0)
nocpg.ind = which(countOverlaps(promoters, cpgi)==0)
heatMatrix(sm, xcoords=c(-1000, 1000), group=list(CpGi=cpg.ind, noCpGi=nocpg.ind))
```
Because the enrichment in windows can have a high dynamic range, it is sometimes
convenient to scale the matrix before plotting.
```{r heatMatrixScales, eval=FALSE, tidy=TRUE, echo=TRUE, fig.width=4, fig.height=3, dpi=300}
sm.scaled = scaleScoreMatrix(sm)
heatMatrix(sm.scaled, xcoords=c(-1000, 1000), group=list(CpGi=cpg.ind, noCpGi=nocpg.ind))
```
Several experiments can be plotted in a side by side fashion using a combination
of ScoreMatrixList and multiHeatMatrix.
```{r multiHeatMatrix1, eval=FALSE, tidy=TRUE, fig.width=4.5, fig.height=3, dpi=300}
cage$tpm = NULL
targets = list(cage=cage, cpgi=cpgi)
sml = ScoreMatrixList(targets=targets, windows=promoters, bin.num=50, strand.aware=TRUE)
multiHeatMatrix(sml, xcoords=c(-1000, 1000))
```
We can put the clustering function clustfun to see whether there are any patterns present in the data.
```{r multiHeatMatrix2, eval=FALSE, tidy=TRUE, fig.width=4.5, fig.height=3, dpi=300}
multiHeatMatrix(sml, xcoords=c(-1000, 1000), clustfun=function(x) kmeans(x, centers=2)$cluster)
```
More advance usage of the ScoreMatrix family of functions and their
visualtization can be found in the specific use-cases at the end of the vignette.
## Annotation of genomic features
-----
Searching for correlation between sets of genomic features is a standard exploratory
method in computational genomics. It is usually done by looking at the overlap between
2 or more sets of ranges and calculating various overlap statistics.
`genomation` contains two sets of functions for annotation of ranges:
the first one is used to facilitate the general annotation of any sets of ranges,
while the second one is used to annotate a given feature with gene structures (promoter,
exon, intron).
### Annotation by generic features
Firstly, we will select the broadPeak files from the `genomatonData` package,
and read in the peaks for the Ctcf transcription factor
```{r getFiles_annotation, eval=TRUE, tidy=TRUE}
library(genomationData)
genomationDataPath = system.file('extdata',package='genomationData')
sampleInfo = read.table(file.path(genomationDataPath, 'SamplesInfo.txt'),
header=TRUE, sep='\t', stringsAsFactors=FALSE)
peak.files=list.files(genomationDataPath, full.names=TRUE, pattern='broadPeak')
names(peak.files) = sampleInfo$sampleName[match(basename(peak.files),sampleInfo$fileName)]
ctcf.peaks = readBroadPeak(peak.files['Ctcf'])
```
Now we will annotate the human the Ctcf binding sites using the CpG islands.
Because the CpG islands are restricted to chromosomes 21 and 22, we will set the
*intersect.chr = TRUE*, which will limit the analysis only to the
chromosomes that are present in both data sets.
```{r readAllPeaks_annotation, eval=TRUE, tidy=TRUE}
data(cpgi)
peak.annot = annotateWithFeature(ctcf.peaks, cpgi, intersect.chr=TRUE)
peak.annot
```
The output of the annotateWithFeature function shows three types of information:
The total number of elements in the target dataset, the percentage of target
dataset that overlaps with the feature dataset. And the percentage of the feature
elements that overlap the target.
### Annotation of genomic features by gene structures
To find the distribution of our designated features around gene structures,
we will first read the transcript features from a file using the
*readTranscriptFeatures* function. *readTranscriptFeatures* reads
a bed12 formatted file and parses the coordinates into a *GRangesList* containing
four elements: exons, introns. promoters and transcription start sites (TSSes).
```{r readTranscriptFeatures, eval=TRUE, tidy=TRUE}
bed.file=system.file("extdata/chr21.refseq.hg19.bed", package = "genomation")
gene.parts = readTranscriptFeatures(bed.file)
```
*annotateWithGeneParts* will give us the overlap statistics between our
CTCF peaks and gene structures. We will again use the *intersect.chr=TRUE*
to limit the analysis.
```{r annotateCtcf, eval=TRUE, tidy=TRUE}
ctcf.annot=annotateWithGeneParts(ctcf.peaks, gene.parts, intersect.chr=TRUE)
ctcf.annot
```
annotateWithGeneParts can also take a set of feature ranges as an argument.
We will use the readGeneric function to load all of the broadPeak files in the
`genomationData`, which we will then annotate.
```{r annotateTargetList, eval=TRUE, tidy=TRUE}
peaks = GRangesList(lapply(peak.files, readGeneric))
names(peaks) = names(peak.files)
annot.list = annotateWithGeneParts(peaks, gene.parts, intersect.chr=TRUE)
```
Gene annotation of multiple feature objects can be visualized in a form of a heatmap,
where rows represent samples, columns the gene structure, and the value is the
percentage of overlap given by priority. If *cluster=TRUE*, then the function
will use hierarhical clustering to order the heatmap.
```{r plotGeneAnnotation, eval=FALSE, tidy=TRUE, dpi=300}
plotGeneAnnotation(annot.list)
plotGeneAnnotation(annot.list, cluster=TRUE)
```
## Use cases for `genomation` package
-----
The `genomation` package provides generalizable functions for genomic data analysis
and visualization. Below we will demonstrate the functionality on specific use cases
### Visualization of ChiP sequencing data
We will visualize the binding profiles of 6 transcription factors around the
Ctcf binding sites.
In the fist step we will select the *.bam files containing mapped reads.
```{r selectBamChipseq, eval=TRUE, tidy=TRUE}
genomationDataPath = system.file('extdata',package='genomationData')
bam.files = list.files(genomationDataPath, full.names=TRUE, pattern='bam$')
bam.files = bam.files[!grepl('Cage', bam.files)]
```
Firstly, we will read in the Ctcf peaks, filter regions from human chromosome 21,
and order them by their signal values. In the end we will resize all ranges to
have a uniform width of 500 bases, fixed on the center of the peak.
```{r readCtcfPeaks, eval=TRUE, tidy=TRUE}
ctcf.peaks = readBroadPeak(file.path(genomationDataPath,
'wgEncodeBroadHistoneH1hescCtcfStdPk.broadPeak.gz'))
ctcf.peaks = ctcf.peaks[seqnames(ctcf.peaks) == 'chr21']
ctcf.peaks = ctcf.peaks[order(-ctcf.peaks$signalValue)]
ctcf.peaks = resize(ctcf.peaks, width=1000, fix='center')
```
In order to extract the coverage values of all transcription factors around
chipseq peaks, we will use the *ScoreMatrixList* function.
*ScoreMatrixList* assign names to each element of the list based on the
names of the bam files. We will use the names of the files to find the
corresponding names of each sample in the SamplesInfo.txt.
Using the *heatmapProfile* on our *ScoreMatrixList*,
we can plot the underlying signal side by side.
```{r ctcfScoreMatrixList, eval=FALSE, tidy=TRUE, fig.width=5, fig.height=3.5, dpi=300}
sml = ScoreMatrixList(bam.files, ctcf.peaks, bin.num=50, type='bam')
sampleInfo = read.table(system.file('extdata/SamplesInfo.txt',
package='genomationData'),header=TRUE, sep='\t')
names(sml) = sampleInfo$sampleName[match(names(sml),sampleInfo$fileName)]
multiHeatMatrix(sml, xcoords=c(-500, 500), col=c('lightgray', 'blue'))
```
Heatmap profile of unscaled coverage shows a slight colocalization of
Ctcf, Rad21 and Znf143.
Because of the large range of signal values in chipseq peaks, the *heatmapProfile*
will not show the true extent of colocalization. To get around this, it is advisable
to independently scale the rows of each element in the ScoreMatrixList.
```{r plotScaledProfile, eval=FALSE, tidy=TRUE, fig.width=5, fig.height=3.5, dpi=300}
sml.scaled = scaleScoreMatrixList(sml)
multiHeatMatrix(sml.scaled, xcoords=c(-500, 500), col=c('lightgray', 'blue'))
```
Heatmap profile of scaled coverage shows much stronger colocalization
of the transcription factors; nevertheless, it is evident that some of the CTCF
peaks have a very weak enrichment.
### Combinatorial binding of transcription factors
In the first step we will read all peak files into a GRanges list. We will use
the SamplesInfo file from the `genomationData` to get he names of the
samples. Four of the peak files are in the Encode broadPeak format, while one is
in the narrowPeak. To read the files, we will use the readGeneric
function. It enables us to select from the files only columns of interest.
As the last step, we will restrict ourselves to peaks that are located on
chromosome 21 and have width 100 and 1000 bp
```{r eadAllPeaks, eval=TRUE, tidy=TRUE}
genomationDataPath = system.file('extdata',package='genomationData')
sampleInfo = read.table(file.path(genomationDataPath, 'SamplesInfo.txt'),
header=TRUE, sep='\t', stringsAsFactors=FALSE)
peak.files = list.files(genomationDataPath, full.names=TRUE, pattern='Peak.gz$')
peaks = list()
for(i in 1:length(peak.files)){
file = peak.files[i]
name=sampleInfo$sampleName[match(basename(file),sampleInfo$fileName)]
message(name)
peaks[[name]] = readGeneric(file, meta.col=list(score=5))
}
peaks = GRangesList(peaks)
peaks = peaks[seqnames(peaks) == 'chr21' & width(peaks) < 1000 & width(peaks) > 100]
```
To find the combination of binding sites we will use the *findFeatureComb*
function. It takes a granges list object, finds the union of the
ranges and designates each range by the combination of overlaps from the original set.
By default, the returned ranges will have a numeric *class*
meta data column, which designates the correponding combination.
If you are interested in the names of the TF which make
the combinations, put the *use.names=TRUE*.
```{r findFeatureComb, eval=TRUE, tidy=TRUE}
tf.comb = findFeatureComb(peaks, width=1000)
```
To visualize the results, we will plot the enricment of resulting regions.
Before doing that we will order the regions by their *class* argument.
```{r visualizeFeatureComb, eval=FALSE, tidy=TRUE, fig.width=5, fig.height=3, dpi=300}
tf.comb = tf.comb[order(tf.comb$class)]
bam.files = list.files(genomationDataPath, full.names=TRUE, pattern='bam$')
bam.files = bam.files[!grepl('Cage', bam.files)]
sml = ScoreMatrixList(bam.files, tf.comb, bin.num=20, type='bam')
names(sml) = sampleInfo$sampleName[match(names(sml),sampleInfo$fileName)]
sml.scaled = scaleScoreMatrixList(sml)
multiHeatMatrix(sml.scaled, xcoords=c(-500, 500), col=c('lightgray', 'blue'))
```
The plot shows perfectly how misleading the peak calling process can be.
Although the plots show that CTFC, Rad21 and Znf143 have almost perfect
colozalization, peak callers have trouble identifying peaks in regions with
lower enrichments and as a result, we get different statistics
when using overlaps.
### Using data from AnnotationHub
We can also use data from AnnotationHub since it can return data as GRanges object. Below we download CpG island and DNAse peak locations from AnnotationHub and get a scoreMatrix on the CpG islands.
```{r annotationHubExample1, eval=FALSE, tidy=TRUE}
library(AnnotationHub)
ah = AnnotationHub()
# retrieve CpG island data from annotationHub
cpgi_query <- query(ah, c("CpG Islands", "UCSC", "hg19"))
cpgi <- ah[[names(cpgi_query)]]
dnase_query <- query(ah, c("wgEncodeOpenChromDnase8988tPk.narrowPeak.gz", "UCSC", "hg19"))
dnaseP <- ah[[names(dnase_query)]]
sm = ScoreMatrixBin(target=dnaseP, windows=cpgi, strand.aware=FALSE)
```
