1简介

转录后修饰可以在rRNA和tRNA中大量发现，并可以通过多种方法进行经典检测。然而，如果要同时确定修饰核苷酸的身份和位置，就会出现困难。经典地说，引物延伸是逆转录(RT)的一种形式，它允许在逆转录过程中被片段访问某些修饰或cDNA序列的变化。其他修饰需要通过化学反应选择性地处理，以影响逆转录的结果。

随着高通量测序技术的发展，这些经典的测序方法也被应用到高通量测序中，从而可以同时检测到更多的RNA分子。随着这些进展，在mRNA上也检测到转录后修饰。在这些高吞吐量技术中，例如Pseudo-Seq(卡莱尔等，2014), RiboMethSeq(Birkedal et al. 2015)和AlkAnilineSeq(Marchand et al. 2018)每一个都能从用所选方法制备的RNA-Seq数据的脚印中检测到特定类型的修饰。

由于类似的模式可以从这些技术中观察到，生物信息管道的重叠已经是并且将会随着新的新兴测序技术变得更加频繁。

RNAmodR实现类和工作流，以检测高通量测序数据中的转录后RNA修饰。它很容易适应新方法，在最初的方法开发阶段和更复杂的筛选阶段都有帮助。

简单地说,从SequenceData，派生特定的子类来访问对齐读取的特定方面，例如5 '端位置或堆数据。这与一个修饰符类可用于检测单个修改类型的特定模式。的SequenceData类可以由不同的修饰符类，允许容易地适应新方法。

# # snapshotDate (): 2022-04-19

库(rtracklayer)库(Rsamtools)库(genomics icfeatures)库(RNAmodR. data)库(RNAmodR)

annotation <- GFF3File(RNAmodR.Data.example.gff3())序列<- RNAmodR.Data.example.fasta()文件<- c(handled = rnamodr . data .example. bam1 ()， handled = rnamodr . data .example. bam2 ()， handled = rnamodr . data .example. bam3 ())

seqdata <- End5SequenceData(files, annotation = annotation, sequences = sequences)

从文件中导入基因组特征作为GRanges对象…准备“元数据”数据帧…##创建TxDb对象…正在从BAM文件中加载5'端位置数据…好吧

seqdata

## End5SequenceData包含60个元素，包含3个数据列和3个元数据列。1 end5.treated。2 end5.treated。3.##    ## - Seqinfo object with 84 sequences from an unspecified genome; no seqlengths:

names(seqdata) #匹配TxDb对象返回的文本名称colnames(seqdata) #返回所有列名的字符列表bamfiles(seqdata) ranges(seqdata) #从TxDb对象生成序列(seqdata) seqinfo(seqdata)

seqdata [1]

## End5SequenceData包含1个元素，包含3个数据列和3个元数据列。1 end5.treated。2 end5.treated。3.##    ## - Seqinfo object with 84 sequences from an unspecified genome; no seqlengths:

SDF <- seqdata[[1]] SDF

## End5SequenceDataFrame与1649行和3列## end5.处理。1 end5.treated。2 end5.treated。3.##    ## 1 1 4 0 ## 2 0 2 0 ## 3 0 0 0 ## 4 0 0 0 ## 5 0 0 0 ## ... ... ... ... ## 1645 0 0 0 ## 1646 0 0 0 ## 1647 0 0 0 ## 1648 0 0 0 ## 1649 0 0 0 ## ## containing a GRanges object with 1 range and 3 metadata columns: ## seqnames ranges strand | exon_id exon_name exon_rank ##    |    ## [1] Q0020_15S_RRNA 1-1649 + | 1 Q0020 1 ## ------- ## seqinfo: 60 sequences from an unspecified genome; no seqlengths ## ## and a 1649-letter RNAString object ## seq: GUAAAAAAUUUAUAAGAAUAUGAUGUUGGUUCAGAU...UGCGGUGGGCUUAUAAAUAUCUUAAAUAUUCUUACA

Names (sdf) #返回数据范围(sdf)序列的列名(sdf)

自卫队[1:2]

## End5SequenceDataFrame与1649行和2列## end5. treatment。1 end5.treated。2##   ## 1 1 4 ## 2 0 2 ## 3 0 0 ## 4 0 0 ## 5 0 0 ## ... ... ... ## 1645 0 0 ## 1646 0 0 ## 1647 0 0 ## 1648 0 0 ## 1649 0 0 ## ## containing a GRanges object with 1 range and 3 metadata columns: ## seqnames ranges strand | exon_id exon_name exon_rank ##    |    ## [1] Q0020_15S_RRNA 1-1649 + | 1 Q0020 1 ## ------- ## seqinfo: 60 sequences from an unspecified genome; no seqlengths ## ## and a 1649-letter RNAString object ## seq: GUAAAAAAUUUAUAAGAAUAUGAUGUUGGUUCAGAU...UGCGGUGGGCUUAUAAAUAUCUUAAAUAUUCUUACA

自卫队[1:3]

##有3行3列的数据帧## end5. treatment。1 end5.treated。2 end5.treated。3.##    ## 1 1 4 0 ## 2 0 2 0 ## 3 0 0 0

mi <- ModInosine(files, annotation = annotation, sequences = sequences)

从文件中导入基因组特征作为GRanges对象…准备“元数据”数据帧…##创建TxDb对象…正在从BAM文件中加载堆积数据…OK ##聚合数据和计算分数…开始搜索“肌苷”…完成了。

names(mi) #匹配TxDb对象返回的文本名称bamfiles(mi) ranges(mi) #从TxDb对象生成sequences(mi) seqinfo(mi) sequenceData(mi) #返回sequenceData

设置(mi)

# # # # $ minCoverage [1] minReplicate美元10 # # # # # #[1]1 # # # # $。mod ## [1] TRUE ## ## $minScore ## [1] 0.4

设置(mi,“minScore”)

# # 0.4 [1]

设置(mi) <- list(minScore = 0.5)设置(mi，"minScore")

# # 0.5 [1]

序列<- RNAmodR.Data.example.AAS.fasta() annotation <- GFF3File(RNAmodR.Data.example.AAS.gff3()) files <- list("SampleSet1" = c(已处理= RNAmodR.Data.example.wt.1()，已处理= RNAmodR.Data.example.wt.2()，已处理= RNAmodR.Data.example.wt.3())， "SampleSet2" = c(已处理= RNAmodR.Data.example.wt.3())， "SampleSet3" = c(已处理= RNAmodR.Data.example.bud23.1()，已处理= RNAmodR.Data.example.bud23.2())， "SampleSet3" = c(已处理= RNAmodR.Data.example.trm8.2())))

msi <- ModSetInosine(files, annotation = annotation, sequences = sequences)

从文件中导入基因组特征作为GRanges对象…准备“元数据”数据帧…##创建TxDb对象…好吧

名称(msi)

# [1] "SampleSet1" "SampleSet2" "SampleSet3"

msi ([1])

包含11个元素的PileupSequenceData的ModInosine对象。## |输入文件:## - handled: /home/biocbuild/。缓存/R/ExperimentHub/74c624bfa535f_2544 ## -处理:/home/biocbuild/。缓存/R/ExperimentHub/74c625044ae71_2546 ##处理:/home/biocbuild/。缓存/R/ExperimentHub/74c625b673527_2548 ## |核苷酸-修改类型(s): RNA - i# # |修改发现:yes(6) ## |设置:## minCoverage minreplication find。mod minScore ##     ## 10 1 TRUE 0.4

bamfiles(msi) ranges(msi) #生成TxDb对象序列(msi) seqinfo(msi)

2检测到的修改分析

检索找到的修改修饰符或ModifierSet对象通过修改()函数。函数返回一个农庄或GRangesList对象，其中包含与所使用的基因组相关的修改的坐标。例如，如果一个转录本从100号位置开始，在转录本的50号位置包含一个修饰过的核苷酸，则返回的坐标为150。

Mod <-修改(msi) Mod [[1]]

GRanges对象有6个范围和5个元数据列:# # seqnames范围链| mod源类型得分# # < Rle > < IRanges > < Rle > | <人物> <人物> <人物> <数字> # # [1]chr2 34 + |我RNAmodR RNAMOD 0.900932 # # [2] chr4 35 + |我RNAmodR RNAMOD 0.899622 # # [3] chr6 34 + |我RNAmodR RNAMOD 0.984035 # # [4] chr7 67 + |我RNAmodR RNAMOD 0.934553 # # [5] chr9 7 + |我RNAmodR RNAMOD 0.709758 # # [6] chr11 35 + |我RNAmodR RNAMOD父# # <人物> 0.874027 # # # # 2 # # [1][2]4 # # 6 # # [3][4]7 # # 9 # # [5][6]11  ## ------- ## seqinfo:来自未知基因组的11个序列;没有seqlengths

要检索与抄本边界相关的坐标，请使用可选参数perTranscript = TRUE．在这里提供的示例中，这将产生相同的坐标，因为示例的映射使用了一个自定义基因组，它不包含负链和每个转录本染色体上的转录本。

mod <- modified (msi, perTranscript = TRUE) mod[[1]]

GRanges对象有6个范围和5个元数据列:# # seqnames范围链| mod源类型得分# # < Rle > < IRanges > < Rle > | <人物> <人物> <人物> <数字> # # [1]chr2 34 * |我RNAmodR RNAMOD 0.900932 # # [2] chr4 35 * |我RNAmodR RNAMOD 0.899622 # # [3] chr6 34 * |我RNAmodR RNAMOD 0.984035 # # [4] chr7 67 * |我RNAmodR RNAMOD 0.934553 # # [5] chr9 7 * |我RNAmodR RNAMOD 0.709758 # # [6] chr11 35 * |我RNAmodR RNAMOD父# # <人物> 0.874027 # # # # 2 # # [1][2]4 # # 6 # # [3][4]7 # # 9 # # [5][6]11  ## ------- ## seqinfo:来自未知基因组的11个序列;没有seqlengths

2.1康普艾的结果

为了比较样本之间的结果，aModifierSet以及需要比较的位置的定义。为了构造一组位置，我们将使用找到的所有修改的交集作为例子。

mod <- modified (msi) coord <- unique(unlist(mod)) coord$score <- NULL coord$sd <- NULL compareByCoord(msi,coord)

##数据帧6行6列## SampleSet1 SampleSet2 SampleSet3名称位置mod ## <数值> <数值> <数值> <因子> <因子> <字符> ## 1 0.900932 0.998134 0.9536b2 <因子> <因子> 2 34 I ## 2 0.899622 0.856241 0.976928 4 35 I ## 3 0.984035 0.992012 0.993128 6 34 I ## 4 0.934553 0.942905 0.684484 0.681451 97 I ## 5 0.874027 0.971474 0.954782 11 35 I

结果也可以用plotCompareByCoord，它接受一个可选参数别名允许将抄本id转换为其他标识符。对于这一步，构建一个TxDb对象，并使用它来构造修饰符/ModifierSet对象作为注释论点。

txdb < - makeTxDbFromGFF(注释)

从文件中导入基因组特征作为GRanges对象…准备“元数据”数据帧…##创建TxDb对象…好吧

Alias <- data.frame(tx_id = names(id2name(txdb))， name = id2name(txdb)))

plotCompareByCoord(msi, coord, alias = alias)

图1:识别肌苷位置的热图

此外，样本集的顺序可以调整，标准化为任何一个样本集和每个文本显示的位置编号。

plotCompareByCoord(msi[c(3,1,2)]， coord, alias = alias, normalize = "SampleSet3"， perTranscript = TRUE)

图2:识别肌苷位置的热图与标准化分数

计算出的分数和数据可以沿着抄本或抄本块显示出来。使用可选参数showSequenceData除了评分数据之外，序列数据的绘图也可以通过将其设置为来触发真正的．

plotData(msi， "2"， from = 10L, to = 45L, alias = alias) # showSequenceData = FALSE

图3:抄本上的分数包含a到G的转换，表明肌苷的存在

plotData(msi[1:2]， "2"， from = 10L, to = 45L, showSequenceData = TRUE, alias = alias)

图4:抄本上的分数包含a到G的转换，表明肌苷的存在
此图包括详细的堆积序列数据。

3.性能度量

由于从高通量测序数据中检测修饰通常依赖于调用修饰的阈值，因此基于所选分数和可用样本分析方法的性能具有相当大的兴趣。分析其性能、功能plotROC ()的功能实现了，它是ROCR包(Sing et al. 2005)(#)的引用。

对于这个小插图中使用的示例数据，所获得的信息是相当有限的，下面的图应该被视为概念的证明。此外，利用发现的修改站点作为输入plotROC是强烈反对的，因为违背了考试的目的。因此，请将下一大块的这方面也视为概念的证明。

coord plotROC (msi)

图5:TPR vs. FPR图

请您看一下plotROC ?额外的细节。的大部分功能ROCR包可以通过附加参数获得，因此plotROC可以大量定制。

4额外的信息

查看用于分析的读取元数据RNAmodR这个函数统计数据()可以使用。它可以与一堆对象类型一起使用:SequenceData，SequenceDataList，SequenceDataSet，修饰符或ModifierSet．为SequenceData *对象时,BamFile要分析的也必须提供，自动完成修饰符*对象。想了解更多细节，请看?统计数据．

Stats <- Stats (msi)统计信息

## name (3): SampleSet1 SampleSet2 SampleSet3

统计[[" SampleSet1 "]]

## DataFrameList长度为3 ## names(3):已处理已处理

统计[[“SampleSet1”]][["治疗"]]

## DataFrame with 12行6列## seqnames seqlength mapped unmapped used used_distro ##       ## 1 chr1 1800 197050 0 159782 83,1252,860，…## 2 chr2 85 5863 0 2459 2,16,16，…## 3 chr3 76 76905 0 63497 35,478,4106，…## 4 chr4 77 8299 0 6554 6,27,36，…## 5 chr5 74 11758 0 8818 520,105,93，…## ... ... ... ... ... ... ...## 8 chr8 75 144293 0 143068 14,44,48，…## 9 chr9 75 13790 0 9753 1,49,43，…## 10 chr10 85 19861 0 17729 35,21,10，…## 11 chr11 77 10532 0 9086 53,92,185，… ## 12 * 0 0 961095 NA NA

返回的数据统计数据()是一个DataFrame，可以包装成DataFrameList或者一个SimpleList取决于输入类型。必须手动进行数据分析，并可用于生成如下图所示的所分析读取长度分布的输出。

图6:分析中使用的读取的长度分布

5进一步发展

的发展RNAmodR将继续下去。分析工作流程的一般方面将在RNAmodR包，而针对检测转录后的新测序技术的附加类将包装在单独的包中。这将允许一般性改进向上游传播，但不会阻碍每个检测策略的个别需求。

举个例子，看看RNAmodR。RiboMethSeq而且RNAmodR。AlkAnilineSeq包。

功能，未来可能会添加:

• 与我们的包交互以进行数据聚合(例如元基因聚合)
• 与我们的包交互，以进行可视化的下游分析

我们欢迎任何形式的捐款。

6Sessioninfo

sessionInfo ()

## R版本4.2.0 RC (22-04-19 r82224) ##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 20.04.4 LTS ## ##矩阵产品:default ## BLAS: /home/biocbuild/bbs-3.15-bio /R/lib/libRblas. ##因此## LAPACK: /home/biocbuild/bbs-3.15-bio /R/lib/libRlapack。因此## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# [3] LC_TIME=en_GB 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 ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ## [8] base ## ##其他附加的包:# # # # [1] RNAmodR_1.10.0 Modstrings_1.12.0 RNAmodR.Data_1.9.0 [4] ExperimentHubData_1.22.0 AnnotationHubData_1.26.0 futile.logger_1.4.3 # # [7] ExperimentHub_2.4.0 AnnotationHub_3.4.0 BiocFileCache_2.4.0 # # [10] dbplyr_2.1.1 GenomicFeatures_1.48.0 AnnotationDbi_1.58.0 # # [13] Biobase_2.56.0 Rsamtools_2.12.0 Biostrings_2.64.0 # # [16] XVector_0.36.0 rtracklayer_1.56.0 GenomicRanges_1.48.0 # # [19] GenomeInfoDb_1.32.0 IRanges_2.30.0 S4Vectors_0.34.0 # # [22] BiocGenerics_0.42.0 BiocStyle_2.24.0 # # # #通过命名空间加载(不附加):# [1] backports_1.4.1 Hmisc_4.7-0 ## [3] splines_4.2.0 BiocParallel_1.30.0 ## [5] ggplot_2 .2.0 biocparallel_1 .6.29 ## [9] integrationbldb_2 .20.0 htmltools_0.5.2 ## [11] magick_2.7.3 fansi_1.0.3 ## [13] checkmate_2.1.0 magrittr_2.0.3 ## [15] memoise_2.0.1 BSgenome_1.64.0 ## [15] cluster_2.1.3 rocr_1 . 11 ## [21] jpeg_0.1-9 colorspace_2.0-3 ## [23] blob_1.2.3 rappdirs_0.3.3 ## [25] xfun_0.30 dplyr_1.0.8 ## [27] crayon_1.5.1rcurl_1.98 - 1.6 # # [29] jsonlite_1.8.0 graph_1.74.0 # # [31] VariantAnnotation_1.42.0 survival_3.3-1 # # [33] glue_1.6.2 gtable_0.3.0 # # [35] zlibbioc_1.42.0 DelayedArray_0.22.0 # # [37] scales_1.2.0 futile.options_1.0.1 # # [39] DBI_1.1.2 Rcpp_1.0.8.3 # # [41] htmlTable_2.4.0 xtable_1.8-4 # # [43] progress_1.2.2 foreign_0.8 - 82 # # [45] bit_4.0.4 OrganismDbi_1.38.0 # # [47] Formula_1.2-4 AnnotationForge_1.38.0 # # [49] htmlwidgets_1.5.4 httr_1.4.2 # # [51] RColorBrewer_1.1-3 ellipsis_0.3.2 # # [53]## [55] xml_3 .1.0 pkgconfig_2.0.3 ## [55] XML_3.99-0.9 nnet_7 -17 ## [57] Gviz_1.40.0 sass_0.4.1 ## [59] utf8_1.2.2 labeling_0.4.2 ## [61] reshape2_1.4.4 tidyselect_1.1.2 ## [63] rlang_1.0.2 later_1.3.0 ## [65] munsell_1 .5.0 biocViews_1.64.0 ## [67] BiocVersion_3.15.2 tools_4.2.0 ## [67] biocversion_1 .0.6 cli_3.3.0 ## [71] generics_0.1.2 RSQLite_2.2.12 ## [73] evaluate_0.15 string_1 .4.0 ## [75] fastmap_1.1.0 yaml_2.3.5 ## [79] purrr_0.3.4 AnnotationFilter_1.20.0 ## [81] KEGGREST_1.36.0 RBGL_1.72.0 ## [83] mime_0.12 formatR_1.12 ## [85] xml2_1.3.3 biomaRt_2.52.0 ## [87] rstudioapi_0.13 compiler_4.2.0 ## [89] filelock_1.0.2 curl_4.3.2 ## [91] png_0.1-7 interactiveDisplayBase_1.34.0 ## [93] tibble_3.1.6 bslib_0.3.1 ## [95] stringi_1.7.6 highr_0.9 ## [97] lattice_0.20-45 ProtGenerics_1.28.0 ## [99] Matrix_1.4-1 vctrs_0.4.1 ## [101] stringdist_0.9.8 pillar_1.7.0 ## [103] BiocCheck_1.32.0 lifecycle_1.0.1 ## [105] RUnit_0.4.32 BiocManager_1.30.17 ## [107] jquerylib_0.1.4 data.table_1.14.2 ## [109] bitops_1.0-7 colorRamps_2.3 ## [111] httpuv_1.6.5 R6_2.5.1 ## [113] BiocIO_1.6.0 latticeExtra_0.6-29 ## [115] bookdown_0.26 promises_1.2.0.1 ## [117] gridExtra_2.3 codetools_0.2-18 ## [119] dichromat_2.0-0 lambda.r_1.2.4 ## [121] assertthat_0.2.1 SummarizedExperiment_1.26.0 ## [123] rjson_0.2.21 withr_2.5.0 ## [125] GenomicAlignments_1.32.0 GenomeInfoDbData_1.2.8 ## [127] parallel_4.2.0 hms_1.1.1 ## [129] rpart_4.1.16 grid_4.2.0 ## [131] rmarkdown_2.14 MatrixGenerics_1.8.0 ## [133] biovizBase_1.44.0 shiny_1.7.1 ## [135] base64enc_0.1-3 restfulr_0.0.13