Genomic regions resulting from next-generation sequencing experiments and bioinformatics pipelines are more meaningful when annotated to genomic features. A SNP occurring in an exon, or an enhancer, is likely of greater interest than one occurring in an inter-genic region. It may be of interest to find that a particular transcription factor overwhelmingly binds in promoters, while another binds mostly in 3’UTRs. Hyper-methylation at promoters containing a CpG island may indicate different regulatory regimes in one condition compared to another.
annotatr provides genomic annotations and a set of functions to read, intersect, summarize, and visualize genomic regions in the context of genomic annotations.
The release version of
annotatr is available via Bioconductor, and can be installed as follows:
if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager") BiocManager::install("annotatr")
The development version of
annotatr can be obtained via the GitHub repository or Bioconductor. It is easiest to install development versions with the
devtools package as follows:
Changelogs for development releases will be detailed on GitHub releases.
There are three types of annotations available to annotatr:
The CpG islands are the basis for all CpG annotations, and are given by the
AnnotationHub package for the given organism. CpG shores are defined as 2Kb upstream/downstream from the ends of the CpG islands, less the CpG islands. CpG shelves are defined as another 2Kb upstream/downstream of the farthest upstream/downstream limits of the CpG shores, less the CpG islands and CpG shores. The remaining genomic regions make up the inter-CGI annotation.
CpG annotations are available for hg19, hg38, mm9, mm10, rn4, rn5, rn6.
The genic annotations are determined by functions from
GenomicFeatures and data from the
org.*.eg.db packages. Genic annotations include 1-5Kb upstream of the TSS, the promoter (< 1Kb upstream of the TSS), 5'UTR, first exons, exons, introns, CDS, 3'UTR, and intergenic regions (the intergenic regions exclude the previous list of annotations). The schematic below illustrates the relationship between the different annotations as extracted from the
TxDb.* packages via
Also included in genic annotations are intronexon and exonintron boundaries. These annotations are 200bp up/down stream of any boundary between an exon and intron. Important to note, is that the boundaries are with respect to the strand of the gene.
Non-intergenic gene annotations include Entrez ID and gene symbol information where it exists. The
org.*.eg.db packages for the appropriate organisms are used to provide gene IDs and gene symbols.
The genic annotations have populated
symbol columns. Respectively they are, the knownGene transcript name, Entrez Gene ID, and gene symbol.
Genic annotations are available for all hg19, hg38, mm9, mm10, rn4, rn5, rn6, dm3, and dm6.
FANTOM5 permissive enhancers were determined from bi-directional CAGE transcription as in Andersson et al. (2014), and are downloaded and processed for hg19 and mm9 from the FANTOM5 resource. Using the
rtracklayer::liftOver() function, enhancers from hg19 are lifted to hg38, and mm9 to mm10.
The long non-coding RNA (lncRNA) annotations are from GENCODE for hg19, hg38, and mm10. The lncRNA transcripts are used, and we eventually plan to include the lncRNA introns/exons at a later date. The lncRNA annotations have populated
symbol columns. Respectively they are, the Ensembl transcript name, Entrez Gene ID, and gene symbol. As per the
transcript_type field in the GENCODE anntotations, the biotypes are given in the
Chromatin states determined by chromHMM (Ernst and Kellis (2012)) in hg19 are available for nine cell lines (Gm12878, H1hesc, Hepg2, Hmec, Hsmm, Huvec, K562, Nhek, and Nhlf) via the UCSC Genome Browser tracks. Annotations for all states can be built using a shortcut like
hg19_Gm12878-chromatin, or specific chromatin states can be accessed via codes like
AnnotationHub Bioconductor package is a client for the AnnotationHub web resource. From the package description:
The AnnotationHub web resource provides a central location where genomic files (e.g., VCF, bed, wig) and other resources from standard locations (e.g., UCSC, Ensembl) can be discovered. The resource includes metadata about each resource, e.g., a textual description, tags, and date of modification. The client creates and manages a local cache of files retrieved by the user, helping with quick and reproducible access.
build_ah_annots() function, users can turn any resource of class
GRanges into an annotation for use in
annotatr. As an example, we create annotations for H3K4me3 ChIP-seq peaks in Gm12878 and H1-hesc cells.
# Create a named vector for the AnnotationHub accession codes with desired names h3k4me3_codes = c('Gm12878' = 'AH23256') # Fetch ah_codes from AnnotationHub and create annotations annotatr understands build_ah_annots(genome = 'hg19', ah_codes = h3k4me3_codes, annotation_class = 'H3K4me3') # The annotations as they appear in annotatr_cache ah_names = c('hg19_H3K4me3_Gm12878') print(annotatr_cache$get('hg19_H3K4me3_Gm12878'))
Users may load their own annotations from BED files using the
read_annotations() function, which uses the
rtracklayer::import() function. The output is a
type. If a user wants to include
symbol in their custom annotations they can be included as extra columns on a BED6 input file.
## Use ENCODE ChIP-seq peaks for EZH2 in GM12878 ## These files contain chr, start, and end columns ezh2_file = system.file('extdata', 'Gm12878_Ezh2_peak_annotations.txt.gz', package = 'annotatr') ## Custom annotation objects are given names of the form genome_custom_name read_annotations(con = ezh2_file, genome = 'hg19', name = 'ezh2', format = 'bed') print(annotatr_cache$get('hg19_custom_ezh2'))
To see what is in the
annotatr_cache environment, do the following:
The following example is based on the results of testing for differential methylation of genomic regions between two conditions using methylSig. The file (
inst/extdata/IDH2mut_v_NBM_multi_data_chr9.txt.gz) contains chromosome locations, as well as categorical and numerical data columns, and provides a good example of the flexibility of
read_regions() uses the
rtracklayer::import() function to read in BED files and convert them to
GRanges objects. The
score columns in a normal BED file can be used for categorical and numeric data, respectively. Additionally, an arbitrary number of categorical and numeric data columns can be appended to a BED6 file. The
extraCols parameter is used for this purpose, and the
rename_score columns allow users to give more descriptive names to these columns.
# This file in inst/extdata represents regions tested for differential # methylation between two conditions. Additionally, there are columns # reporting the p-value on the test for differential meth., the # meth. difference between the two groups, and the group meth. rates. dm_file = system.file('extdata', 'IDH2mut_v_NBM_multi_data_chr9.txt.gz', package = 'annotatr') extraCols = c(diff_meth = 'numeric', mu0 = 'numeric', mu1 = 'numeric') dm_regions = read_regions(con = dm_file, genome = 'hg19', extraCols = extraCols, format = 'bed', rename_name = 'DM_status', rename_score = 'pval') # Use less regions to speed things up dm_regions = dm_regions[1:2000] print(dm_regions)
Users may select annotations a la carte via the accessors listed with
builtin_annotations(), shortcuts, or use custom annotations as described above. The
hg19_cpgs shortcut annotates regions to CpG islands, CpG shores, CpG shelves, and inter-CGI. The
hg19_basicgenes shortcut annotates regions to 1-5Kb, promoters, 5'UTRs, exons, introns, and 3'UTRs. Shortcuts for other
builtin_genomes() are accessed in a similar way.
annotate_regions() requires a
GRanges object (either the result of
read_regions() or an existing object), a
GRanges object of the
annotations, and a logical value indicating whether to
ignore.strand when calling
GenomicRanges::findOverlaps(). The positive integer
minoverlap is also passed to
GenomicRanges::findOverlaps() and specifies the minimum overlap required for a region to be assigned to an annotation.
Before annotating regions, they must be built with
build_annotations() which requires a character vector of desired annotation codes.
# Select annotations for intersection with regions # Note inclusion of custom annotation, and use of shortcuts annots = c('hg19_cpgs', 'hg19_basicgenes', 'hg19_genes_intergenic', 'hg19_genes_intronexonboundaries', 'hg19_custom_ezh2', 'hg19_H3K4me3_Gm12878') # Build the annotations (a single GRanges object) annotations = build_annotations(genome = 'hg19', annotations = annots) # Intersect the regions we read in with the annotations dm_annotated = annotate_regions( regions = dm_regions, annotations = annotations, ignore.strand = TRUE, quiet = FALSE) # A GRanges object is returned print(dm_annotated)
annotate_regions() function returns a
GRanges, but it may be more convenient to manipulate a coerced
data.frame. For example,
# Coerce to a data.frame df_dm_annotated = data.frame(dm_annotated) # See the GRanges column of dm_annotaed expanded print(head(df_dm_annotated)) # Subset based on a gene symbol, in this case NOTCH1 notch1_subset = subset(df_dm_annotated, annot.symbol == 'NOTCH1') print(head(notch1_subset))
Given a set of annotated regions, it is important to know how the annotations compare to those of a randomized set of regions. The
randomize_regions() function is a wrapper of
regioneR::randomizeRegions() from the
regioneR package that creates a set of random regions given a
GRanges object. After creating the random set, they must be annotated with
annotate_regions() for later use. Only
builtin_genomes() can be used in our wrapper function. Downstream functions that support using random region annotations are
It is important to note that if the regions to be randomized have a particular property, for example they are CpGs, the
randomize_regions() wrapper will not preserve that property! Instead, we recommend using
universe being the superset of the data regions you want to sample from.
# Randomize the input regions dm_random_regions = randomize_regions( regions = dm_regions, allow.overlaps = TRUE, per.chromosome = TRUE) # Annotate the random regions using the same annotations as above # These will be used in later functions dm_random_annotated = annotate_regions( regions = dm_random_regions, annotations = annotations, ignore.strand = TRUE, quiet = TRUE)
When there is no categorical or numerical information associated with the regions,
summarize_annotations() is the only possible summarization function to use. It gives the counts of regions in each annotation type (see example below). If there is categorical and/or numerical information, then
summarize_categorical() may be used. Using random region annotations is only available for
# Find the number of regions per annotation type dm_annsum = summarize_annotations( annotated_regions = dm_annotated, quiet = TRUE) print(dm_annsum) # Find the number of regions per annotation type # and the number of random regions per annotation type dm_annsum_rnd = summarize_annotations( annotated_regions = dm_annotated, annotated_random = dm_random_annotated, quiet = TRUE) print(dm_annsum_rnd) # Take the mean of the diff_meth column across all regions # occurring in an annotation. dm_numsum = summarize_numerical( annotated_regions = dm_annotated, by = c('annot.type', 'annot.id'), over = c('diff_meth'), quiet = TRUE) print(dm_numsum) # Count the occurrences of classifications in the DM_status # column across the annotation types. dm_catsum = summarize_categorical( annotated_regions = dm_annotated, by = c('annot.type', 'DM_status'), quiet = TRUE) print(dm_catsum)
The 5 plot functions described below are to be used on the object returned by
annotate_regions(). The plot functions return an object of type
ggplot that can be viewed (
ggsave), or modified with additional
# View the number of regions per annotation. This function # is useful when there is no classification or data # associated with the regions. annots_order = c( 'hg19_custom_ezh2', 'hg19_H3K4me3_Gm12878', 'hg19_genes_1to5kb', 'hg19_genes_promoters', 'hg19_genes_5UTRs', 'hg19_genes_exons', 'hg19_genes_intronexonboundaries', 'hg19_genes_introns', 'hg19_genes_3UTRs', 'hg19_genes_intergenic') dm_vs_kg_annotations = plot_annotation( annotated_regions = dm_annotated, annotation_order = annots_order, plot_title = '# of Sites Tested for DM annotated on chr9', x_label = 'knownGene Annotations', y_label = 'Count') print(dm_vs_kg_annotations)
plot_annotation() can also use the annotated random regions in the
annotated_random argument to plot the number of random regions per annotation type next to the number of input data regions.
# View the number of regions per annotation and include the annotation # of randomized regions annots_order = c( 'hg19_custom_ezh2', 'hg19_H3K4me3_Gm12878', 'hg19_genes_1to5kb', 'hg19_genes_promoters', 'hg19_genes_5UTRs', 'hg19_genes_exons', 'hg19_genes_intronexonboundaries', 'hg19_genes_introns', 'hg19_genes_3UTRs', 'hg19_genes_intergenic') dm_vs_kg_annotations_wrandom = plot_annotation( annotated_regions = dm_annotated, annotated_random = dm_random_annotated, annotation_order = annots_order, plot_title = 'Dist. of Sites Tested for DM (with rndm.)', x_label = 'Annotations', y_label = 'Count') print(dm_vs_kg_annotations_wrandom)
# View a heatmap of regions occurring in pairs of annotations annots_order = c( 'hg19_custom_ezh2', 'hg19_H3K4me3_Gm12878', 'hg19_genes_promoters', 'hg19_genes_5UTRs', 'hg19_genes_exons', 'hg19_genes_introns', 'hg19_genes_3UTRs', 'hg19_genes_intergenic') dm_vs_coannotations = plot_coannotations( annotated_regions = dm_annotated, annotation_order = annots_order, axes_label = 'Annotations', plot_title = 'Regions in Pairs of Annotations') print(dm_vs_coannotations)
With numerical data, the
plot_numerical() function plots a single variable (histogram) or two variables (scatterplot) at the region level, faceting over the categorical variable of choice. It is possible to include two categorical variables to facet over (see below). Note, when the plot is a histogram, the distribution over all regions is plotted within each facet.
dm_vs_regions_annot = plot_numerical( annotated_regions = dm_annotated, x = 'mu0', facet = 'annot.type', facet_order = c('hg19_genes_1to5kb','hg19_genes_promoters', 'hg19_genes_5UTRs','hg19_genes_3UTRs', 'hg19_custom_ezh2', 'hg19_genes_intergenic', 'hg19_cpg_islands'), bin_width = 5, plot_title = 'Group 0 Region Methylation In Genes', x_label = 'Group 0') print(dm_vs_regions_annot)
dm_vs_regions_annot2 = plot_numerical( annotated_regions = dm_annotated, x = 'diff_meth', facet = c('annot.type','DM_status'), facet_order = list(c('hg19_genes_promoters','hg19_genes_5UTRs','hg19_cpg_islands'), c('hyper','hypo','none')), bin_width = 5, plot_title = 'Group 0 Region Methylation In Genes', x_label = 'Methylation Difference') print(dm_vs_regions_annot2)
dm_vs_regions_name = plot_numerical( annotated_regions = dm_annotated, x = 'mu0', y = 'mu1', facet = 'annot.type', facet_order = c('hg19_genes_1to5kb','hg19_genes_promoters', 'hg19_genes_5UTRs','hg19_genes_3UTRs', 'hg19_custom_ezh2', 'hg19_genes_intergenic', 'hg19_cpg_islands', 'hg19_cpg_shores'), plot_title = 'Region Methylation: Group 0 vs Group 1', x_label = 'Group 0', y_label = 'Group 1') print(dm_vs_regions_name)
plot_numerical_coannotations() shows the distribution of numerical data for regions occurring in any two annotations, as well as in one or the other annotation. For example, the following example shows CpG methylation rates for CpGs occurring in just promoters, just CpG islands, and both promoters and CpG islands.
dm_vs_num_co = plot_numerical_coannotations( annotated_regions = dm_annotated, x = 'mu0', annot1 = 'hg19_cpg_islands', annot2 = 'hg19_genes_promoters', bin_width = 5, plot_title = 'Group 0 Perc. Meth. in CpG Islands and Promoters', x_label = 'Percent Methylation') print(dm_vs_num_co)
# View the counts of CpG annotations in data classes # The orders for the x-axis labels. This is also a subset # of the labels (hyper, hypo, none). x_order = c( 'hyper', 'hypo') # The orders for the fill labels. Can also use this # parameter to subset annotation types to fill. fill_order = c( 'hg19_cpg_islands', 'hg19_cpg_shores', 'hg19_cpg_shelves', 'hg19_cpg_inter') # Make a barplot of the data class where each bar # is composed of the counts of CpG annotations. dm_vs_cpg_cat1 = plot_categorical( annotated_regions = dm_annotated, x='DM_status', fill='annot.type', x_order = x_order, fill_order = fill_order, position='stack', plot_title = 'DM Status by CpG Annotation Counts', legend_title = 'Annotations', x_label = 'DM status', y_label = 'Count') print(dm_vs_cpg_cat1)
# Use the same order vectors as the previous code block, # but use proportional fill instead of counts. # Make a barplot of the data class where each bar # is composed of the *proportion* of CpG annotations. dm_vs_cpg_cat2 = plot_categorical( annotated_regions = dm_annotated, x='DM_status', fill='annot.type', x_order = x_order, fill_order = fill_order, position='fill', plot_title = 'DM Status by CpG Annotation Proportions', legend_title = 'Annotations', x_label = 'DM status', y_label = 'Proportion') print(dm_vs_cpg_cat2)
plot_annotation() one may add annotations for random regions to the
annotated_random parameter of
plot_categorical(). The result is a Random Regions bar representing the distribution of random regions for the categorical variable used for
fill. NOTE: Random regions can only be added when
fill = 'annot.type'.
# Add in the randomized annotations for "Random Regions" bar # Make a barplot of the data class where each bar # is composed of the *proportion* of CpG annotations, and # includes "All" regions tested for DM and "Random Regions" # regions consisting of randomized regions. dm_vs_cpg_cat_random = plot_categorical( annotated_regions = dm_annotated, annotated_random = dm_random_annotated, x='DM_status', fill='annot.type', x_order = x_order, fill_order = fill_order, position='fill', plot_title = 'DM Status by CpG Annotation Proportions', legend_title = 'Annotations', x_label = 'DM status', y_label = 'Proportion') print(dm_vs_cpg_cat_random)
# View the proportions of data classes in knownGene annotations # The orders for the x-axis labels. x_order = c( 'hg19_custom_ezh2', 'hg19_genes_1to5kb', 'hg19_genes_promoters', 'hg19_genes_5UTRs', 'hg19_genes_exons', 'hg19_genes_introns', 'hg19_genes_3UTRs', 'hg19_genes_intergenic') # The orders for the fill labels. fill_order = c( 'hyper', 'hypo', 'none') dm_vs_kg_cat = plot_categorical( annotated_regions = dm_annotated, x='annot.type', fill='DM_status', x_order = x_order, fill_order = fill_order, position='fill', legend_title = 'DM Status', x_label = 'knownGene Annotations', y_label = 'Proportion') print(dm_vs_kg_cat)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.