# RNAmodR: AlkAnilineSeq

Felix G.M. Ernst and Denis L.J. Lafontaine

#### 2021-05-19

#### 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`.

    ## snapshotDate(): 2021-05-05

    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 ... OK
    ## 2. ModAlkAnilineSeq analysis 'Bud23del':
    ## Loading normalized 5'-end position data from BAM files ... OK
    ## Loading Pileup data from BAM files ... OK
    ## Aggregating data and calculating scores ... Starting to search for '3-Methylcytidine', '7-Methylguanosine', 'Dihydrouridine' ... done.
    ## 3. ModAlkAnilineSeq analysis 'Trm8del':
    ## Loading normalized 5'-end position data from BAM files ... OK
    ## Loading Pileup data from BAM files ... OK
    ## Aggregating data and calculating scores ... Starting to search for '3-Methylcytidine', '7-Methylguanosine', 'Dihydrouridine' ... done.
    ## 1. ModAlkAnilineSeq analysis 'wt':
    ## Loading normalized 5'-end position data from BAM files ... OK
    ## Loading Pileup data from BAM files ... OK
    ## Aggregating data and calculating scores ... Starting to search for '3-Methylcytidine', '7-Methylguanosine', 'Dihydrouridine' ... done.

    msaas

    ## 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 seqlengths

# 3 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.1.0 RC (2021-05-10 r80283)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server x64 (build 17763)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] Rsamtools_2.8.0             RNAmodR.Data_1.5.0
    ##  [3] ExperimentHubData_1.18.0    AnnotationHubData_1.22.0
    ##  [5] futile.logger_1.4.3         ExperimentHub_2.0.0
    ##  [7] AnnotationHub_3.0.0         BiocFileCache_2.0.0
    ##  [9] dbplyr_2.1.1                RNAmodR.AlkAnilineSeq_1.6.0
    ## [11] RNAmodR_1.6.0               Modstrings_1.8.0
    ## [13] Biostrings_2.60.0           XVector_0.32.0
    ## [15] rtracklayer_1.52.0          GenomicRanges_1.44.0
    ## [17] GenomeInfoDb_1.28.0         IRanges_2.26.0
    ## [19] S4Vectors_0.30.0            BiocGenerics_0.38.0
    ## [21] BiocStyle_2.20.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] snow_0.4-3                    backports_1.2.1
    ##   [3] Hmisc_4.5-0                   plyr_1.8.6
    ##   [5] lazyeval_0.2.2                splines_4.1.0
    ##   [7] BiocParallel_1.26.0           ggplot2_3.3.3
    ##   [9] digest_0.6.27                 ensembldb_2.16.0
    ##  [11] htmltools_0.5.1.1             magick_2.7.2
    ##  [13] fansi_0.4.2                   magrittr_2.0.1
    ##  [15] checkmate_2.0.0               memoise_2.0.0
    ##  [17] BSgenome_1.60.0               cluster_2.1.2
    ##  [19] ROCR_1.0-11                   matrixStats_0.58.0
    ##  [21] prettyunits_1.1.1             jpeg_0.1-8.1
    ##  [23] colorspace_2.0-1              blob_1.2.1
    ##  [25] rappdirs_0.3.3                xfun_0.23
    ##  [27] dplyr_1.0.6                   crayon_1.4.1
    ##  [29] RCurl_1.98-1.3                jsonlite_1.7.2
    ##  [31] graph_1.70.0                  survival_3.2-11
    ##  [33] VariantAnnotation_1.38.0      glue_1.4.2
    ##  [35] gtable_0.3.0                  zlibbioc_1.38.0
    ##  [37] DelayedArray_0.18.0           scales_1.1.1
    ##  [39] futile.options_1.0.1          DBI_1.1.1
    ##  [41] Rcpp_1.0.6                    xtable_1.8-4
    ##  [43] progress_1.2.2                htmlTable_2.2.1
    ##  [45] foreign_0.8-81                bit_4.0.4
    ##  [47] OrganismDbi_1.34.0            Formula_1.2-4
    ##  [49] AnnotationForge_1.34.0        getopt_1.20.3
    ##  [51] htmlwidgets_1.5.3             httr_1.4.2
    ##  [53] RColorBrewer_1.1-2            ellipsis_0.3.2
    ##  [55] farver_2.1.0                  pkgconfig_2.0.3
    ##  [57] XML_3.99-0.6                  Gviz_1.36.0
    ##  [59] nnet_7.3-16                   sass_0.4.0
    ##  [61] utf8_1.2.1                    labeling_0.4.2
    ##  [63] tidyselect_1.1.1              rlang_0.4.11
    ##  [65] reshape2_1.4.4                later_1.2.0
    ##  [67] AnnotationDbi_1.54.0          biocViews_1.60.0
    ##  [69] munsell_0.5.0                 BiocVersion_3.13.1
    ##  [71] tools_4.1.0                   cachem_1.0.5
    ##  [73] generics_0.1.0                RSQLite_2.2.7
    ##  [75] evaluate_0.14                 stringr_1.4.0
    ##  [77] fastmap_1.1.0                 yaml_2.2.1
    ##  [79] knitr_1.33                    bit64_4.0.5
    ##  [81] purrr_0.3.4                   KEGGREST_1.32.0
    ##  [83] AnnotationFilter_1.16.0       RBGL_1.68.0
    ##  [85] mime_0.10                     formatR_1.9
    ##  [87] biomaRt_2.48.0                compiler_4.1.0
    ##  [89] rstudioapi_0.13               filelock_1.0.2
    ##  [91] curl_4.3.1                    png_0.1-7
    ##  [93] interactiveDisplayBase_1.30.0 tibble_3.1.2
    ##  [95] bslib_0.2.5.1                 stringi_1.6.2
    ##  [97] highr_0.9                     GenomicFeatures_1.44.0
    ##  [99] lattice_0.20-44               ProtGenerics_1.24.0
    ## [101] Matrix_1.3-3                  vctrs_0.3.8
    ## [103] stringdist_0.9.6.3            BiocCheck_1.28.0
    ## [105] pillar_1.6.1                  lifecycle_1.0.0
    ## [107] optparse_1.6.6                RUnit_0.4.32
    ## [109] BiocManager_1.30.15           jquerylib_0.1.4
    ## [111] data.table_1.14.0             bitops_1.0-7
    ## [113] httpuv_1.6.1                  colorRamps_2.3
    ## [115] R6_2.5.0                      BiocIO_1.2.0
    ## [117] latticeExtra_0.6-29           bookdown_0.22
    ## [119] promises_1.2.0.1              gridExtra_2.3
    ## [121] codetools_0.2-18              lambda.r_1.2.4
    ## [123] dichromat_2.0-0               assertthat_0.2.1
    ## [125] SummarizedExperiment_1.22.0   rjson_0.2.20
    ## [127] withr_2.4.2                   GenomicAlignments_1.28.0
    ## [129] GenomeInfoDbData_1.2.6        hms_1.1.0
    ## [131] grid_4.1.0                    rpart_4.1-15
    ## [133] rmarkdown_2.8                 MatrixGenerics_1.4.0
    ## [135] biovizBase_1.40.0             Biobase_2.52.0
    ## [137] shiny_1.6.0                   base64enc_0.1-3
    ## [139] restfulr_0.0.13

# References

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.