# RNAmodR.ML: detecting patterns of post-transcriptional modifications using machine learning

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

#### 2020-10-27

#### Package

RNAmodR.ML 1.4.0

# 1 Introduction

Post-transcriptional modifications can be found abundantly in rRNA and tRNA and can be detected classically via several strategies. However, difficulties arise if the identity and the position of the modified nucleotides is to be determined at the same time. Classically, a primer extension, a form of reverse transcription (RT), would allow certain modifications to be accessed by blocks during the RT changes or changes in the cDNA sequences. Other modification would need to be selectively treated by chemical reactions to influence the outcome of the reverse transcription.

With the increased availability of high throughput sequencing, these classical methods were adapted to high throughput methods allowing more RNA molecules to be accessed at the same time. However, patterns of some modification cannot be detected by accessing small number of parameters. For these cases machine learning models can be trained on data from positions known to be modified in order to detect additional modified positions.

To extend the functionality of the `RNAmodR` package and classical detection strategies used for RiboMethSeq or AlkAnilineSeq (see `RNAmodR.RiboMethSeq` and `RNAmodR.AlkAnilineSeq` packages) towards detection through machine learning models, `RNAmodR.ML` provides classes and an example workflow.

# 2 Using RNAmodR.ML

    ## snapshotDate(): 2020-10-02

    library(rtracklayer)
    library(GenomicRanges)
    library(RNAmodR.ML)
    library(RNAmodR.Data)

The `ModifierML` class extends the `Modifier` class from the `RNAmodR` package and adds one slot, `mlModel`, a getter/setter `getMLModel`/`setMLModel`, an additional `useMLModel` function to be called from the `aggregate` function.

The slot `mlModel` can either be an empty character or contain a name of a `ModifierMLModel` class, which is loaded upon creation of a `ModifierML` object, and serves as a wrapper around a machine learning model. For different types of machine learning models `ModifierMLModel` derived classes are available, which currently are:

- `ModifierMLranger` for models generated with the `ranger` package [Wright and Ziegler (2017)]
- `ModifierMLkeras` for models generated with the `keras` package [Allaire and Chollet (2018)]

To illustrate how to develop a machine learning model with help from the `RNAmodR.ML` package, an example is given below.

# 3 Development of new `Modifier` class

As an example for this vignette, we will try to detect D positions in AlkAnilineSeq data. First define a specific `ModifierML` class loading pileup and coverage data. In this example, the RNA specific `RNAModifierML` class is used.

    setClass("ModMLExample",
             contains = c("RNAModifierML"),
             prototype = list(mod = c("D"),
                              score = "score",
                              dataType = c("PileupSequenceData",
                                           "CoverageSequenceData"),
                              mlModel = character(0)))
    # constructor function for ModMLExample
    ModMLExample <- function(x, annotation = NA, sequences = NA, seqinfo = NA,
                               ...){
      RNAmodR:::Modifier("ModMLExample", x = x, annotation = annotation,
                         sequences = sequences, seqinfo = seqinfo, ...)
    }

    setClass("ModSetMLExample",
             contains = "ModifierSet",
             prototype = list(elementType = "ModMLExample"))
    # constructor function for ModSetMLExample
    ModSetMLExample <- function(x, annotation = NA, sequences = NA, seqinfo = NA,
                                  ...){
      RNAmodR:::ModifierSet("ModMLExample", x, annotation = annotation,
                            sequences = sequences, seqinfo = seqinfo, ...)
    }

Since the `mlModel` slot contains an empty character, the creation of the object will not automatically trigger a search for modifications. However, it will aggregate the data in the format we want to use. The `aggregate_example` function is just an example and the aggregation of the data is part of the model building. (See (#Summary))

    setMethod(
      f = "aggregateData",
      signature = signature(x = "ModMLExample"),
      definition =
        function(x){
          aggregate_example(x)
        }
    )

# 4 Getting training data

To gather training data, we just create a `ModMLExample` object and let it do its job.

    me <-  ModMLExample(files[[1]], annotation, sequences)

    ## Import genomic features from the file as a GRanges object ... OK
    ## Prepare the 'metadata' data frame ... OK
    ## Make the TxDb object ... OK
    ## Loading Pileup data from BAM files ... OK
    ## Loading Coverage data from BAM files ... OK
    ## Aggregating data and calculating scores ... Starting to search for 'Dihydrouridine' ...

    ## Warning: ML model not set. Skipped search for modifications ...

    ## done.

Afterwards we need to load/create coordinates for positions known to be modified as well as positions known to be unmodified.

    data("dmod",package = "RNAmodR.ML")
    # we just select the next U position from known positions
    nextUPos <- function(gr){
      nextU <- lapply(seq.int(1L,2L),
                      function(i){
                        subseq <- subseq(sequences(me)[dmod$Parent], start(dmod)+3L)
                        pos <- start(dmod) + 2L +
                          vapply(strsplit(as.character(subseq),""),
                        function(y){which(y == "U")[i]},integer(1))
                        ans <- dmod
                        ranges(ans) <- IRanges(start = pos, width = 1L)
                        ans
                      })
      nextU <- do.call(c,nextU)
      nextU$mod <- NULL
      unique(nextU)
    }
    nondmod <- nextUPos(dmod)
    nondmod <- nondmod[!(nondmod %in% dmod)]
    coord <- unique(c(dmod,nondmod))
    coord <- coord[order(as.integer(coord$Parent))]

With these coordinates the aggregated data of the `ModMLExample` can be subset to a training data set using the function `trainingData`.

    trainingData <- trainingData(me,coord)
    trainingData <- unlist(trainingData, use.names = FALSE)
    # converting logical labels to integer
    trainingData$labels <- as.integer(trainingData$labels)

# 5 Training a model

How a specific model can be trained or what kind of strategies can be employed to successfully train a model, is out of scope for the vignette. For this example a random forest is trained using the functionality provided by the `ranger` package.

    library(ranger)
    model <- ranger(labels ~ ., data.frame(trainingData))

# 6 Constructing a ‘ModifierMLModel’

Now, the trained model can be used to create a new `ModifierMLModel` class and object.

    setClass("ModifierMLexample",
             contains = c("ModifierMLranger"),
             prototype = list(model = model))
    ModifierMLexample <- function(...){
      new("ModifierMLexample")
    }
    mlmodel <- ModifierMLexample()

To be able to use the model via the `ModifierMLModel` class, we also need to define an accessor to the predictions made by the model. This function is called `useModel` and is already prefined for the `ModifierMLModel` classes mentioned above in secion Using RNAmodR.ML.

    getMethod("useModel", c("ModifierMLranger","ModifierML"))

    ## Method Definition:
    ##
    ## function (x, y)
    ## {
    ##     data <- getAggregateData(y)
    ##     model <- x@model
    ##     if (!is(model, "ranger")) {
    ##         stop("model is not a ranger model")
    ##     }
    ##     unlisted_data <- unlist(data, use.names = FALSE)
    ##     p <- predict(x@model, data.frame(unlisted_data))
    ##     relist(p$predictions, data)
    ## }
    ## <bytecode: 0x000000003da3c588>
    ## <environment: namespace:RNAmodR.ML>
    ##
    ## Signatures:
    ##         x                  y
    ## target  "ModifierMLranger" "ModifierML"
    ## defined "ModifierMLranger" "ModifierML"

If the results of these function is not usable for a specific model, it can redefined for your needs. As defined by `RNAmodR.ML` the function returns a `NumericList` along the aggregated data of the `ModifierML` object.

# 7 Setting and using the model

The generated `ModifierMLexample` wrapper can now be set for the `ModifierML` object using the `setMLModel` function. (If a model is saved as part of a package, this step ist not necessary, because it can be part of the class definition)

    setMLModel(me) <- mlmodel

In order for the prediction data to be usable, another function has to be implemented to save the predictions in the aggregated data. The function is called `useMLModel`.

    setMethod(f = "useMLModel",
              signature = signature(x = "ModMLExample"),
              definition =
                function(x){
                  predictions <- useModel(getMLModel(x), x)
                  data <- getAggregateData(x)
                  unlisted_data <- unlist(data, use.names = FALSE)
                  unlisted_data$score <- unlist(predictions)
                  x@aggregate <- relist(unlisted_data,data)
                  x
                }
    )

By re-running the `aggregate` function and force an update of the data, the predictions are made and used to populate the `score` column as outlined above.

    me <- aggregate(me, force = TRUE)

# 8 Performance

During the model building phase some form of a performance measurement usually is included. In addition to these model specific measurements, `RNAmodR.ML` includes the functionality from the `ROCR` package inherited from the `RNAmodR` package. With this the performance of a model can evaluted over the training set or any coordinates.

    plotROC(me, dmod)

Figure 1: Performance of the maching learning model to distinguish unmodified from modified nucleotides

# 9 Using a `ModifierML` class

Since we want to use the `ModifierML` object to detect modifications, we also need to define the `findMod` function. Based on the information on the performance, we set a threshold of `0.8` for the prediction score for detecting D modifications. In the example below this threshold is hardcoded in the `find_mod_example` function, but can also be implemented using the `settings` function.

    setMethod(
      f = "findMod",
      signature = signature(x = "ModMLExample"),
      definition =
        function(x){
          find_mod_example(x, 25L)
        }
    )

Now we can redfine the `ModMLExample` class with the slot `mlModel` already set to the `ModifierMLexample` class. The `ModMLExample` is now complete.

    rm(me)
    setClass("ModMLExample",
             contains = c("RNAModifierML"),
             prototype = list(mod = c("D"),
                              score = "score",
                              dataType = c("PileupSequenceData",
                                           "CoverageSequenceData"),
                              mlModel = "ModifierMLexample"))
    me <-  ModMLExample(files[[1]], annotation, sequences)

    ## Import genomic features from the file as a GRanges object ... OK
    ## Prepare the 'metadata' data frame ... OK
    ## Make the TxDb object ... OK
    ## Loading Pileup data from BAM files ... OK
    ## Loading Coverage data from BAM files ... OK
    ## Aggregating data and calculating scores ... Starting to search for 'Dihydrouridine' ... done.

The detected modifications can be access using the `modifications` function.

    mod <- modifications(me)
    mod <- split(mod, factor(mod$Parent,levels = unique(mod$Parent)))
    mod

    ## GRangesList object of length 36:
    ## $`1`
    ## GRanges object with 2 ranges and 5 metadata columns:
    ##             seqnames    ranges strand |         mod      source        type
    ##                <Rle> <IRanges>  <Rle> | <character> <character> <character>
    ##   [1] Q0020_15S_RRNA        48      + |           D  RNAmodR.ML      RNAMOD
    ##   [2] Q0020_15S_RRNA       323      + |           D  RNAmodR.ML      RNAMOD
    ##           score      Parent
    ##       <numeric> <character>
    ##   [1]  0.918833           1
    ##   [2]  0.809783           1
    ##   -------
    ##   seqinfo: 60 sequences from an unspecified genome; no seqlengths
    ##
    ## $`3`
    ## GRanges object with 1 range and 5 metadata columns:
    ##       seqnames    ranges strand |         mod      source        type     score
    ##          <Rle> <IRanges>  <Rle> | <character> <character> <character> <numeric>
    ##   [1]  RDN18-1       229      + |           D  RNAmodR.ML      RNAMOD    0.8062
    ##            Parent
    ##       <character>
    ##   [1]           3
    ##   -------
    ##   seqinfo: 60 sequences from an unspecified genome; no seqlengths
    ##
    ## $`4`
    ## GRanges object with 4 ranges and 5 metadata columns:
    ##       seqnames    ranges strand |         mod      source        type     score
    ##          <Rle> <IRanges>  <Rle> | <character> <character> <character> <numeric>
    ##   [1]  RDN25-1         2      + |           D  RNAmodR.ML      RNAMOD  0.842100
    ##   [2]  RDN25-1       748      + |           D  RNAmodR.ML      RNAMOD  0.811533
    ##   [3]  RDN25-1      1720      + |           D  RNAmodR.ML      RNAMOD  0.837467
    ##   [4]  RDN25-1      1818      + |           D  RNAmodR.ML      RNAMOD  0.803400
    ##            Parent
    ##       <character>
    ##   [1]           4
    ##   [2]           4
    ##   [3]           4
    ##   [4]           4
    ##   -------
    ##   seqinfo: 60 sequences from an unspecified genome; no seqlengths
    ##
    ## ...
    ## <33 more elements>

# 10 Refining a model

Some of the modification found look reasonable. However, some of the positions seem to be noise.

    options(ucscChromosomeNames=FALSE)

    plotDataByCoord(sequenceData(me),mod[["4"]][1])

Figure 2: Visualization of sequence data

Several options exist to improve the model: Either the threshold applied to the prediction score can be raised to a higher value, like `0.9` or the model can maybe retrained by including these position in another training data set. In addition, the training data might be improved in general by higher sequencing depth.

    nonValidMod <- mod[c("1","4")]
    nonValidMod[["18"]] <- nonValidMod[["18"]][2]
    nonValidMod[["26"]] <- nonValidMod[["26"]][2]
    nonValidMod <- unlist(nonValidMod)
    nonValidMod <- nonValidMod[,"Parent"]
    coord <- unique(c(dmod,nondmod,nonValidMod))
    coord <- coord[order(as.integer(coord$Parent))]

As an example, a new model is trained including the wrongly identified positions from the previous model as unmodified positions.

    trainingData <- trainingData(me,coord)
    trainingData <- unlist(trainingData, use.names = FALSE)
    trainingData$labels <- as.integer(trainingData$labels)

    model2 <- ranger(labels ~ ., data.frame(trainingData), num.trees = 2000)
    setClass("ModifierMLexample2",
             contains = c("ModifierMLranger"),
             prototype = list(model = model2))
    ModifierMLexample2 <- function(...){
      new("ModifierMLexample2")
    }
    mlmodel2 <- ModifierMLexample2()
    me2 <- me
    setMLModel(me2) <- mlmodel2
    me2 <- aggregate(me2, force = TRUE)

After updating the `ModifierMLexample` class and aggregating the data again the performance looks a bit better…

    plotROC(me2, dmod, score="score")

Figure 3: Performance aggregation of multiple samples and strategies

… and leads to a better result for detecting D modifications. Some positions are not detected anymore, which suggest that this model is still not an optimal solution and other factors could and should be improved upon as suggested above.

    setMethod(
      f = "findMod",
      signature = signature(x = "ModMLExample"),
      definition =
        function(x){
          find_mod_example(x, 25L)
        }
    )
    me2 <- modify(me2, force = TRUE)
    modifications(me2)

    ## GRanges object with 45 ranges and 5 metadata columns:
    ##         seqnames    ranges strand |         mod      source        type
    ##            <Rle> <IRanges>  <Rle> | <character> <character> <character>
    ##    [1]  tA-AGC-D        47      + |           D  RNAmodR.ML      RNAMOD
    ##    [2]  tA-TGC-A        47      + |           D  RNAmodR.ML      RNAMOD
    ##    [3]  tC-GCA-B        19      + |           D  RNAmodR.ML      RNAMOD
    ##    [4]  tC-GCA-B        46      + |           D  RNAmodR.ML      RNAMOD
    ##    [5]  tE-CTC-D        20      + |           D  RNAmodR.ML      RNAMOD
    ##    ...       ...       ...    ... .         ...         ...         ...
    ##   [41]  tV-CAC-D        47      + |           D  RNAmodR.ML      RNAMOD
    ##   [42] tW-CCA-G1        16      + |           D  RNAmodR.ML      RNAMOD
    ##   [43] tW-CCA-G1        46      + |           D  RNAmodR.ML      RNAMOD
    ##   [44]  tY-GTA-D        21      + |           D  RNAmodR.ML      RNAMOD
    ##   [45]  tY-GTA-D        22      + |           D  RNAmodR.ML      RNAMOD
    ##            score      Parent
    ##        <numeric> <character>
    ##    [1]  0.970375           6
    ##    [2]  0.956458           7
    ##    [3]  0.910459           8
    ##    [4]  0.930042           8
    ##    [5]  0.818417          11
    ##    ...       ...         ...
    ##   [41]  0.916308          57
    ##   [42]  0.801900          59
    ##   [43]  0.836733          59
    ##   [44]  0.931858          60
    ##   [45]  0.815908          60
    ##   -------
    ##   seqinfo: 60 sequences from an unspecified genome; no seqlengths

In addition to training a single model, several models can be trained and combined to a `ModifierSet`.

    mse <- ModSetMLExample(list(one = me, two = me2))

An overall performance over several models can be analyzed or the individual performance compaired.

    plotROC(mse, dmod, score= "score",
            plot.args = list(avg = "threshold", spread.estimate = "stderror"))

Figure 4: Performance average across models

If several models are trained and each provides useful information, these can be package into a single `ModifierMLModel` class to combine the output of several models. Some of the functions outlined above need, e.g. `useMLModel` and/or `useModel`, to be modified to provide one or more scores for detecting a modification.

# 11 Packaging

If the created models can be saved to file, they can also be included in a package. This would allow others to use the models and the models can potentially be improved upon with increasing version numbers.

# 12 Summary

`RNAmodR.ML` provides the interface for machine learning models to be used with `RNAmodR` to detect modified nucleotides in high throughput sequencing data. For your own work defining a working model might take some time. We hope that by using `RNAmodR.ML` the steps surrounding this crucial step might become a bit easier.

However, if some steps or design choices made for `RNAmodR.ML` do not suit your need, let us know. Contributions are always welcome as well.

# 13 Hints

We also want to provide some additional hints and suggestions for developing machine learning models.

1. the aggregate function is used in the example above as a feature engineering tool. You might want to skip this step, if you want to use a deep learning model for example with `keras`.
2. If you don’t want to engage in a feature enginering step and just want to aggregate the sequence data as is, just do a custom `cbind` on the data from the `SequenceData` objects (`cbind` works on `SequenceData` objects, if they are of the same type. Convert each of them to a `CompressedSplitDataFrameList` using `as(x,"CompressedSplitDataFrameList")`).
3. in a deep learning context, if a coordinate is selected without a flanking value, e.g. when using `trainingData`, a 2D tensor is returned (sample, values). This can be converted into a 3D tensor by providing a flanking value.

# 14 Sessioninfo

    sessionInfo()

    ## R version 4.0.3 (2020-10-10)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## 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] ranger_0.12.1            Rsamtools_2.6.0          RNAmodR.Data_1.3.0
    ##  [4] ExperimentHubData_1.16.0 AnnotationHubData_1.20.0 futile.logger_1.4.3
    ##  [7] ExperimentHub_1.16.0     AnnotationHub_2.22.0     BiocFileCache_1.14.0
    ## [10] dbplyr_1.4.4             RNAmodR.ML_1.4.0         RNAmodR_1.4.0
    ## [13] Modstrings_1.6.0         Biostrings_2.58.0        XVector_0.30.0
    ## [16] rtracklayer_1.50.0       GenomicRanges_1.42.0     GenomeInfoDb_1.26.0
    ## [19] IRanges_2.24.0           S4Vectors_0.28.0         BiocGenerics_0.36.0
    ## [22] 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                    gplots_3.1.0
    ##  [53] RColorBrewer_1.1-2            ellipsis_0.3.1
    ##  [55] pkgconfig_2.0.3               XML_3.99-0.5
    ##  [57] Gviz_1.34.0                   nnet_7.3-14
    ##  [59] tidyselect_1.1.0              rlang_0.4.8
    ##  [61] reshape2_1.4.4                later_1.1.0.1
    ##  [63] AnnotationDbi_1.52.0          munsell_0.5.0
    ##  [65] biocViews_1.58.0              BiocVersion_3.12.0
    ##  [67] tools_4.0.3                   generics_0.0.2
    ##  [69] RSQLite_2.2.1                 evaluate_0.14
    ##  [71] stringr_1.4.0                 fastmap_1.0.1
    ##  [73] yaml_2.2.1                    knitr_1.30
    ##  [75] bit64_4.0.5                   caTools_1.18.0
    ##  [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                KernSmooth_2.23-17
    ## [113] gridExtra_2.3                 codetools_0.2-16
    ## [115] lambda.r_1.2.4                dichromat_2.0-0
    ## [117] gtools_3.8.2                  assertthat_0.2.1
    ## [119] SummarizedExperiment_1.20.0   openssl_1.4.3
    ## [121] rBiopaxParser_2.30.0          GenomicAlignments_1.26.0
    ## [123] GenomeInfoDbData_1.2.4        hms_0.5.3
    ## [125] grid_4.0.3                    rpart_4.1-15
    ## [127] rmarkdown_2.5                 MatrixGenerics_1.2.0
    ## [129] biovizBase_1.38.0             Biobase_2.50.0
    ## [131] shiny_1.5.0                   base64enc_0.1-3

# References

Allaire, JJ, and François Chollet. 2018. Keras: R Interface to ’Keras’. https://CRAN.R-project.org/package=keras.

Wright, Marvin N., and Andreas Ziegler. 2017. “ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R.” Journal of Statistical Software 77 (1):1–17. https://doi.org/10.18637/jss.v077.i01.