7-methyl guanosine (m7G), 3-methyl cytidine (m3C) and Dihydrouridine (D) are commonly found in rRNA and tRNA and can be detected classically by primer extension analysis. However, since the modifications do not interfere with Watson-Crick base pairing, a specific chemical treatment needs to be employed to cause strand breaks specifically at the modified positions. Initially, this involved a sodium borhydride treatment to create abasic sites and cleaving the RNA at abasic sites with aniline.
This classical protocol was converted to a high throughput sequencing method call AlkAnilineSeq and allows modified position be detected by an accumulation of 5’-ends at the N+1 position [Marchand et al. (2018)]. It was found, that m3C is susceptible to this treatment, which allows m7G, m3C and D to be detected by the same method from the same data sets, since the identify of the unmodified nucleotide informs about the three modified nucleotides.
The ModAlkAnilineSeq
class uses the the NormEnd5SequenceData
class to store
and aggregate data along the transcripts. The calculated scores follow the
nomenclature of [Marchand et al. (2018)] with the names scoreNC
(default) and scoreSR
.
## snapshotDate(): 2020-10-02
library(rtracklayer)
library(GenomicRanges)
library(RNAmodR.AlkAnilineSeq)
library(RNAmodR.Data)
The example workflow is limited to 18S rRNA and some tRNA from S.cerevisiae.
As 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.AAS.gff3())
sequences <- RNAmodR.Data.example.AAS.fasta()
files <- list("wt" = c(treated = RNAmodR.Data.example.wt.1(),
treated = RNAmodR.Data.example.wt.2(),
treated = RNAmodR.Data.example.wt.3()),
"Bud23del" = c(treated = RNAmodR.Data.example.bud23.1(),
treated = RNAmodR.Data.example.bud23.2()),
"Trm8del" = c(treated = RNAmodR.Data.example.trm8.1(),
treated = RNAmodR.Data.example.trm8.2()))
The analysis is triggered by the construction of a ModSetAlkAnilineSeq
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).
msaas <- ModSetAlkAnilineSeq(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
msaas
## ModSetAlkAnilineSeq of length 3
## names(3): wt Bud23del Trm8del
## | Modification type(s): m7G / m3C / D
## wt Bud23del Trm8del
## | Modifications found: yes (9) yes (8) yes (7)
## | Settings:
## minCoverage minReplicate find.mod minLength minSignal minScoreNC
## <integer> <integer> <logical> <integer> <integer> <integer>
## wt 10 1 TRUE 9 10 50
## Bud23del 10 1 TRUE 9 10 50
## Trm8del 10 1 TRUE 9 10 50
## minScoreSR minScoreBaseScore scoreOperator
## <numeric> <numeric> <character>
## wt 0.5 0.9 &
## Bud23del 0.5 0.9 &
## Trm8del 0.5 0.9 &
As expected the m7G1575 is missing from the Bud23del samples.
mod <- modifications(msaas)
lapply(mod,head, n = 2L)
## $wt
## GRanges object with 2 ranges and 6 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr1 1575 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## [2] chr3 46 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## score scoreSR Parent
## <numeric> <numeric> <character>
## [1] 162.228 0.984209 1
## [2] 373.773 0.841166 3
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengths
##
## $Bud23del
## GRanges object with 2 ranges and 6 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr3 46 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## [2] chr5 50 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## score scoreSR Parent
## <numeric> <numeric> <character>
## [1] 254.6403 0.858101 3
## [2] 86.3556 0.605249 5
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengths
##
## $Trm8del
## GRanges object with 2 ranges and 6 metadata columns:
## seqnames ranges strand | mod source type
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr1 1575 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## [2] chr3 37 + | m7G RNAmodR.AlkAnilineSeq RNAMOD
## score scoreSR Parent
## <numeric> <numeric> <character>
## [1] 117.2479 0.98729 1
## [2] 69.9604 0.97953 3
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengths
As outlined in the RNAmodR
package we can compare the samples using the
plotCompareByCoord
to prepare a heatmap. For this we select some position
from the found modifications. In addition we prepare an alias table.
coord <- mod[[1L]]
alias <- data.frame(tx_id = c(1L,3L,5L,6L,7L,8L,10L,11L),
name = c("18S rRNA","tF(GAA)B","tG(GCC)B","tT(AGT)B",
"tQ(TTG)B","tC(GCA)B","tS(CGA)C","tV(AAC)E1"),
stringsAsFactors = FALSE)
plotCompareByCoord(msaas, coord, score = "scoreSR", alias = alias,
normalize = TRUE)
plotCompareByCoord(msaas, coord[1L], score = "scoreSR", alias = alias)