RNAmodR: AlkAnilineSeq
2020-10-27
Abstract
detection of m7G, m3C and D modification by AlkAnilineSeq1 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.
## snapshotDate(): 2020-10-02library(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.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-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] Rsamtools_2.6.0 RNAmodR.Data_1.3.0
## [3] ExperimentHubData_1.16.0 AnnotationHubData_1.20.0
## [5] futile.logger_1.4.3 ExperimentHub_1.16.0
## [7] AnnotationHub_2.22.0 BiocFileCache_1.14.0
## [9] dbplyr_1.4.4 RNAmodR.AlkAnilineSeq_1.4.0
## [11] RNAmodR_1.4.0 Modstrings_1.6.0
## [13] Biostrings_2.58.0 XVector_0.30.0
## [15] rtracklayer_1.50.0 GenomicRanges_1.42.0
## [17] GenomeInfoDb_1.26.0 IRanges_2.24.0
## [19] S4Vectors_0.28.0 BiocGenerics_0.36.0
## [21] BiocStyle_2.18.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.10 Hmisc_4.4-1
## [3] plyr_1.8.6 lazyeval_0.2.2
## [5] splines_4.0.3 BiocParallel_1.24.0
## [7] ggplot2_3.3.2 digest_0.6.27
## [9] ensembldb_2.14.0 htmltools_0.5.0
## [11] magick_2.5.0 magrittr_1.5
## [13] checkmate_2.0.0 memoise_1.1.0
## [15] BSgenome_1.58.0 cluster_2.1.0
## [17] ROCR_1.0-11 matrixStats_0.57.0
## [19] askpass_1.1 prettyunits_1.1.1
## [21] jpeg_0.1-8.1 colorspace_1.4-1
## [23] blob_1.2.1 rappdirs_0.3.1
## [25] xfun_0.18 dplyr_1.0.2
## [27] jsonlite_1.7.1 crayon_1.3.4
## [29] RCurl_1.98-1.2 graph_1.68.0
## [31] survival_3.2-7 VariantAnnotation_1.36.0
## [33] glue_1.4.2 gtable_0.3.0
## [35] zlibbioc_1.36.0 DelayedArray_0.16.0
## [37] scales_1.1.1 futile.options_1.0.1
## [39] DBI_1.1.0 Rcpp_1.0.5
## [41] xtable_1.8-4 progress_1.2.2
## [43] htmlTable_2.1.0 foreign_0.8-80
## [45] bit_4.0.4 OrganismDbi_1.32.0
## [47] Formula_1.2-4 AnnotationForge_1.32.0
## [49] getopt_1.20.3 htmlwidgets_1.5.2
## [51] httr_1.4.2 RColorBrewer_1.1-2
## [53] ellipsis_0.3.1 farver_2.0.3
## [55] pkgconfig_2.0.3 XML_3.99-0.5
## [57] Gviz_1.34.0 nnet_7.3-14
## [59] labeling_0.4.2 tidyselect_1.1.0
## [61] rlang_0.4.8 reshape2_1.4.4
## [63] later_1.1.0.1 AnnotationDbi_1.52.0
## [65] munsell_0.5.0 biocViews_1.58.0
## [67] BiocVersion_3.12.0 tools_4.0.3
## [69] generics_0.0.2 RSQLite_2.2.1
## [71] evaluate_0.14 stringr_1.4.0
## [73] fastmap_1.0.1 yaml_2.2.1
## [75] knitr_1.30 bit64_4.0.5
## [77] purrr_0.3.4 AnnotationFilter_1.14.0
## [79] RBGL_1.66.0 mime_0.9
## [81] formatR_1.7 xml2_1.3.2
## [83] biomaRt_2.46.0 compiler_4.0.3
## [85] rstudioapi_0.11 curl_4.3
## [87] png_0.1-7 interactiveDisplayBase_1.28.0
## [89] tibble_3.0.4 stringi_1.5.3
## [91] highr_0.8 GenomicFeatures_1.42.0
## [93] lattice_0.20-41 ProtGenerics_1.22.0
## [95] Matrix_1.2-18 vctrs_0.3.4
## [97] stringdist_0.9.6.3 BiocCheck_1.26.0
## [99] pillar_1.4.6 lifecycle_0.2.0
## [101] optparse_1.6.6 RUnit_0.4.32
## [103] BiocManager_1.30.10 data.table_1.13.2
## [105] bitops_1.0-6 httpuv_1.5.4
## [107] colorRamps_2.3 R6_2.4.1
## [109] latticeExtra_0.6-29 bookdown_0.21
## [111] promises_1.1.1 gridExtra_2.3
## [113] codetools_0.2-16 lambda.r_1.2.4
## [115] dichromat_2.0-0 assertthat_0.2.1
## [117] SummarizedExperiment_1.20.0 openssl_1.4.3
## [119] rBiopaxParser_2.30.0 GenomicAlignments_1.26.0
## [121] GenomeInfoDbData_1.2.4 hms_0.5.3
## [123] grid_4.0.3 rpart_4.1-15
## [125] rmarkdown_2.5 MatrixGenerics_1.2.0
## [127] biovizBase_1.38.0 Biobase_2.50.0
## [129] shiny_1.5.0 base64enc_0.1-3References
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.