annotatr: Making sense of genomic regionsAnnotationHub AnnotationsGenomic 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 is available from Bioconductor devel, which needs the devel version of R and Bioconductor:
BiocManager::install("annotatr", version = "devel")
The source is on GitHub. Changes in each version are listed in the package’s NEWS, news(package = "annotatr").
There are three kinds of annotations available to annotatr:
AnnotationHub: any GRanges resource in the Bioconductor AnnotationHub web resource.TxDb or EnsDb of gene models.Built-in annotations are named by codes of the form [genome]_[group]_[type], e.g. hg38_genes_promoters or mm10_cpg_islands, and all codes are listed by builtin_annotations(). builtin_annotations_table() shows which groups are available for which genome builds:
| genome | genes | mane | canonical | cpgs | enhancers | chromatin | lncrna | ccres |
|---|---|---|---|---|---|---|---|---|
| dm3 | ✓ | |||||||
| dm6 | ✓ | ✓ | ||||||
| danRer10 | ✓ | ✓ | ||||||
| danRer11 | ✓ | ✓ | ✓ | |||||
| galGal5 | ✓ | ✓ | ||||||
| hg19 | ✓ | ✓ | ✓ | ✓ | ✓ | |||
| hg38 | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | |
| mm9 | ✓ | ✓ | ✓ | |||||
| mm10 | ✓ | ✓ | ✓ | ✓ | ✓ | |||
| mm39 | ✓ | ✓ | ✓ | |||||
| oviariramb2 | ✓ | ✓ | ✓ | |||||
| rn4 | ✓ | ✓ | ||||||
| rn5 | ✓ | ✓ | ||||||
| rn6 | ✓ | ✓ | ||||||
| rn7 | ✓ | ✓ | ✓ |
Shortcuts build a whole group: [genome]_basicgenes, hg38_basicmane, [genome]_basiccanonical, [genome]_cpgs, [genome]_ccres, and hg19_[cell line]-chromatin. The gene shortcuts build 1-5Kb upstream of the TSS, promoters, 5’UTRs, exons, introns, and 3’UTRs.
The CpG islands are the basis for all CpG annotations. They are the UCSC Genome Browser’s CpG island track for the genome, from AnnotationHub for hg19, mm9, rn4, and rn5, from the UCSC database for most other genomes, and from the UCSC GenArk assembly hub for sheep (oviariramb2). 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.
The genic annotations are determined by functions from GenomicFeatures and data from the TxDb.* and 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 GenomicFeatures functions.
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.
For sheep (oviariramb2), which has no TxDb.* or org.*.eg.db packages, gene annotations are built from the Ensembl EnsDb in AnnotationHub. Its CpG islands come from the UCSC GenArk assembly hub.
Non-intergenic gene annotations have populated tx_id, gene_id, symbol, entrez_id, and ensembl_id columns. The tx_id and gene_id are the transcript and gene IDs of the source of the gene models: Entrez Gene IDs for most genomes, Ensembl IDs for rn4 and sheep, and FlyBase IDs for dm3 and dm6. Whatever the source’s IDs are, the symbol, entrez_id, and ensembl_id columns give the gene symbol, Entrez Gene ID, and Ensembl gene ID (without version), from the org.*.eg.db package for the organism, or the EnsDb for sheep. When a gene ID maps to more than one symbol or ID, the first is used.
Genic annotations are available for every genome in builtin_genomes().
The genic annotations above include every transcript of a gene, so a region in an exon of one isoform may be in an intron of another. To annotate regions to one main isoform per gene, use the MANE Select transcripts for hg38 (Morales et al., 2022). MANE (Matched Annotation from NCBI and EMBL-EBI) chooses one representative transcript for each human protein-coding gene, which is identical in RefSeq and Ensembl/GENCODE.
The MANE annotations have the same types as the genic annotations, except intergenic: MANE only covers protein-coding genes, so the regions outside MANE transcripts would include non-coding genes. The codes are of the form hg38_mane_promoters, and the hg38_basicmane shortcut builds the same types as hg38_basicgenes. The tx_id column is the Ensembl transcript ID, and gene_id is the Entrez Gene ID, as for hg38_genes_*. The annotations are built from the MANE GTF from NCBI, which needs the txdbmaker package, and are cached like the others. MANE Plus Clinical transcripts, which some clinical genes have in addition to MANE Select, are not included.
# One transcript per gene
mane_annotations = build_annotations(genome = 'hg38', annotations = 'hg38_basicmane')
# Or alongside all transcripts, e.g. to compare them
annotations = build_annotations(genome = 'hg38',
annotations = c('hg38_genes_promoters', 'hg38_mane_promoters'))
Ensembl chooses one canonical transcript for every gene, of every biotype (protein-coding, lncRNA, etc.). The [genome]_canonical_* annotations are the genic annotations built from only those transcripts, for the current assemblies hg38, mm39, rn7, danRer11, dm6, and oviariramb2, from the Ensembl 113 EnsDb in AnnotationHub. The [genome]_basiccanonical shortcut builds the same types as basicgenes, and intergenic regions are available too, since every gene has a canonical transcript.
How Ensembl chooses the canonical transcript depends on the species. For human, it is the MANE Select transcript when there is one (see above), and otherwise is chosen from conservation, expression, APPRIS and UniProt annotations, CDS length, and clinical variants. For other species, Ensembl chooses by biotype (protein-coding first) and then the longest combined exon length, so the canonical transcript is roughly the longest protein-coding transcript.
The gene_id and ensembl_id columns are Ensembl gene IDs, and symbol and entrez_id come from the EnsDb. Genes on alternate haplotypes and fix patches (e.g. chr6_GL000251v2_alt), which Ensembl annotates as separate genes, are left out, so there is one canonical transcript per gene on the reference assembly.
# One transcript per gene, for all genes, in mouse
canonical_annotations = build_annotations(genome = 'mm39',
annotations = c('mm39_basiccanonical', 'mm39_canonical_intergenic'))
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. For hg38 and mm10, the enhancer-like ENCODE cCREs (below) are a more recent alternative, identified directly in those genome builds.
The ENCODE registry of candidate cis-regulatory elements (cCREs), from SCREEN (ENCODE Project Consortium et al., 2020), is available for hg38 and mm10, version 4. cCREs are 150–350 bp regions of accessible chromatin, classified by their histone marks, CTCF, and transcription factor binding:
| Code | Class |
|---|---|
[genome]_ccre_PLS |
promoter-like |
[genome]_ccre_pELS |
proximal enhancer-like (within 2 kb of a TSS) |
[genome]_ccre_dELS |
distal enhancer-like |
[genome]_ccre_CA-H3K4me3 |
chromatin accessible, with H3K4me3 |
[genome]_ccre_CA-CTCF |
chromatin accessible, with CTCF |
[genome]_ccre_CA-TF |
chromatin accessible, with TF binding |
[genome]_ccre_CA |
chromatin accessible only |
[genome]_ccre_TF |
TF binding only |
The [genome]_ccres shortcut builds all 8 classes. The id column is ENCODE’s ID for each cCRE, its accession (e.g. EH38E2776516), which can be looked up in SCREEN. There are about 2.3 million cCREs for hg38 and 0.9 million for mm10, most of them distal enhancer-like, so annotations with all classes are large, and plots of all 8 classes can be busy.
ccre_annotations = build_annotations(genome = 'hg38', annotations = 'hg38_ccres')
The long non-coding RNA (lncRNA) transcripts are from GENCODE for hg19 (release 19), hg38 (release 31), and mm10 (release M6). The tx_id column is the Ensembl transcript ID, ensembl_id is the Ensembl gene ID, symbol is the gene symbol, and gene_id and entrez_id are the Entrez Gene ID matched by symbol. The transcript biotype, from GENCODE’s transcript_type field, is at the start of the id column.
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 hg19_chromatin_Gm12878-StrongEnhancer or hg19_chromatin_Gm12878-Repressed.
AnnotationHub AnnotationsThe 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.
Using the 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'))
## GRanges object with 57476 ranges and 7 metadata columns:
## seqnames ranges strand | id tx_id
## <Rle> <IRanges> <Rle> | <character> <logical>
## [1] chr1 713208-713477 * | H3K4me3_Gm12878:1 <NA>
## [2] chr1 713874-714056 * | H3K4me3_Gm12878:2 <NA>
## [3] chr1 714474-714750 * | H3K4me3_Gm12878:3 <NA>
## [4] chr1 715069-715388 * | H3K4me3_Gm12878:4 <NA>
## [5] chr1 724097-724311 * | H3K4me3_Gm12878:5 <NA>
## ... ... ... ... . ... ...
## [57472] chrX 154996923-154997189 * | H3K4me3_Gm12878:57472 <NA>
## [57473] chrX 154997422-154997785 * | H3K4me3_Gm12878:57473 <NA>
## [57474] chrX 155100454-155128015 * | H3K4me3_Gm12878:57474 <NA>
## [57475] chrX 155148379-155155444 * | H3K4me3_Gm12878:57475 <NA>
## [57476] chrX 155227027-155228269 * | H3K4me3_Gm12878:57476 <NA>
## gene_id symbol entrez_id ensembl_id type
## <logical> <logical> <character> <character> <character>
## [1] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [2] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [3] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [4] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [5] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## ... ... ... ... ... ...
## [57472] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [57473] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [57474] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [57475] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## [57476] <NA> <NA> <NA> <NA> hg19_H3K4me3_Gm12878
## -------
## seqinfo: 298 sequences (2 circular) from hg19 genome
Users may load their own annotations from BED files using the read_annotations() function, which uses the rtracklayer::import() function. The output is a GRanges with the same mcols() as built-in annotations: id, tx_id, gene_id, symbol, entrez_id, ensembl_id, and type. Any of tx_id, gene_id, symbol, entrez_id, and ensembl_id can be included as extra columns of a BED6 input file, with the extraCols argument, and the others are NA.
## 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')
## read_annotations() returns the custom annotations
ezh2_annotations = read_annotations(con = ezh2_file, genome = 'hg19', name = 'ezh2', format = 'bed')
print(ezh2_annotations)
## GRanges object with 2472 ranges and 7 metadata columns:
## seqnames ranges strand | id tx_id
## <Rle> <IRanges> <Rle> | <character> <character>
## [1] chr1 860063-860382 * | ezh2:1 <NA>
## [2] chr1 934911-935230 * | ezh2:2 <NA>
## [3] chr1 3573321-3573640 * | ezh2:3 <NA>
## [4] chr1 6301401-6301720 * | ezh2:4 <NA>
## [5] chr1 6301996-6302315 * | ezh2:5 <NA>
## ... ... ... ... . ... ...
## [2468] chrX 99880950-99881269 * | ezh2:2468 <NA>
## [2469] chrX 108514101-108514420 * | ezh2:2469 <NA>
## [2470] chrX 111981673-111981992 * | ezh2:2470 <NA>
## [2471] chrX 118109216-118109535 * | ezh2:2471 <NA>
## [2472] chrX 136114771-136115090 * | ezh2:2472 <NA>
## gene_id symbol entrez_id ensembl_id type
## <character> <character> <character> <character> <character>
## [1] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [2] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [3] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [4] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [5] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## ... ... ... ... ... ...
## [2468] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [2469] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [2470] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [2471] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## [2472] <NA> <NA> <NA> <NA> hg19_custom_ezh2
## -------
## seqinfo: 298 sequences (2 circular) from hg19 genome
Custom annotations can be combined with other annotations using c(), e.g. c(ezh2_annotations, build_annotations(genome = 'hg19', annotations = 'hg19_cpgs')). read_annotations() also stores them in the annotatr_cache environment for the current R session, with names of the form genome_custom_name, so they can be included by name in build_annotations(), as in the example below. To see what is in the annotatr_cache environment, do the following:
print(annotatr_cache$list_env())
## [1] "hg19_H3K4me3_Gm12878" "hg19_custom_ezh2"
The built-in genic annotations come from the TxDb.* packages, MANE, and Ensembl for the genomes annotatr knows about. To build the same annotation types from other gene models, e.g. GENCODE or RefSeq, or for a genome without built-in annotations, use build_txdb_annotations() with any TxDb or EnsDb object. A GTF or GFF3 file can be made into a TxDb with txdbmaker::makeTxDbFromGFF().
The genome and group arguments are required. The annotation codes are of the form [genome]_[group]_[type], e.g. mm39_gencode_promoters, so group names where the gene models came from, and annotations from different gene models can be used side by side. The annotations argument takes the gene annotation types, or the basicgenes shortcut (the default). Gene symbols come from an EnsDb, or from an OrgDb given with orgdb, with keytype saying what kind of gene IDs the TxDb has.
# A TxDb from the GENCODE mouse GTF, with Ensembl gene IDs
gencode_txdb = txdbmaker::makeTxDbFromGFF('gencode.vM38.annotation.gtf.gz')
gencode_annots = build_txdb_annotations(gencode_txdb, genome = 'mm39', group = 'gencode',
annotations = c('basicgenes', 'intergenic'),
orgdb = org.Mm.eg.db::org.Mm.eg.db, keytype = 'ENSEMBL')
# Combine with built-in annotations, and use as any other annotations
annotations = c(gencode_annots, build_annotations(genome = 'mm39', annotations = 'mm39_cpgs'))
These annotations are not cached on disk, since they are built from a TxDb already in memory. summarize_genes() and the plotting functions work with them as with built-in gene annotations.
Building annotations can take a while, especially gene annotations, which are built from the TxDb.* packages or an EnsDb. So build_annotations() saves each annotation it builds, and each file it downloads, in a cache on disk, and later calls load them from the cache in seconds. The cache persists between R sessions. A cached annotation is rebuilt automatically after annotatr is updated, and for gene annotations, after the TxDb.* or org.*.eg.db package is updated.
# See what is cached, with sizes and the versions each annotation was built with
list_cached_annotations()
# Build without reading or writing the cache
annotations = build_annotations(genome = 'hg19', annotations = 'hg19_cpgs', cache = FALSE)
# Remove the cached hg19 annotations and downloads, so they are built again
clear_cached_annotations(genome = 'hg19')
# Keep the cache somewhere else, e.g. on a cluster with a small home directory
options(annotatr.cache = '/path/with/space')
If annotations seem out of date or a cached file causes errors, see the troubleshooting section of ?"cached-annotations". This cache is separate from annotatr_cache, which holds custom annotations for the current R session only.
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 annotatr.
read_regions() uses the rtracklayer::import() function to read in BED files and convert them to GRanges objects. BED files are 0-based and GRanges are 1-based, so the BED line chr1 99 100 becomes the single base chr1:100. A single CpG written with a 1-based start, chr1 100 100, would become an empty range that overlaps no annotations, so check the coordinates of files from other tools. The name and 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_name and 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)
## GRanges object with 2000 ranges and 5 metadata columns:
## seqnames ranges strand | DM_status pval diff_meth
## <Rle> <IRanges> <Rle> | <character> <numeric> <numeric>
## [1] chr9 10850-10948 * | none 0.5045502 -10.7329047
## [2] chr9 10950-11048 * | none 0.2227126 8.7195270
## [3] chr9 28950-29048 * | none 0.5530958 0.0700847
## [4] chr9 72850-72948 * | hyper 0.0116294 44.8753244
## [5] chr9 72950-73048 * | none 0.1752872 17.7606626
## ... ... ... ... . ... ... ...
## [1996] chr9 35605150-35605248 * | none 0.274255 -0.0539158
## [1997] chr9 35605250-35605348 * | none 0.918064 0.0329283
## [1998] chr9 35605350-35605448 * | none 0.614312 -0.0977500
## [1999] chr9 35605450-35605548 * | none 1.000000 0.0000000
## [2000] chr9 35605550-35605648 * | none 0.814567 0.0349967
## mu0 mu1
## <numeric> <numeric>
## [1] 79.981920 90.7148252
## [2] 86.704015 77.9844878
## [3] 0.124081 0.0539963
## [4] 72.455413 27.5800883
## [5] 28.440368 10.6797057
## ... ... ...
## [1996] 0.000000 0.0539158
## [1997] 0.328024 0.2950959
## [1998] 0.130184 0.2279345
## [1999] 0.000000 0.0000000
## [2000] 0.118272 0.0832756
## -------
## seqinfo: 298 sequences (2 circular) from hg19 genome
Every GRanges has a Seqinfo describing its genome: the chromosome names, their lengths, whether they are circular, and the genome (e.g. hg19). Setting genome in read_regions() fills this in, and annotations from build_annotations() have it for their genome.
When regions are annotated, annotate_regions() checks the regions and annotations match:
hg38 and annotations from hg19, it gives an error, rather than annotating regions to the wrong places.2 in the regions and chr2 in the annotations, it gives an error. Annotations from build_annotations() use UCSC-style names (chr1, chr2, etc.), and GenomeInfoDb::seqlevelsStyle(regions) = 'UCSC' changes the regions to match.genome in read_regions().The Seqinfo of regions or annotations can be seen with Seqinfo::seqinfo(), and the genome with Seqinfo::genome().
# The genome of the regions read in above
print(unique(Seqinfo::genome(dm_regions)))
## [1] "hg19"
Users may select annotations a la carte via the codes 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. The other shortcuts, and the genomes they are available for, are described in the Annotations section above.
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)
## GRanges object with 24517 ranges and 6 metadata columns:
## seqnames ranges strand | DM_status pval diff_meth
## <Rle> <IRanges> <Rle> | <character> <numeric> <numeric>
## [1] chr9 10850-10948 * | none 0.50455 -10.7329
## [2] chr9 10850-10948 * | none 0.50455 -10.7329
## [3] chr9 10850-10948 * | none 0.50455 -10.7329
## [4] chr9 10850-10948 * | none 0.50455 -10.7329
## [5] chr9 10850-10948 * | none 0.50455 -10.7329
## ... ... ... ... . ... ... ...
## [24513] chr9 35605550-35605648 * | none 0.814567 0.0349967
## [24514] chr9 35605550-35605648 * | none 0.814567 0.0349967
## [24515] chr9 35605550-35605648 * | none 0.814567 0.0349967
## [24516] chr9 35605550-35605648 * | none 0.814567 0.0349967
## [24517] chr9 35605550-35605648 * | none 0.814567 0.0349967
## mu0 mu1 annot
## <numeric> <numeric> <GRanges>
## [1] 79.9819 90.7148 chr9:6917-10916:+
## [2] 79.9819 90.7148 chr9:6935-10934:+
## [3] 79.9819 90.7148 chr9:6937-10936:+
## [4] 79.9819 90.7148 chr9:6937-10936:+
## [5] 79.9819 90.7148 chr9:6937-10936:+
## ... ... ... ...
## [24513] 0.118272 0.0832756 chr9:35605185-35606184:-
## [24514] 0.118272 0.0832756 chr9:35605259-35605616:+
## [24515] 0.118272 0.0832756 chr9:35605259-35605835:+
## [24516] 0.118272 0.0832756 chr9:35605488-35605835:+
## [24517] 0.118272 0.0832756 chr9:35603969-35605991:*
## -------
## seqinfo: 298 sequences (2 circular) from hg19 genome
The 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_annotated expanded
print(head(df_dm_annotated))
## seqnames start end width strand DM_status pval diff_meth mu0
## 1 chr9 10850 10948 99 * none 0.5045502 -10.7329 79.98192
## 2 chr9 10850 10948 99 * none 0.5045502 -10.7329 79.98192
## 3 chr9 10850 10948 99 * none 0.5045502 -10.7329 79.98192
## 4 chr9 10850 10948 99 * none 0.5045502 -10.7329 79.98192
## 5 chr9 10850 10948 99 * none 0.5045502 -10.7329 79.98192
## 6 chr9 10850 10948 99 * none 0.5045502 -10.7329 79.98192
## mu1 annot.seqnames annot.start annot.end annot.width annot.strand
## 1 90.71483 chr9 6917 10916 4000 +
## 2 90.71483 chr9 6935 10934 4000 +
## 3 90.71483 chr9 6937 10936 4000 +
## 4 90.71483 chr9 6937 10936 4000 +
## 5 90.71483 chr9 6937 10936 4000 +
## 6 90.71483 chr9 6941 10940 4000 +
## annot.id annot.tx_id annot.gene_id annot.symbol annot.entrez_id
## 1 1to5kb:179973 ENST00000806493.1_1 100287596 DDX11L5 100287596
## 2 1to5kb:179974 ENST00000806495.1_1 100287596 DDX11L5 100287596
## 3 1to5kb:179975 ENST00000806497.1_1 100287596 DDX11L5 100287596
## 4 1to5kb:179976 ENST00000806496.1_1 <NA> <NA> <NA>
## 5 1to5kb:179977 ENST00000806494.1_1 100287596 DDX11L5 100287596
## 6 1to5kb:179978 ENST00000806498.1_1 100287596 DDX11L5 100287596
## annot.ensembl_id annot.type
## 1 ENSG00000304823 hg19_genes_1to5kb
## 2 ENSG00000304823 hg19_genes_1to5kb
## 3 ENSG00000304823 hg19_genes_1to5kb
## 4 <NA> hg19_genes_1to5kb
## 5 ENSG00000304823 hg19_genes_1to5kb
## 6 ENSG00000304823 hg19_genes_1to5kb
# Subset based on a gene symbol, in this case NOTCH1
notch1_subset = subset(df_dm_annotated, annot.symbol == 'NOTCH1')
print(head(notch1_subset))
## [1] seqnames start end width
## [5] strand DM_status pval diff_meth
## [9] mu0 mu1 annot.seqnames annot.start
## [13] annot.end annot.width annot.strand annot.id
## [17] annot.tx_id annot.gene_id annot.symbol annot.entrez_id
## [21] annot.ensembl_id annot.type
## <0 rows> (or 0-length row.names)
Knowing how many regions fall in each annotation is more useful when compared to what would be expected. The right comparison is a background: the set of regions your data could have come from. For differentially methylated (DM) regions, that is all regions tested for differential methylation. For differentially bound ChIP-seq peaks, it might be a consensus set of all peaks. The background carries the same biases as the data (e.g. mappability, or the CpGs an assay measures), so differences between the two reflect the biology of interest rather than those biases.
Here the regions read in above are all regions tested for DM, and the DM_status column says which are DM. So the data are the DM regions, and the background is all tested regions, which are already annotated. Downstream functions that support a background are summarize_annotations(), plot_annotation(), and plot_categorical(), via their annotated_random argument.
# The data are the DM regions, and the background is all tested regions
dm_sig_annotated = dm_annotated[dm_annotated$DM_status != 'none']
NOTE: randomize_regions() is deprecated. Placing regions uniformly at random across the genome ignores those biases, so nearly any data look “enriched” in genic and CpG annotations compared to random regions. When there is no natural background, the permutation test regioneR::permTest() in the regioneR package, with an appropriate mask, gives a null distribution and a p-value. If the data are drawn from a fixed set of candidate regions, regioneR::resampleRegions() samples from that set.
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_numerical() and/or summarize_categorical() may be used. Using a background is only available for summarize_annotations().
# Find the number of regions per annotation type
dm_annsum = summarize_annotations(
annotated_regions = dm_annotated,
quiet = TRUE)
print(dm_annsum)
## # A tibble: 14 × 2
## annot.type n
## <chr> <int>
## 1 hg19_H3K4me3_Gm12878 747
## 2 hg19_cpg_inter 905
## 3 hg19_cpg_islands 848
## 4 hg19_cpg_shelves 46
## 5 hg19_cpg_shores 341
## 6 hg19_custom_ezh2 7
## 7 hg19_genes_1to5kb 595
## 8 hg19_genes_3UTRs 69
## 9 hg19_genes_5UTRs 320
## 10 hg19_genes_exons 790
## 11 hg19_genes_intergenic 181
## 12 hg19_genes_intronexonboundaries 571
## 13 hg19_genes_introns 1452
## 14 hg19_genes_promoters 783
# Find the number of DM regions per annotation type
# and the number of background regions per annotation type
dm_annsum_bg = summarize_annotations(
annotated_regions = dm_sig_annotated,
annotated_random = dm_annotated,
quiet = TRUE)
print(dm_annsum_bg)
## # A tibble: 28 × 3
## # Groups: data_type [2]
## data_type annot.type n
## <chr> <chr> <int>
## 1 Background hg19_H3K4me3_Gm12878 747
## 2 Background hg19_cpg_inter 905
## 3 Background hg19_cpg_islands 848
## 4 Background hg19_cpg_shelves 46
## 5 Background hg19_cpg_shores 341
## 6 Background hg19_custom_ezh2 7
## 7 Background hg19_genes_1to5kb 595
## 8 Background hg19_genes_3UTRs 69
## 9 Background hg19_genes_5UTRs 320
## 10 Background hg19_genes_exons 790
## # ℹ 18 more rows
# 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)
## # A tibble: 8,535 × 5
## # Groups: annot.type [14]
## annot.type annot.id n mean sd
## <chr> <chr> <int> <dbl> <dbl>
## 1 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27530 1 0.0701 NA
## 2 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27531 8 1.28 3.78
## 3 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27532 2 13.4 5.11
## 4 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27534 10 0.526 0.975
## 5 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27535 8 0.407 0.923
## 6 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27543 2 -0.0530 0.0749
## 7 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27544 11 0.192 0.427
## 8 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27545 2 2.80 10.1
## 9 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27549 2 -0.811 1.75
## 10 hg19_H3K4me3_Gm12878 H3K4me3_Gm12878:27555 1 -1.50 NA
## # ℹ 8,525 more rows
# 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)
## # A tibble: 40 × 3
## # Groups: annot.type [14]
## annot.type DM_status n
## <chr> <chr> <int>
## 1 hg19_H3K4me3_Gm12878 hyper 78
## 2 hg19_H3K4me3_Gm12878 hypo 8
## 3 hg19_H3K4me3_Gm12878 none 661
## 4 hg19_cpg_inter hyper 32
## 5 hg19_cpg_inter hypo 90
## 6 hg19_cpg_inter none 783
## 7 hg19_cpg_islands hyper 151
## 8 hg19_cpg_islands hypo 4
## 9 hg19_cpg_islands none 693
## 10 hg19_cpg_shelves hyper 2
## # ℹ 30 more rows
summarize_genes() gives one row per gene, from the gene annotations with a gene ID (e.g. promoters, 1-5Kb, UTRs, exons, and introns). A region counts once toward each gene it is annotated to, and once toward each of the gene’s annotation types. The over columns are summarized with their mean, median, and standard deviation, and the categories of the by column are counted. Use format = 'long' for one row per gene and annotation type instead.
# One row per gene, with the number of regions in each gene annotation
# type, the number of regions in each DM classification, and the mean,
# median, and standard deviation of the methylation difference
dm_genesum = summarize_genes(
annotated_regions = dm_annotated,
over = 'diff_meth',
by = 'DM_status',
quiet = TRUE)
print(dm_genesum)
## # A tibble: 166 × 18
## gene_id symbol entrez_id ensembl_id n_regions n_1to5kb n_promoters n_5UTRs
## <chr> <chr> <chr> <chr> <int> <int> <int> <int>
## 1 23081 KDM4C 23081 ENSG00000… 41 4 10 7
## 2 5789 PTPRD 5789 ENSG00000… 41 2 8 5
## 3 1271 CNTFR 1271 ENSG00000… 40 15 9 2
## 4 4507 MTAP 4507 ENSG00000… 40 3 5 1
## 5 23189 KANK1 23189 ENSG00000… 37 3 8 5
## 6 169792 GLIS3 169792 ENSG00000… 35 5 15 8
## 7 4781 NFIB 4781 ENSG00000… 34 15 17 3
## 8 158038 LINGO2 158038 ENSG00000… 33 7 9 4
## 9 92949 ADAMTS… 92949 ENSG00000… 31 0 3 1
## 10 105375968 LINC02… 105375968 ENSG00000… 31 1 9 0
## # ℹ 156 more rows
## # ℹ 10 more variables: n_exons <int>, n_intronexonboundaries <int>,
## # n_introns <int>, n_3UTRs <int>, n_hyper <int>, n_hypo <int>, n_none <int>,
## # diff_meth_mean <dbl>, diff_meth_median <dbl>, diff_meth_sd <dbl>
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 (print), saved (ggsave), or modified with additional ggplot2 code.
# 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)
Figure 1: Number of DM regions per annotation
The plot_annotation() can also use an annotated background in the annotated_random argument to plot the number of background regions per annotation type next to the number of data regions.
# View the number of DM regions per annotation and include the
# annotation of the background of all tested 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_wbg = plot_annotation(
annotated_regions = dm_sig_annotated,
annotated_random = dm_annotated,
annotation_order = annots_order,
plot_title = 'Dist. of DM Sites (with background)',
x_label = 'Annotations',
y_label = 'Count')
print(dm_vs_kg_annotations_wbg)
Figure 2: Number of DM regions per annotation with the background of all tested regions
# 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)
Figure 3: Number of DM regions per pair of annotations
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)
Figure 4: Methylation Rates in Group 0 for Regions Over DM Status
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 = 'Methylation Differences by Annotation and DM Status',
x_label = 'Methylation Difference')
print(dm_vs_regions_annot2)
Figure 5: Methylation Differences for Regions Over DM Status and Annotation Type
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)
Figure 6: Methylation Rates in Regions Over DM Status in Group 0 vs Group 1
The 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)
Figure 7: Group 0 methylation Rates in Regions in promoters, CpG islands, and both
# 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)
Figure 8: Differential methylation classification with counts of CpG annotations
# 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)
Figure 9: Differential methylation classification with proportion of CpG annotations
As with plot_annotation() one may add an annotated background to the annotated_random parameter of plot_categorical(). The result is a Background bar representing the distribution of the background regions for the categorical variable used for fill. NOTE: A background can only be added when fill = 'annot.type'. In the plots above the data include all tested regions, so the All bar already plays this role. Below, the data are only the DM regions, so the All bar is all DM regions and the Background bar is all tested regions.
# Add in the background of all tested regions for the "Background" bar
# Make a barplot of the DM classes where each bar
# is composed of the *proportion* of CpG annotations, and
# includes "All" DM regions and the "Background" of all
# regions tested for DM.
dm_vs_cpg_cat_bg = plot_categorical(
annotated_regions = dm_sig_annotated, annotated_random = dm_annotated,
x='DM_status', fill='annot.type',
x_order = c('hyper', 'hypo'), 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_bg)
Figure 10: Differential methylation classification with proportion of CpG annotations and the background
# 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)
Figure 11: Basic gene annotations with proportions of DM classification
sessionInfo()
## R version 4.6.1 Patched (2026-06-24 r90190)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.7.8
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.81.1
## [3] BiocIO_1.23.3 Biostrings_2.81.9
## [5] XVector_0.53.0 rtracklayer_1.73.0
## [7] GenomicRanges_1.65.4 Seqinfo_1.3.2
## [9] IRanges_2.47.5 S4Vectors_0.51.10
## [11] BiocGenerics_0.59.12 generics_0.1.4
## [13] annotatr_1.39.17 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1
## [2] dplyr_1.2.1
## [3] farver_2.1.2
## [4] blob_1.3.0
## [5] filelock_1.0.3
## [6] S7_0.2.2
## [7] bitops_1.1-0
## [8] fastmap_1.2.0
## [9] RCurl_1.98-1.20
## [10] BiocFileCache_3.3.0
## [11] GenomicAlignments_1.49.2
## [12] XML_3.99-0.24
## [13] digest_0.6.39
## [14] lifecycle_1.0.5
## [15] KEGGREST_1.53.6
## [16] RSQLite_3.53.3
## [17] magrittr_2.0.5
## [18] compiler_4.6.1
## [19] rlang_1.3.0
## [20] sass_0.4.10
## [21] tools_4.6.1
## [22] utf8_1.2.6
## [23] yaml_2.3.12
## [24] knitr_1.52
## [25] labeling_0.4.3
## [26] bit_4.6.0
## [27] curl_8.0.0
## [28] plyr_1.8.9
## [29] regioneR_1.45.0
## [30] RColorBrewer_1.1-3
## [31] BiocParallel_1.47.0
## [32] withr_3.0.3
## [33] purrr_1.2.2
## [34] grid_4.6.1
## [35] ggplot2_4.0.3
## [36] scales_1.4.0
## [37] tinytex_0.61
## [38] dichromat_2.0-1
## [39] cli_3.6.6
## [40] rmarkdown_2.32
## [41] crayon_1.5.3
## [42] otel_0.2.0
## [43] reshape2_1.4.5
## [44] httr_1.4.9
## [45] tzdb_0.5.0
## [46] rjson_0.2.23
## [47] BiocBaseUtils_1.15.1
## [48] DBI_1.3.0
## [49] cachem_1.1.0
## [50] stringr_1.6.0
## [51] parallel_4.6.1
## [52] AnnotationDbi_1.75.2
## [53] BiocManager_1.30.27
## [54] restfulr_0.0.17
## [55] matrixStats_1.5.0
## [56] vctrs_0.7.3
## [57] jsonlite_2.0.0
## [58] bookdown_0.48
## [59] hms_1.1.4
## [60] bit64_4.8.6
## [61] magick_2.9.1
## [62] GenomicFeatures_1.65.0
## [63] jquerylib_0.1.4
## [64] glue_1.8.1
## [65] codetools_0.2-20
## [66] stringi_1.8.9
## [67] gtable_0.3.6
## [68] BiocVersion_3.24.0
## [69] GenomeInfoDb_1.49.1
## [70] UCSC.utils_1.9.0
## [71] tibble_3.3.1
## [72] pillar_1.11.1
## [73] rappdirs_0.3.4
## [74] htmltools_0.5.9
## [75] R6_2.6.1
## [76] dbplyr_2.6.0
## [77] httr2_1.3.0
## [78] evaluate_1.0.5
## [79] Biobase_2.73.2
## [80] readr_2.2.0
## [81] AnnotationHub_4.3.2
## [82] png_0.1-9
## [83] Rsamtools_2.29.0
## [84] cigarillo_1.3.1
## [85] TxDb.Hsapiens.UCSC.hg19.knownGene_3.22.1
## [86] memoise_2.0.1
## [87] bslib_0.12.0
## [88] Rcpp_1.1.2
## [89] org.Hs.eg.db_3.23.1
## [90] xfun_0.61
## [91] pkgconfig_2.0.3