EpiTxDb: Storing and accessing epitranscriptomic information using the AnnotationDbi interface
2020-04-28
Package
EpiTxDb 1.0.0
1 Installation
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("EpiTxDb","EpiTxDb.Hs.hg38"))2 Introduction
The epitranscriptome includes all post-transcriptional modifications of the RNA and describes and additional layer of information encoded on RNA. Like the term epigenome it is not about a change in nucleotide sequences, but the addition of functional elements through modifications.
With the development of high throughput detection strategies for specific RNA modifications, such as miCLIP and Pseudo-Seq amongst other, a large number of modified positions have been identified and were summarized via the RMBase project [Xuan et al. (2017);Sun et al. (2015)] project.
To make these information avaialble within the Bioconductor universe EpiTxDb
was developed, which facilitates the storage of epitranscriptomic information.
More specifically, it can keep track of modification identity, position, the
enzyme for introducing it on the RNA, a specifier which determines the position
on the RNA to be modified and the literature references each modification is
associated with.
3 Getting started
library(EpiTxDb)
library(EpiTxDb.Hs.hg38)The class EpiTxDb is the class for storing the epitranscriptomic data. It
inherits the inner workings of AnnotationDb class from the AnnotationDbi
package.
As an example for the vignette the snoRNAdb data [Lestrade and Weber (2006)]
from the EpiTxDb.Hs.hg38 package will be used. The data is stored in the
AnnotationHub and is downloaded and cached upon the first request.
etdb <- EpiTxDb.Hs.hg38.snoRNAdb()## snapshotDate(): 2020-03-31## loading from cacheetdb## EpiTxDb object:
## # Db type: EpiTxDb
## # Supporting package: EpiTxDb
## # Data source: snoRNAdb
## # Organism: Homo sapiens
## # Genome: hg38
## # Coordinates: per Transcript
## # Nb of modifications: 235
## # Db created by: EpiTxDb package from Bioconductor
## # Creation time: 2020-02-26 10:34:30 +0100 (Wed, 26 Feb 2020)
## # EpiTxDb version at creation time: 0.99.0
## # RSQLite version at creation time: 2.2.0
## # DBSCHEMAVERSION: 1.0As expected for an AnnotationDb class the general accessors are available.
keytypes(etdb)## [1] "MODID" "MODNAME" "MODSTRAND" "MODTYPE"
## [5] "REF" "REFTYPE" "RXENSEMBL" "RXENSEMBLTRANS"
## [9] "RXENTREZID" "RXGENENAME" "SNID" "SNNAME"
## [13] "SPECENSEMBL" "SPECENSEMBLTRANS" "SPECENTREZID" "SPECGENENAME"
## [17] "SPECTYPE"columns(etdb)## [1] "MODEND" "MODID" "MODNAME" "MODSTART"
## [5] "MODSTRAND" "MODTYPE" "REF" "REFID"
## [9] "REFTYPE" "RXENSEMBL" "RXENSEMBLTRANS" "RXENTREZID"
## [13] "RXGENENAME" "RXID" "RXRANK" "SNID"
## [17] "SNNAME" "SPECENSEMBL" "SPECENSEMBLTRANS" "SPECENTREZID"
## [21] "SPECGENENAME" "SPECID" "SPECTYPE"head(keys(etdb, "MODID"))## [1] "1" "2" "3" "4" "5" "6"select(etdb, keys = "1",
columns = c("MODNAME","MODTYPE","MODSTART","MODSTRAND","SNNAME",
"RXGENENAME","SPECTYPE","SPECGENENAME"),
keytype = "MODID")## 'select()' returned 1:1 mapping between keys and columnsThe columns with the prefix RX or SPEC reference the reaction enzyme and the
location specifier. This can be the same information, but for ribosomal
modifications from the snoRNAdb it is of course fibrillarin and a snoRNA.
In addition the following accessor for metadata are available as well.
species(etdb)## [1] "Homo sapiens"organism(etdb)## [1] "Homo sapiens"seqlevels(etdb)## [1] "NR_003285" "NR_003286" "NR_003287" "NR_004430" "NR_002716" "NR_003925"
## [7] "NR_002756" "NR_004394" "NR_029422"4 Accessing RNA modifications
The specialized accessors are modifications() and modificationsBy().
modifications() allows for filtering results, whereas modificationsBy()
returns all the modifications in batches separated by certain information.
modifications(etdb, columns = c("mod_id","mod_type","mod_name",
"rx_genename","spec_genename",
"ref_type","ref"),
filter = list(mod_id = 1:3))## GRanges object with 3 ranges and 7 metadata columns:
## seqnames ranges strand | mod_id mod_type mod_name
## <Rle> <IRanges> <Rle> | <integer> <character> <character>
## [1] NR_003285 14 + | 1 Um Um_14
## [2] NR_003285 55 + | 2 Y Y_55
## [3] NR_003285 69 + | 3 Y Y_69
## rx_genename spec_genename ref_type
## <CharacterList> <CharacterList> <CharacterList>
## [1] fibrillarin SNORD71 PMID
## [2] dyskerin pseudouridine synthase 1 SNORA72 PMID
## [3] dyskerin pseudouridine synthase 1 SNORA69 PMID
## ref
## <CharacterList>
## [1] 16381836
## [2] 16381836
## [3] 16381836
## -------
## seqinfo: 9 sequences from hg38 genome; no seqlengths# split by sequence name, usually a transcipt identifier
modificationsBy(etdb, by = "seqnames")## GRangesList object of length 9:
## $NR_003285
## GRanges object with 4 ranges and 3 metadata columns:
## seqnames ranges strand | mod_id mod mod_name
## <Rle> <IRanges> <Rle> | <integer> <character> <character>
## [1] NR_003285 14 + | 1 Um Um_14
## [2] NR_003285 55 + | 2 Y Y_55
## [3] NR_003285 69 + | 3 Y Y_69
## [4] NR_003285 75 + | 4 Gm Gm_75
## -------
## seqinfo: 9 sequences from hg38 genome; no seqlengths
##
## ...
## <8 more elements># split modification type
modificationsBy(etdb, by = "modtype")## GRangesList object of length 5:
## $Am
## GRanges object with 39 ranges and 3 metadata columns:
## seqnames ranges strand | mod mod_id mod_name
## <Rle> <IRanges> <Rle> | <character> <integer> <character>
## [1] NR_003286 27 + | Am 5 Am_27
## [2] NR_003286 99 + | Am 9 Am_99
## [3] NR_003286 159 + | Am 15 Am_159
## [4] NR_003286 166 + | Am 16 Am_166
## [5] NR_003286 468 + | Am 25 Am_468
## ... ... ... ... . ... ... ...
## [35] NR_002716 30 + | Am 207 Am_30
## [36] NR_003925 65 + | Am 216 Am_65
## [37] NR_004394 47 + | Am 226 Am_47
## [38] NR_004394 53 + | Am 227 Am_53
## [39] NR_004394 53 + | Am 228 Am_53
## -------
## seqinfo: 9 sequences from hg38 genome; no seqlengths
##
## ...
## <4 more elements>5 Shifting coordinates from genomic to transcritomic
Since epitranscriptomic modifications by their nature refere to each individual it might be of interest to shift the coordinates transctipt coordinates keeping in mind that with transcript variants multiple options exist for each of the transcript maturation process.
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(BSgenome.Hsapiens.UCSC.hg38)txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
seqlevels(txdb) <- "chr1"
bs <- BSgenome.Hsapiens.UCSC.hg38
etdb <- EpiTxDb.Hs.hg38.RMBase()## snapshotDate(): 2020-03-31## loading from cachetx <- exonsBy(txdb)
mod <- modifications(etdb, filter = list(sn_name = "chr1"))
length(mod)## [1] 47275In the following example we will focus on shifting the coordinates to individual
mature transcripts. However, keep in mind, that premature transcript might be
of interest as well and this can be controlled via the tx arguments of
shiftGenomicToTranscript()
mod_tx <- shiftGenomicToTranscript(mod, tx)## Warning: Coordinates for 905 ranges of 'subject' not found:
## 'chr1:14662:-','chr1:14668:-','chr1:14766:-','chr1:139005:-','chr1:139020:-','chr1:139040:-','chr1:139127:-','chr1:629454:+','chr1:629456:+','chr1:629465:+'
## and more ...length(mod_tx)## [1] 157694Due to multiple matches for each transcript variant the number of modifications has increased.
With the we can plot the relative positions of modifications by type on
chr1 transcripts.
mod_tx <- split(mod_tx,seqnames(mod_tx))
names <- Reduce(intersect,list(names(mod_tx),names(tx)))
# Getting the corresponding 5'-UTR and 3'-UTR annotations
fp <- fiveUTRsByTranscript(txdb)
tp <- threeUTRsByTranscript(txdb)
tx <- tx[names]
mod_tx <- mod_tx[names]
fp_m <- match(names,names(fp))
fp_m <- fp_m[!is.na(fp_m)]
tp_m <- match(names,names(tp))
tp_m <- tp_m[!is.na(tp_m)]
fp <- fp[fp_m]
tp <- tp[tp_m]
# Getting lengths of transcripts, 5'-UTR and 3'-UTR
tx_lengths <- sum(width(tx))
fp_lengths <- rep(0L,length(tx))
names(fp_lengths) <- names
fp_lengths[names(fp)] <- sum(width(fp))
tp_lengths <- rep(0L,length(tx))
names(tp_lengths) <- names
tp_lengths[names(tp)] <- sum(width(tp))
# Rescale modifications
# CDS start is at position 1L and cds end at position 1000L
from <- IRanges(fp_lengths+1L, tx_lengths - tp_lengths)
to <- IRanges(1L,1000L)
mod_rescale <- rescale(mod_tx, to, from)
# Construct result data.frame
rel_pos <- data.frame(mod_type = unlist(mcols(mod_rescale,level="within")[,"mod_type"]),
rel_pos = unlist(start(mod_rescale)))
rel_pos <- rel_pos[rel_pos$rel_pos < 1500 & rel_pos$rel_pos > -500,]library(ggplot2)
ggplot(rel_pos[rel_pos$mod_type %in% c("m6A","m1A","Y"),],
aes(x = rel_pos, colour = mod_type)) +
geom_density()
6 Session info
sessionInfo()## R version 4.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.11-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.11-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggplot2_3.3.0
## [2] BSgenome.Hsapiens.UCSC.hg38_1.4.3
## [3] BSgenome_1.56.0
## [4] rtracklayer_1.48.0
## [5] TxDb.Hsapiens.UCSC.hg38.knownGene_3.10.0
## [6] GenomicFeatures_1.40.0
## [7] EpiTxDb.Hs.hg38_0.99.3
## [8] AnnotationHub_2.20.0
## [9] BiocFileCache_1.12.0
## [10] dbplyr_1.4.3
## [11] EpiTxDb_1.0.0
## [12] Modstrings_1.4.0
## [13] Biostrings_2.56.0
## [14] XVector_0.28.0
## [15] AnnotationDbi_1.50.0
## [16] Biobase_2.48.0
## [17] GenomicRanges_1.40.0
## [18] GenomeInfoDb_1.24.0
## [19] IRanges_2.22.0
## [20] S4Vectors_0.26.0
## [21] BiocGenerics_0.34.0
## [22] BiocStyle_2.16.0
##
## loaded via a namespace (and not attached):
## [1] assertive.base_0.0-7 colorspace_1.4-1
## [3] ellipsis_0.3.0 farver_2.0.3
## [5] assertive.sets_0.0-3 bit64_0.9-7
## [7] interactiveDisplayBase_1.26.0 xml2_1.3.2
## [9] assertive.data.uk_0.0-2 codetools_0.2-16
## [11] knitr_1.28 jsonlite_1.6.1
## [13] assertive_0.3-5 assertive.data.us_0.0-2
## [15] Rsamtools_2.4.0 shiny_1.4.0.2
## [17] BiocManager_1.30.10 compiler_4.0.0
## [19] httr_1.4.1 assertthat_0.2.1
## [21] Matrix_1.2-18 fastmap_1.0.1
## [23] tRNA_1.6.0 later_1.0.0
## [25] htmltools_0.4.0 prettyunits_1.1.1
## [27] tools_4.0.0 gtable_0.3.0
## [29] glue_1.4.0 GenomeInfoDbData_1.2.3
## [31] dplyr_0.8.5 rappdirs_0.3.1
## [33] Rcpp_1.0.4.6 vctrs_0.2.4
## [35] assertive.files_0.0-2 assertive.datetimes_0.0-2
## [37] assertive.models_0.0-2 tRNAdbImport_1.6.0
## [39] xfun_0.13 stringr_1.4.0
## [41] Structstrings_1.4.0 mime_0.9
## [43] lifecycle_0.2.0 XML_3.99-0.3
## [45] zlibbioc_1.34.0 scales_1.1.0
## [47] hms_0.5.3 promises_1.1.0
## [49] SummarizedExperiment_1.18.0 assertive.matrices_0.0-2
## [51] assertive.strings_0.0-3 yaml_2.2.1
## [53] curl_4.3 memoise_1.1.0
## [55] biomaRt_2.44.0 stringi_1.4.6
## [57] RSQLite_2.2.0 BiocVersion_3.11.1
## [59] BiocParallel_1.22.0 rlang_0.4.5
## [61] pkgconfig_2.0.3 matrixStats_0.56.0
## [63] bitops_1.0-6 evaluate_0.14
## [65] lattice_0.20-41 assertive.data_0.0-3
## [67] purrr_0.3.4 labeling_0.3
## [69] GenomicAlignments_1.24.0 assertive.properties_0.0-4
## [71] bit_1.1-15.2 tidyselect_1.0.0
## [73] assertive.code_0.0-3 magrittr_1.5
## [75] bookdown_0.18 R6_2.4.1
## [77] magick_2.3 DelayedArray_0.14.0
## [79] DBI_1.1.0 withr_2.2.0
## [81] pillar_1.4.3 assertive.numbers_0.0-2
## [83] RCurl_1.98-1.2 tibble_3.0.1
## [85] crayon_1.3.4 assertive.types_0.0-3
## [87] rmarkdown_2.1 progress_1.2.2
## [89] grid_4.0.0 blob_1.2.1
## [91] digest_0.6.25 xtable_1.8-4
## [93] httpuv_1.5.2 openssl_1.4.1
## [95] munsell_0.5.0 assertive.reflection_0.0-4
## [97] askpass_1.1References
Lestrade, Laurent, and Michel J. Weber. 2006. “snoRNA-LBME-db, a comprehensive database of human H/ACA and C/D box snoRNAs.” Nucleic Acids Research 34 (January):D158–D162. https://doi.org/10.1093/nar/gkj002.
Sun, Wen-Ju, Jun-Hao Li, Shun Liu, Jie Wu, Hui Zhou, Liang-Hu Qu, and Jian-Hua Yang. 2015. “RMBase: a resource for decoding the landscape of RNA modifications from high-throughput sequencing data.” Nucleic Acids Research 44 (D1):D259–D265. https://doi.org/10.1093/nar/gkv1036.
Xuan, Jia-Jia, Wen-Ju Sun, Peng-Hui Lin, Ke-Ren Zhou, Shun Liu, Ling-Ling Zheng, Liang-Hu Qu, and Jian-Hua Yang. 2017. “RMBase v2.0: deciphering the map of RNA modifications from epitranscriptome sequencing data.” Nucleic Acids Research 46 (D1):D327–D334. https://doi.org/10.1093/nar/gkx934.