RNAmodR: detecting patterns of post-transcriptional modifications using machine learning
2019-10-29
Package
RNAmodR 1.0.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(): 2019-10-22library(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:
ModifierMLrangerfor models generated with therangerpackage [Wright and Ziegler (2017)]ModifierMLkerasfor models generated with thekeraspackage [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 in .local(x, ...): 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: 0x558ee3248f58>
## <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) <- mlmodelIn 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 35:
## $`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.927066666666667 1
## [2] 0.812133333333333 1
## -------
## seqinfo: 60 sequences from an unspecified genome; no seqlengths
##
## $`3`
## GRanges object with 2 ranges and 5 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] RDN18-1 229 + | D RNAmodR.ML RNAMOD
## [2] RDN18-1 1669 + | D RNAmodR.ML RNAMOD
## score Parent
## <numeric> <character>
## [1] 0.815866666666667 3
## [2] 0.83 3
## -------
## seqinfo: 60 sequences from an unspecified genome; no seqlengths
##
## $`4`
## GRanges object with 5 ranges and 5 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] RDN25-1 2 + | D RNAmodR.ML RNAMOD
## [2] RDN25-1 748 + | D RNAmodR.ML RNAMOD
## [3] RDN25-1 1720 + | D RNAmodR.ML RNAMOD
## [4] RDN25-1 2297 + | D RNAmodR.ML RNAMOD
## [5] RDN25-1 2571 + | D RNAmodR.ML RNAMOD
## score Parent
## <numeric> <character>
## [1] 0.886333333333334 4
## [2] 0.804266666666667 4
## [3] 0.8189 4
## [4] 0.8014 4
## [5] 0.8118 4
## -------
## seqinfo: 60 sequences from an unspecified genome; no seqlengths
##
## ...
## <32 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 44 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
## ... ... ... ... . ... ... ...
## [40] tV-CAC-D 47 + | D RNAmodR.ML RNAMOD
## [41] tW-CCA-G1 16 + | D RNAmodR.ML RNAMOD
## [42] tW-CCA-G1 46 + | D RNAmodR.ML RNAMOD
## [43] tY-GTA-D 21 + | D RNAmodR.ML RNAMOD
## [44] tY-GTA-D 22 + | D RNAmodR.ML RNAMOD
## score Parent
## <numeric> <character>
## [1] 0.960133333333333 6
## [2] 0.953575 7
## [3] 0.927171428571428 8
## [4] 0.93705 8
## [5] 0.833783333333333 11
## ... ... ...
## [40] 0.931816666666666 57
## [41] 0.803391666666666 59
## [42] 0.825916666666667 59
## [43] 0.902133333333333 60
## [44] 0.8046 60
## -------
## seqinfo: 60 sequences from an unspecified genome; no seqlengthsIn 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.
- 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. - 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
cbindon the data from theSequenceDataobjects (cbindworks onSequenceDataobjects, if they are of the same type. Convert each of them to aCompressedSplitDataFrameListusingas(x,"CompressedSplitDataFrameList")). - 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 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-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] ranger_0.11.2 Rsamtools_2.2.0
## [3] RNAmodR.Data_0.99.6 ExperimentHubData_1.12.0
## [5] AnnotationHubData_1.16.0 futile.logger_1.4.3
## [7] ExperimentHub_1.12.0 AnnotationHub_2.18.0
## [9] BiocFileCache_1.10.0 dbplyr_1.4.2
## [11] RNAmodR.ML_1.0.0 RNAmodR_1.0.0
## [13] Modstrings_1.2.0 Biostrings_2.54.0
## [15] XVector_0.26.0 rtracklayer_1.46.0
## [17] GenomicRanges_1.38.0 GenomeInfoDb_1.22.0
## [19] IRanges_2.20.0 S4Vectors_0.24.0
## [21] BiocGenerics_0.32.0 BiocStyle_2.14.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.5 Hmisc_4.2-0
## [3] plyr_1.8.4 assertive.files_0.0-2
## [5] lazyeval_0.2.2 splines_3.6.1
## [7] BiocParallel_1.20.0 ggplot2_3.2.1
## [9] digest_0.6.22 ensembldb_2.10.0
## [11] htmltools_0.4.0 gdata_2.18.0
## [13] magrittr_1.5 checkmate_1.9.4
## [15] memoise_1.1.0 BSgenome_1.54.0
## [17] assertive.datetimes_0.0-2 assertive.numbers_0.0-2
## [19] cluster_2.1.0 ROCR_1.0-7
## [21] matrixStats_0.55.0 askpass_1.1
## [23] prettyunits_1.0.2 colorspace_1.4-1
## [25] blob_1.2.0 rappdirs_0.3.1
## [27] assertive.strings_0.0-3 xfun_0.10
## [29] dplyr_0.8.3 jsonlite_1.6
## [31] crayon_1.3.4 RCurl_1.95-4.12
## [33] graph_1.64.0 zeallot_0.1.0
## [35] survival_2.44-1.1 VariantAnnotation_1.32.0
## [37] glue_1.3.1 gtable_0.3.0
## [39] zlibbioc_1.32.0 DelayedArray_0.12.0
## [41] scales_1.0.0 futile.options_1.0.1
## [43] DBI_1.0.0 assertive.data.uk_0.0-2
## [45] assertive.models_0.0-2 Rcpp_1.0.2
## [47] assertive.code_0.0-3 xtable_1.8-4
## [49] progress_1.2.2 htmlTable_1.13.2
## [51] foreign_0.8-72 bit_1.1-14
## [53] OrganismDbi_1.28.0 assertive.data.us_0.0-2
## [55] Formula_1.2-3 AnnotationForge_1.28.0
## [57] getopt_1.20.3 htmlwidgets_1.5.1
## [59] httr_1.4.1 gplots_3.0.1.1
## [61] RColorBrewer_1.1-2 acepack_1.4.1
## [63] pkgconfig_2.0.3 XML_3.98-1.20
## [65] Gviz_1.30.0 nnet_7.3-12
## [67] tidyselect_0.2.5 rlang_0.4.1
## [69] reshape2_1.4.3 later_1.0.0
## [71] AnnotationDbi_1.48.0 biocViews_1.54.0
## [73] munsell_0.5.0 BiocVersion_3.10.1
## [75] tools_3.6.1 RSQLite_2.1.2
## [77] assertive.reflection_0.0-4 fastmap_1.0.1
## [79] evaluate_0.14 stringr_1.4.0
## [81] yaml_2.2.0 knitr_1.25
## [83] bit64_0.9-7 assertive.matrices_0.0-2
## [85] caTools_1.17.1.2 purrr_0.3.3
## [87] AnnotationFilter_1.10.0 assertive.sets_0.0-3
## [89] RBGL_1.62.0 mime_0.7
## [91] formatR_1.7 biomaRt_2.42.0
## [93] compiler_3.6.1 rstudioapi_0.10
## [95] curl_4.2 interactiveDisplayBase_1.24.0
## [97] tibble_2.1.3 stringi_1.4.3
## [99] highr_0.8 GenomicFeatures_1.38.0
## [101] lattice_0.20-38 assertive.base_0.0-7
## [103] ProtGenerics_1.18.0 Matrix_1.2-17
## [105] assertive.data_0.0-3 vctrs_0.2.0
## [107] stringdist_0.9.5.5 BiocCheck_1.22.0
## [109] pillar_1.4.2 optparse_1.6.4
## [111] RUnit_0.4.32 BiocManager_1.30.9
## [113] data.table_1.12.6 bitops_1.0-6
## [115] httpuv_1.5.2 colorRamps_2.3
## [117] assertive.types_0.0-3 R6_2.4.0
## [119] latticeExtra_0.6-28 bookdown_0.14
## [121] assertive.properties_0.0-4 promises_1.1.0
## [123] KernSmooth_2.23-16 gridExtra_2.3
## [125] codetools_0.2-16 lambda.r_1.2.4
## [127] dichromat_2.0-0 gtools_3.8.1
## [129] assertthat_0.2.1 SummarizedExperiment_1.16.0
## [131] openssl_1.4.1 rBiopaxParser_2.26.0
## [133] GenomicAlignments_1.22.0 GenomeInfoDbData_1.2.2
## [135] hms_0.5.1 grid_3.6.1
## [137] rpart_4.1-15 rmarkdown_1.16
## [139] biovizBase_1.34.0 Biobase_2.46.0
## [141] shiny_1.4.0 base64enc_0.1-3
## [143] assertive_0.3-5References
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.