rRMSAnalyzer: analysis script of RiboMethSeq | Human by Allyson Moureaux inspired by Virginie Marcel the 19/06/2025 thanks to rRMSAnalyzer & rRMSReports packages & rRMSFonctions
For more specific help: https://ribosomecrcl.github.io/rRMSAnalyzer/
Before starting analysis
Create a new Project and prepare working environment
go to
/file/new project/ create a new folder
In this folder, copy and
paste files of directory “rRMSAnalyser_template”:
metadata.csv +
rRMSAnalyzer_template.R
For Human analysis take
“human.methylated.rda” and “human.suspected.rda”
Don’t forget to put
your 5’end counts file in this directory.
open metadata.csv in
Excel, fill in the file and save it under .csv with “;” as
separator.
Package update RMSAnalyzer and RMSReport
library(devtools)
devtools::install_github("RibosomeCRCL/rRMSAnalyzer")
library(rRMSAnalyzer)
library(rRMSReports)RiboClass creation
Model data to get used to the package
data("ribo_toy")Personal data preparation
ribo <- load_ribodata(
#data & metadata files path
count_path = "~/Desktop/RMS/path/counts/",
metadata = "~/Desktop/RMS/path/metadata.csv",# data & metadata files separator
count_sep = "\t",
metadata_sep = ";",
# count data parameters :
count_header = FALSE,
count_value = 3,
count_rnaid = 1,
count_pos = 2,
# Metadata parameters :
metadata_key = "filename",
metadata_id = "samplename",
# c-score parameters :
flanking = 6,
method = "median",
ncores = 1)QC step
QC_report for all data
report_qc(ribo = ribo, specie = "human", library_col = "run", project_name = "my_project") #ribo_toy lib_col = "run"adjustment ComBat-Seq if necessary
plot_pca(ribo, "lib") #ribo_toy lib_col = "run"
ribo_adj <- adjust_bias(ribo_toy,"run")
plot_pca(ribo_adj,"lib") #ribo_toy lib_col = "run" QC_report for adjusted data
report_qc(ribo = ribo_adj, specie = "human", library_col = "run", project_name = "name", comments = "./comment_QC.Rmd") #ribo_toy lib_col = "run"Remove problematic or not used samples (optionnal)
Keep a small number of samples
ribo_adj_small <- keep_ribo_samples(ribo_adj,c("sample1","sample2","sample3","sample4"))Remove a small number of samples
ribo_adj_small <- remove_ribo_samples(ribo_adj,c("xxxx","xxxccc","RNA_ref-1","RNA_ref-2"))Uniformisation of Riboclass name
ribo_adj <- ribo_adj_small Human rRNA 2’Ome sites annotation
You must have download the file on git and saved it locally
For known sites ….
Choose between this section or the next one. If both of known and suspected sites analysis are desired, script has to be run twice from here.
data("human_methylated")
cat("human_methylated's rna names: ", unique(human_methylated$rRNA),"\n")
cat("ribo's rna names: ", as.character(ribo$rna_names$current_name)) # change name of RiboClass if necessary
ribo_adj_name <- rename_rna(ribo, #change name of RiboClass if necessary for ribo_adj if you adjusted the data
new_names = c("5S", "5.8S", "18S", "28S"))
ribo_adj_annot <- annotate_site(ribo_adj_name,
annot = human_methylated,
anno_rna = "rRNA",
anno_pos = "Position",
anno_value = "Nomenclature")…. or for suspected sites
If you choose to analyse known site first pass this chunk and come back to it when first analysis is done.
redo all the analysis for adjusted data if necessary
load(file = "/home/bioinfo/Desktop/RMS/human.suspected.rda")
ribo_adj_annot <- annotate_site(ribo_adj_annot, human.suspected) # verify rRNA names
ribo_adj_annot <- annotate_site(ribo_adj_annot,
annot = human.suspected,
anno_rna = "rRNA",
anno_pos = "Position",
anno_value = "Nomenclature")Reports
In this section we run the two reports for analysis.
To launch as many time as the number of comparisons asked(A, B, C …).
Ex of metadata:
| filename | samplename | variableA | variableB | variableC | treatment | condition | lib |
|---|---|---|---|---|---|---|---|
| S26_R1.5_counts.csv | variableB_N2_P2 | NA | P2 | NA | variableB | P2 | L2 |
| S27_R1.5_counts.csv | variableA_N2_P1 | P1 | NA | NA | variableA | P1 | L2 |
| S28_R1.5_counts.csv | variableC_N1_ctrl | NA | NA | ctrl | variableC | ctrl | L2 |
| S30_R1.5_counts.csv | variableB_N1_P2 | NA | P2 | NA | variableB | P2 | L2 |
| S31_R1.5_counts.csv | variableC_N2_P2 | NA | NA | P2 | variableC | P2 | L2 |
| S32_R1.5_counts.csv | variableA_N3_P1 | P1 | NA | NA | variableA | P1 | L2 |
| S33_R1.5_counts.csv | variableA_N3_ctrl | ctrl | NA | NA | variableA | ctrl | L2 |
| S34_R1.5_counts.csv | variableC_N1_ctrl | NA | NA | ctrl | variableC | ctrl | L2 |
| S35_R1.5_counts.csv | variableA_N1_P2 | P2 | NA | NA | variableA | P2 | L2 |
| S36_R1.5_counts.csv | variableA_N2_ctrl | ctrl | NA | NA | variableA | ctrl | L2 |
| S37_R1.5_counts.csv | variableC_N2_ctrl | NA | NA | ctrl | variableC | ctrl | L2 |
| S38_R1.5_counts.csv | variableA_N2_P2 | P2 | NA | NA | variableA | P2 | L2 |
| S39_R1.5_counts.csv | variableC_N1_P2 | NA | NA | P2 | variableC | P2 | L2 |
| S40_R1.5_counts.csv | variableB_N1_ctrl | NA | ctrl | NA | variableB | ctrl | L2 |
| S41_R1.5_counts.csv | variableA_N1_P1 | P1 | NA | NA | variableA | P1 | L2 |
| S42_R1.5_counts.csv | variableA_N1_ctrl | ctrl | NA | NA | variableA | ctrl | L2 |
| S43_R1.5_counts.csv | variableC_N3_P2 | NA | NA | P2 | variableC | P2 | L2 |
| S44_R1.5_counts.csv | variableB_N3_ctrl | NA | ctrl | NA | variableB | ctrl | L2 |
| S45_R1.5_counts.csv | variableA_N3_P2 | P2 | NA | NA | variableA | P2 | L2 |
If only one column condition, it is still the same thing. Here in the metadata example we can see 3 columns variable. They are the different comparison that is wanted to be analysed. You can replace the name for your custom ones. You can also run the analysis on “condition” column. In this case analysis will be done camparing as many modalities as you have in the same report. (The column treatment is optionnal)
————————————–variableA——————————–
filter data
(to do only if you want to compare a subset of your samples and not all the samples together otherwise pass this chunk)
Data are filtered with metadata to keep only samples present in the column “variableA” (in the case you have some NA in the column). Samples with NA will be erased from subset for the analysis.
kept_samples <- ribo_adj_annot$metadata %>%
dplyr::filter(!is.na(variableA)) %>% # keep lines that are not "NA"
dplyr::pull(samplename)
ribo_adj_annot_variableA <- keep_ribo_samples(ribo_adj_annot,kept_samples) # create a new riboclass name for the subdata variableAcreate mandatory table for the diff_sites report
The last report needs this table to know what comparisons has to be run. Here in variable A we have 3 modalities : P1, P2 and ctrl. We want to compare P1 to the ctrl and P2 to the ctrl. The order in ctrl and cases lines is important. The first argument of the “ctrl” line will be compared to the first argument of the “cases” line. Same for the second third etc. You can add as much argument as you want in each lines as long as they exist in your metadata column (here variableA) and the line ctrl and cases have the same length. For example this :
comparisons <- tibble::tibble(
comp = c(“comp1”, “comp2”),
ctrl = c(“ctrl”),
cases = c(“P1”, “P2”)
)
can’t work because
ctrl has 1 argument and cases has 2.
This :
comparisons <-
tibble::tibble(
comp = c(“comp1”, “comp2”),
ctrl = c(“ctrl”,
“ctrl”),
cases = c(“P1”, “P3”)
)
can’t work because P3
doesn’t exist in variableA column in the metadata.
global profile analysis of rRNA 2’Ome
report_2ome_sites(ribo = ribo_adj_annot_variableA, specie = "human", condition_col = "variableA", project_name = "name", comments = "./comment_2ome_variableA.Rmd") Analysis of rRNAs 2’Ome level differences
Here we need the comparison table created previously.
report_diff_sites(ribo = ribo_adj_annot_variableA, specie = "human", condition_col = "variableA", project_name = "name", comparisons = comparisons, comments = "./comment_diff_site_variableA.Rmd") ————————————–B—————————————-
Now we do the exact same thing but for the variableB column.
filter with metadata to keep only samples present in the colum “variableB”
kept_samples <- ribo_adj_annot$metadata %>%
dplyr::filter(!is.na(variableB)) %>%
dplyr::pull(samplename)
ribo_adj_annot_variableB <- keep_ribo_samples(ribo_adj_annot,kept_samples)create the necessary comparison table for the diff_sites report
2’Ome of rRNAs analysis of the global profile
report_2ome_sites(ribo = ribo_adj_annot_variableB, specie = "human", condition_col = "variableB", project_name = "name", comments = "./comment_2ome_variableB.Rmd") Analysis of rRNAs 2’Ome level differences
report_diff_sites(ribo = ribo_adj_annot_variableB, specie = "human", condition_col = "variableB", project_name = "name", comparisons = comparisons, comments = "./comment_diff_site_variableB.Rmd") ————————————–C—————————————-
Do it again for the last column variableC
filter with metadata to keep only samples present in the colum “variableC”
kept_samples <- ribo_adj_annot$metadata %>%
dplyr::filter(!is.na(variableC)) %>%
dplyr::pull(samplename)
ribo_adj_annot_variableC <- keep_ribo_samples(ribo_adj_annot,kept_samples)create the necessary comparison table for the diff_sites report
Check local environments of sites of interest (optionnal)
If you want to verify a site that shows up as significant in the last report
For all samples
plot_counts_env(ribo_adj_annot, rna= "18S", pos = 1442, samples = "all", condition = "condition")For some samples
plot_counts_env(ribo_adj_annot, rna = "18S", pos = 1442, samples = c("sample1", "sample2"))data exporting
ribo_df <- extract_data(ribo_adj_annot,
col = "cscore",
only_annotated = TRUE)
write.csv(ribo_df, "/path/to/folder/cscore_df.csv") #give pathSession info
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> loaded via a namespace (and not attached):
#> [1] digest_0.6.37 desc_1.4.3 R6_2.6.1 fastmap_1.2.0
#> [5] xfun_0.53 cachem_1.1.0 knitr_1.50 htmltools_0.5.8.1
#> [9] rmarkdown_2.29 lifecycle_1.0.4 cli_3.6.5 sass_0.4.10
#> [13] pkgdown_2.1.3 textshaping_1.0.3 jquerylib_0.1.4 systemfonts_1.2.3
#> [17] compiler_4.5.1 tools_4.5.1 ragg_1.5.0 evaluate_1.0.5
#> [21] bslib_0.9.0 yaml_2.3.10 jsonlite_2.0.0 rlang_1.1.6
#> [25] fs_1.6.6 htmlwidgets_1.6.4