1 Introduction

Among the various post-transcriptional RNA modifications, 2’-O methylations are commonly found in rRNA and tRNA. They promote the endo conformation of the ribose and confere resistance to alkaline degradation by preventing a nucleophilic attack on the 3’-phosphate especially in flexible RNA, which is fascilitated by high pH conditions. This property can be queried using a method called RiboMethSeq [Birkedal et al. (2015)] for which RNA is treated in alkaline conditions and RNA fragments are used to prepare a sequencing library [Marchand et al. (2017)].

At position containing a 2’-O methylations, read ends are less frequent, which is used to detect and score the 2’-O methylations.

The ModRiboMethSeq class uses the the ProtectedEndSequenceData class to store and aggregate data along the transcripts. The calculated scores follow the nomenclature of [Birkedal et al. (2015);Galvanin et al. (2019)] with the names scoreRMS (default), scoreA, scoreB and scoreMean.

## snapshotDate(): 2019-10-22
library(rtracklayer)
library(GenomicRanges)
library(RNAmodR.RiboMethSeq)
library(RNAmodR.Data)

2 Example workflow

The example workflow is limited to two 2’-O methylated position on 5.8S rRNA, since the size of the raw data is limited. For 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.RMS.gff3())
sequences <- RNAmodR.Data.example.RMS.fasta()
files <- list("Sample1" = c(treated = RNAmodR.Data.example.RMS.1()),
              "Sample2" = c(treated = RNAmodR.Data.example.RMS.2()))

3 Analysis of data

The analysis is triggered by the construction of a ModSetRiboMethSeq 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).

msrms <- ModSetRiboMethSeq(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
msrms
## ModSetRiboMethSeq of length 2
## names(2): Sample1 Sample2
## | Modification type(s):  Am / Cm / Gm / Um                                      
##                        Sample1 Sample2
## | Modifications found:      no yes (1)
## | Settings:
##         minCoverage minReplicate  find.mod maxLength minSignal
##           <integer>    <integer> <logical> <integer> <integer>
## Sample1          10            1      TRUE        50        10
## Sample2          10            1      TRUE        50        10
##         flankingRegion
##              <integer>
## Sample1              6
## Sample2              6
##         minScoreA minScoreB minScoreRMS minScoreMean flankingRegionMean
##         <numeric> <numeric>   <numeric>    <numeric>          <integer>
## Sample1       0.6       3.6        0.75         0.75                  2
## Sample2       0.6       3.6        0.75         0.75                  2
##               weights
##         <NumericList>
## Sample1   0.9,1,0,...
## Sample2   0.9,1,0,...
##         scoreOperator
##           <character>
## Sample1             &
## Sample2             &

4 Visualizing the results

To compare samples, we need to know, which positions should be part of the comparison. This can either be done by aggregating the detect over all samples and use the union or intersect or by using publish data. We want to assemble a GRanges object from the latter by utilising the infomation from the snoRNAdb [Lestrade and Weber (2006)].

In this specific example only information for the 5.8S RNA is used, since the example data would be to big otherwise. The information regarding the parent and seqname must match the information used as the annotation data. Check that it matches the output of ranges() on a SequenceData, Modifier oder ModifierSet object.

table <- read.csv2(RNAmodR.Data.snoRNAdb(), stringsAsFactors = FALSE)
## snapshotDate(): 2019-10-22
## see ?RNAmodR.Data and browseVignettes('RNAmodR.Data') for documentation
## loading from cache
table <- table[table$hgnc_id == "53533",] # Subset to RNA5.8S
# keep only the current coordinates
table <- table[,1:7]
snoRNAdb <- GRanges(seqnames = "chr1",
              ranges = IRanges(start = table$position,
                               width = 1),
              strand = "+",
              type = "RNAMOD",
              mod = table$modification,
              Parent = "1", #this is the transcript id
              Activity = IRanges::CharacterList(strsplit(table$guide,",")))
coord <- split(snoRNAdb,snoRNAdb$Parent)

In addition to the coordinates of published, we also want to include more meaningful names for the transcripts. For this we provide a data.frame with two columns, tx_id and name. All values in the first column have to match transcript IDs.

ranges(msrms)
## GRangesList object of length 1:
## $`1`
## GRanges object with 1 range and 3 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name exon_rank
##          <Rle> <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1     1-157      + |         1 NR_003285.2         1
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
alias <- data.frame(tx_id = "1", name = "5.8S rRNA", stringsAsFactors = FALSE)
plotCompareByCoord(msrms[c(2,1)], coord, alias = alias)
Heatmap showing RiboMethSeq scores for 2'-O methylated positions on the 5.8S rRNA.

Figure 1: Heatmap showing RiboMethSeq scores for 2’-O methylated positions on the 5.8S rRNA

Results can also be compared on a sequence level, by selecting specific coordinates to compare.

singleCoord <- coord[[1]][1,]
plotDataByCoord(msrms, singleCoord)
RiboMethSeq scores around Um(14) on 5.8S rRNA.

Figure 2: RiboMethSeq scores around Um(14) on 5.8S rRNA

By default only the RiboMethSeq score and the ScoreMean are shown. The raw sequence data can be inspected as well

singleCoord <- coord[[1]][1,]
plotDataByCoord(msrms, singleCoord, showSequenceData = TRUE)
RiboMethSeq scores around Um(14) on 5.8S rRNA. Sequence data is shown by setting `showSequenceData = TRUE`.

Figure 3: RiboMethSeq scores around Um(14) on 5.8S rRNA
Sequence data is shown by setting showSequenceData = TRUE.

5 Performance

To access the performance of the method in combination with samples used, use the plotROC function.

plotROC(msrms,coord)
TPR versus FPR plot.

Figure 4: TPR versus FPR plot

The example given here should be regarded as a proof of concept. Based on the results, minimal scores for calling modified positions can be adjusted to the individual requirements.

settings(msrms) <- list(minScoreMean = 0.7)
msrms
## ModSetRiboMethSeq of length 2
## names(2): Sample1 Sample2
## | Modification type(s):  Am / Cm / Gm / Um                                      
##                        Sample1 Sample2
## | Modifications found:      no yes (1)
## | Settings:
##         minCoverage minReplicate  find.mod maxLength minSignal
##           <integer>    <integer> <logical> <integer> <integer>
## Sample1          10            1      TRUE        50        10
## Sample2          10            1      TRUE        50        10
##         flankingRegion
##              <integer>
## Sample1              6
## Sample2              6
##         minScoreA minScoreB minScoreRMS minScoreMean flankingRegionMean
##         <numeric> <numeric>   <numeric>    <numeric>          <integer>
## Sample1       0.6       3.6        0.75          0.7                  2
## Sample2       0.6       3.6        0.75          0.7                  2
##               weights
##         <NumericList>
## Sample1   0.9,1,0,...
## Sample2   0.9,1,0,...
##         scoreOperator
##           <character>
## Sample1             &
## Sample2             &
## Warning: Settings were changed after data aggregation or modification search.
## Rerun with modify(x,force = TRUE) to update with current settings.

As the warning suggested, after modifying the settings the results should be updated by running modify(x,force = TRUE).

msrms2 <- modify(msrms,force = TRUE)

6 Session info

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] Rsamtools_2.2.0           RNAmodR.Data_0.99.6      
##  [3] ExperimentHubData_1.12.0  AnnotationHubData_1.16.0 
##  [5] futile.logger_1.4.3       ExperimentHub_1.12.0     
##  [7] AnnotationHub_2.18.0      BiocFileCache_1.10.0     
##  [9] dbplyr_1.4.2              RNAmodR.RiboMethSeq_1.0.0
## [11] RNAmodR_1.0.0             Modstrings_1.2.0         
## [13] Biostrings_2.54.0         XVector_0.26.0           
## [15] rtracklayer_1.46.0        GenomicRanges_1.38.0     
## [17] GenomeInfoDb_1.22.0       IRanges_2.20.0           
## [19] S4Vectors_0.24.0          BiocGenerics_0.32.0      
## [21] 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] labeling_0.3                  tidyselect_0.2.5             
##  [69] rlang_0.4.1                   reshape2_1.4.3               
##  [71] later_1.0.0                   AnnotationDbi_1.48.0         
##  [73] biocViews_1.54.0              munsell_0.5.0                
##  [75] BiocVersion_3.10.1            tools_3.6.1                  
##  [77] RSQLite_2.1.2                 assertive.reflection_0.0-4   
##  [79] fastmap_1.0.1                 evaluate_0.14                
##  [81] stringr_1.4.0                 yaml_2.2.0                   
##  [83] knitr_1.25                    bit64_0.9-7                  
##  [85] assertive.matrices_0.0-2      caTools_1.17.1.2             
##  [87] purrr_0.3.3                   AnnotationFilter_1.10.0      
##  [89] assertive.sets_0.0-3          RBGL_1.62.0                  
##  [91] mime_0.7                      formatR_1.7                  
##  [93] biomaRt_2.42.0                compiler_3.6.1               
##  [95] rstudioapi_0.10               curl_4.2                     
##  [97] interactiveDisplayBase_1.24.0 tibble_2.1.3                 
##  [99] stringi_1.4.3                 highr_0.8                    
## [101] GenomicFeatures_1.38.0        lattice_0.20-38              
## [103] assertive.base_0.0-7          ProtGenerics_1.18.0          
## [105] Matrix_1.2-17                 assertive.data_0.0-3         
## [107] vctrs_0.2.0                   stringdist_0.9.5.5           
## [109] BiocCheck_1.22.0              pillar_1.4.2                 
## [111] optparse_1.6.4                RUnit_0.4.32                 
## [113] BiocManager_1.30.9            data.table_1.12.6            
## [115] bitops_1.0-6                  httpuv_1.5.2                 
## [117] colorRamps_2.3                assertive.types_0.0-3        
## [119] R6_2.4.0                      latticeExtra_0.6-28          
## [121] bookdown_0.14                 assertive.properties_0.0-4   
## [123] promises_1.1.0                KernSmooth_2.23-16           
## [125] gridExtra_2.3                 codetools_0.2-16             
## [127] lambda.r_1.2.4                dichromat_2.0-0              
## [129] gtools_3.8.1                  assertthat_0.2.1             
## [131] SummarizedExperiment_1.16.0   openssl_1.4.1                
## [133] rBiopaxParser_2.26.0          GenomicAlignments_1.22.0     
## [135] GenomeInfoDbData_1.2.2        hms_0.5.1                    
## [137] grid_3.6.1                    rpart_4.1-15                 
## [139] rmarkdown_1.16                biovizBase_1.34.0            
## [141] Biobase_2.46.0                shiny_1.4.0                  
## [143] base64enc_0.1-3               assertive_0.3-5

References

Birkedal, Ulf, Mikkel Christensen-Dalsgaard, Nicolai Krogh, Radhakrishnan Sabarinathan, Jan Gorodkin, and Henrik Nielsen. 2015. “Profiling of Ribose Methylations in Rna by High-Throughput Sequencing.” Angewandte Chemie (International Ed. In English) 54 (2):451–55. https://doi.org/10.1002/anie.201408362.

Galvanin, Adeline, Lilia Ayadi, Mark Helm, Yuri Motorin, and Virginie Marchand. 2019. “Mapping and Quantification of tRNA 2’-O-Methylation by Ribomethseq.” Edited by Narendra Wajapeyee and Romi Gupta. New York, NY: Springer New York, 273–95. https://doi.org/10.1007/978-1-4939-8808-2_21.

Lestrade, Laurent, and Michel J. Weber. 2006. “snoRNA-LBME-db, a comprehensive database of human H/ACA and C/D box snoRNAs.” Nucleic Acids Research 34 (January):D158–D162. https://doi.org/10.1093/nar/gkj002.

Marchand, Virginie, Lilia Ayadi, Aseel El Hajj, Florence Blanloeil-Oillo, Mark Helm, and Yuri Motorin. 2017. “High-Throughput Mapping of 2’-O-Me Residues in Rna Using Next-Generation Sequencing (Illumina Ribomethseq Protocol).” Methods in Molecular Biology (Clifton, N.J.) 1562:171–87. https://doi.org/10.1007/978-1-4939-6807-7_12.