RNAmodR: AlkAnilineSeq
2023-04-25
Abstract
detection of m7G, m3C and D modification by AlkAnilineSeq
1 Introduction
7-methyl guanosine (m7G), 3-methyl cytidine (m3C) and Dihydrouridine (D) are commonly found in rRNA and tRNA and can be detected classically by primer extension analysis. However, since the modifications do not interfere with Watson-Crick base pairing, a specific chemical treatment needs to be employed to cause strand breaks specifically at the modified positions. Initially, this involved a sodium borhydride treatment to create abasic sites and cleaving the RNA at abasic sites with aniline.
This classical protocol was converted to a high throughput sequencing method call AlkAnilineSeq and allows modified position be detected by an accumulation of 5’-ends at the N+1 position (Marchand et al. 2018). It was found, that m3C is susceptible to this treatment, which allows m7G, m3C and D to be detected by the same method from the same data sets, since the identify of the unmodified nucleotide informs about the three modified nucleotides.
The ModAlkAnilineSeq class uses the the NormEnd5SequenceData class to store
and aggregate data along the transcripts. The calculated scores follow the
nomenclature of (Marchand et al. 2018) with the names scoreNC
(default) and scoreSR.
## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationDbi'## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationForge'## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'ExperimentHubData'library(rtracklayer)
library(GenomicRanges)
library(RNAmodR.AlkAnilineSeq)
library(RNAmodR.Data)2 Example workflow
The example workflow is limited to 18S rRNA and some tRNA from S.cerevisiae.
As annotation data either a gff file or a TxDb object and for sequence data
a fasta file or a BSgenome object can be used. The data is provided as bam
files.
annotation <- GFF3File(RNAmodR.Data.example.AAS.gff3())
sequences <- RNAmodR.Data.example.AAS.fasta()
files <- list("wt" = c(treated = RNAmodR.Data.example.wt.1(),
treated = RNAmodR.Data.example.wt.2(),
treated = RNAmodR.Data.example.wt.3()),
"Bud23del" = c(treated = RNAmodR.Data.example.bud23.1(),
treated = RNAmodR.Data.example.bud23.2()),
"Trm8del" = c(treated = RNAmodR.Data.example.trm8.1(),
treated = RNAmodR.Data.example.trm8.2()))The analysis is triggered by the construction of a ModSetAlkAnilineSeq object.
Internally parallelization is used via the BiocParallel package, which would
allow optimization depending on number/size of input files (number of samples,
number of replicates, number of transcripts, etc).
msaas <- ModSetAlkAnilineSeq(files, annotation = annotation, sequences = sequences)## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ... OKmsaas## ModSetAlkAnilineSeq of length 3
## names(3): wt Bud23del Trm8del
## | Modification type(s): m7G / m3C / D
## wt Bud23del Trm8del
## | Modifications found: yes (9) yes (8) yes (7)
## | Settings:
## minCoverage minReplicate find.mod minLength minSignal minScoreNC
## <integer> <integer> <logical> <integer> <integer> <integer>
## wt 10 1 TRUE 9 10 50
## Bud23del 10 1 TRUE 9 10 50
## Trm8del 10 1 TRUE 9 10 50
## minScoreSR minScoreBaseScore scoreOperator
## <numeric> <numeric> <character>
## wt 0.5 0.9 &
## Bud23del 0.5 0.9 &
## Trm8del 0.5 0.9 &As expected the m7G1575 is missing from the Bud23del samples.
mod <- modifications(msaas)
lapply(mod,head, n = 2L)## $wt
## GRanges object with 2 ranges and 6 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr1 1575 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## [2] chr3 46 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## score scoreSR Parent
## <numeric> <numeric> <character>
## [1] 162.228 0.984209 1
## [2] 373.773 0.841166 3
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengths
##
## $Bud23del
## GRanges object with 2 ranges and 6 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr3 46 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## [2] chr5 50 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## score scoreSR Parent
## <numeric> <numeric> <character>
## [1] 254.6403 0.858101 3
## [2] 86.3556 0.605249 5
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengths
##
## $Trm8del
## GRanges object with 2 ranges and 6 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr1 1575 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## [2] chr3 37 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## score scoreSR Parent
## <numeric> <numeric> <character>
## [1] 117.2479 0.98729 1
## [2] 69.9604 0.97953 3
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengths3 Visualizing the results
As outlined in the RNAmodR package we can compare the samples using the
plotCompareByCoord to prepare a heatmap. For this we select some position
from the found modifications. In addition we prepare an alias table.
coord <- mod[[1L]]
alias <- data.frame(tx_id = c(1L,3L,5L,6L,7L,8L,10L,11L),
name = c("18S rRNA","tF(GAA)B","tG(GCC)B","tT(AGT)B",
"tQ(TTG)B","tC(GCA)B","tS(CGA)C","tV(AAC)E1"),
stringsAsFactors = FALSE)plotCompareByCoord(msaas, coord, score = "scoreSR", alias = alias,
normalize = TRUE)
Figure 1: Heatmap showing Stop ratio scores for detected m7G, m3C and D positions
plotCompareByCoord(msaas, coord[1L], score = "scoreSR", alias = alias)
Figure 2: Heatmap showing Stop ratio scores for detected m7G1575 on the 18S rRNA
In addition, the aggregate data along the transcript visualized as well.
plotData(msaas, "1", from = 1550L, to = 1600L)
Figure 3: Stop ratio scores for detected m7G1575 on the 18S rRNA plotted as bar plots along the sequence
This includes raw data as well.
plotData(msaas[1L:2L], "1", from = 1550L, to = 1600L, showSequenceData = TRUE)
Figure 4: Stop ratio scores for detected m7G1575 on the 18S rRNA plotted as bar plots along the sequence
The raw sequence data is shown by setting `showSequenceData = TRUE.
4 Session info
sessionInfo()## R version 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] Rsamtools_2.16.0 RNAmodR.Data_1.13.0
## [3] ExperimentHubData_1.26.0 AnnotationHubData_1.30.0
## [5] futile.logger_1.4.3 ExperimentHub_2.8.0
## [7] AnnotationHub_3.8.0 BiocFileCache_2.8.0
## [9] dbplyr_2.3.2 RNAmodR.AlkAnilineSeq_1.14.0
## [11] RNAmodR_1.14.0 Modstrings_1.16.0
## [13] Biostrings_2.68.0 XVector_0.40.0
## [15] rtracklayer_1.60.0 GenomicRanges_1.52.0
## [17] GenomeInfoDb_1.36.0 IRanges_2.34.0
## [19] S4Vectors_0.38.0 BiocGenerics_0.46.0
## [21] BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] later_1.3.0 BiocIO_1.10.0
## [3] bitops_1.0-7 filelock_1.0.2
## [5] tibble_3.2.1 graph_1.78.0
## [7] XML_3.99-0.14 rpart_4.1.19
## [9] lifecycle_1.0.3 lattice_0.21-8
## [11] ensembldb_2.24.0 OrganismDbi_1.42.0
## [13] backports_1.4.1 magrittr_2.0.3
## [15] Hmisc_5.0-1 sass_0.4.5
## [17] rmarkdown_2.21 jquerylib_0.1.4
## [19] yaml_2.3.7 httpuv_1.6.9
## [21] RUnit_0.4.32 Gviz_1.44.0
## [23] DBI_1.1.3 RColorBrewer_1.1-3
## [25] zlibbioc_1.46.0 purrr_1.0.1
## [27] AnnotationFilter_1.24.0 biovizBase_1.48.0
## [29] RCurl_1.98-1.12 nnet_7.3-18
## [31] VariantAnnotation_1.46.0 rappdirs_0.3.3
## [33] GenomeInfoDbData_1.2.10 AnnotationForge_1.42.0
## [35] codetools_0.2-19 DelayedArray_0.26.0
## [37] xml2_1.3.3 tidyselect_1.2.0
## [39] farver_2.1.1 matrixStats_0.63.0
## [41] base64enc_0.1-3 GenomicAlignments_1.36.0
## [43] jsonlite_1.8.4 ellipsis_0.3.2
## [45] Formula_1.2-5 tools_4.3.0
## [47] progress_1.2.2 stringdist_0.9.10
## [49] Rcpp_1.0.10 glue_1.6.2
## [51] gridExtra_2.3 xfun_0.39
## [53] MatrixGenerics_1.12.0 dplyr_1.1.2
## [55] withr_2.5.0 formatR_1.14
## [57] BiocManager_1.30.20 fastmap_1.1.1
## [59] latticeExtra_0.6-30 fansi_1.0.4
## [61] digest_0.6.31 R6_2.5.1
## [63] mime_0.12 colorspace_2.1-0
## [65] jpeg_0.1-10 dichromat_2.0-0.1
## [67] biomaRt_2.56.0 RSQLite_2.3.1
## [69] utf8_1.2.3 generics_0.1.3
## [71] data.table_1.14.8 prettyunits_1.1.1
## [73] httr_1.4.5 htmlwidgets_1.6.2
## [75] pkgconfig_2.0.3 gtable_0.3.3
## [77] blob_1.2.4 htmltools_0.5.5
## [79] bookdown_0.33 RBGL_1.76.0
## [81] ProtGenerics_1.32.0 scales_1.2.1
## [83] Biobase_2.60.0 png_0.1-8
## [85] colorRamps_2.3.1 knitr_1.42
## [87] lambda.r_1.2.4 rstudioapi_0.14
## [89] reshape2_1.4.4 rjson_0.2.21
## [91] checkmate_2.1.0 curl_5.0.0
## [93] biocViews_1.68.0 cachem_1.0.7
## [95] stringr_1.5.0 BiocVersion_3.17.1
## [97] parallel_4.3.0 foreign_0.8-84
## [99] AnnotationDbi_1.62.0 restfulr_0.0.15
## [101] pillar_1.9.0 grid_4.3.0
## [103] vctrs_0.6.2 promises_1.2.0.1
## [105] xtable_1.8-4 cluster_2.1.4
## [107] htmlTable_2.4.1 evaluate_0.20
## [109] magick_2.7.4 GenomicFeatures_1.52.0
## [111] cli_3.6.1 compiler_4.3.0
## [113] futile.options_1.0.1 rlang_1.1.0
## [115] crayon_1.5.2 labeling_0.4.2
## [117] interp_1.1-4 plyr_1.8.8
## [119] stringi_1.7.12 deldir_1.0-6
## [121] BiocParallel_1.34.0 BiocCheck_1.36.0
## [123] munsell_0.5.0 lazyeval_0.2.2
## [125] Matrix_1.5-4 BSgenome_1.68.0
## [127] hms_1.1.3 bit64_4.0.5
## [129] ggplot2_3.4.2 KEGGREST_1.40.0
## [131] shiny_1.7.4 highr_0.10
## [133] SummarizedExperiment_1.30.0 interactiveDisplayBase_1.38.0
## [135] ROCR_1.0-11 memoise_2.0.1
## [137] bslib_0.4.2 bit_4.0.5References
Marchand, Virginie, Lilia Ayadi, Felix G. M. Ernst, Jasmin Hertler, Valérie Bourguignon-Igel, Adeline Galvanin, Annika Kotter, Mark Helm, Denis L. J. Lafontaine, and Yuri Motorin. 2018. “AlkAniline-Seq: Profiling of m7G and m3C Rna Modifications at Single Nucleotide Resolution.” Angewandte Chemie International Edition 57 (51): 16785–90. https://doi.org/10.1002/anie.201810946.