RNAmodR: analyzing high throughput sequencing data for post-transcriptional RNA modification footprints
2023-04-25
Package
RNAmodR 1.14.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 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. With these advances post-transcriptional modifications were also detected on mRNA. Among these high throughput techniques are for example Pseudo-Seq (Carlile et al. 2014), RiboMethSeq (Birkedal et al. 2015) and AlkAnilineSeq (Marchand et al. 2018) each able to detect a specific type of modification from footprints in RNA-Seq data prepared with the selected methods.
Since similar pattern can be observed from some of these techniques, overlaps of the bioinformatical pipeline already are and will become more frequent with new emerging sequencing techniques.
RNAmodR implements classes and a workflow to detect post-transcriptional RNA
modifications in high throughput sequencing data. It is easily adaptable to new
methods and can help during the phase of initial method development as well as
more complex screenings.
Briefly, from the SequenceData, specific subclasses are derived for accessing
specific aspects of aligned reads, e.g. 5’-end positions or pileup data. With
this a Modifier class can be used to detect specific patterns for individual
types of modifications. The SequenceData classes can be shared by different
Modifier classes allowing easy adaptation to new methods.
## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationDbi'## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationForge'## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'ExperimentHubData'library(rtracklayer)
library(Rsamtools)
library(GenomicFeatures)
library(RNAmodR.Data)
library(RNAmodR)1.1 SequenceData
Each SequenceData object is created with a named character vector, which can
be coerced to a BamFileList, or named BamFileList. The names must be either
“treated” or “control” describing the condition the data file belongs to.
Multiple files can be given per condition and are used as replicates.
annotation <- GFF3File(RNAmodR.Data.example.gff3())
sequences <- RNAmodR.Data.example.fasta()
files <- c(Treated = RNAmodR.Data.example.bam.1(),
Treated = RNAmodR.Data.example.bam.2(),
Treated = RNAmodR.Data.example.bam.3())For annotation and sequences several input are accepted. annotation can
be a GRangesList, a GFF3File or a TxDb object. Internally, a GFF3File
is converted to a TxDb object and a GRangesList is retrieved using the
exonsBy function.
seqdata <- End5SequenceData(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
## Loading 5'-end position data from BAM files ... OKseqdata## End5SequenceData with 60 elements containing 3 data columns and 3 metadata columns
## - Data columns:
## end5.treated.1 end5.treated.2 end5.treated.3
## <integer> <integer> <integer>
## - Seqinfo object with 84 sequences from an unspecified genome; no seqlengths:SequenceData extends from a CompressedSplitDataFrameList and contains the
data per transcript alongside the annotation information and the sequence. The
additional data stored within the SequenceData can be accessed by several
functions.
names(seqdata) # matches the transcript names as returned by a TxDb object
colnames(seqdata) # returns a CharacterList of all column names
bamfiles(seqdata)
ranges(seqdata) # generate from a TxDb object
sequences(seqdata)
seqinfo(seqdata)Currently the following SequenceData classes are implemented:
End5SequenceDataEnd3SequenceDataEndSequenceDataProtectedEndSequenceDataCoverageSequenceDataPileupSequenceDataNormEnd5SequenceDataNormEnd3SequenceData
The data types and names of the columns are different for most of the
SequenceData classes. As a naming convenction a descriptor is combined with
the condition as defined in the files input and the replicate number. For more
details please have a look at the man pages, e.g. ?End5SequenceData.
SequenceData objects can be subset like a CompressedSplitDataFrameList.
Elements are returned as a SequenceDataFrame dependent of the type of
SequenceData used. For each SequenceData class a matching
SequenceDataFrame is implemented.
seqdata[1]## End5SequenceData with 1 elements containing 3 data columns and 3 metadata columns
## - Data columns:
## end5.treated.1 end5.treated.2 end5.treated.3
## <integer> <integer> <integer>
## - Seqinfo object with 84 sequences from an unspecified genome; no seqlengths:sdf <- seqdata[[1]]
sdf## End5SequenceDataFrame with 1649 rows and 3 columns
## end5.treated.1 end5.treated.2 end5.treated.3
## <integer> <integer> <integer>
## 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
## <Rle> <IRanges> <Rle> | <integer> <character> <integer>
## [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...UGCGGUGGGCUUAUAAAUAUCUUAAAUAUUCUUACAThe SequenceDataFrame objects retains some accessor functions from the
SequenceData class.
names(sdf) # this returns the columns names of the data
ranges(sdf)
sequences(sdf)Subsetting of a SequenceDataFrame returns a SequenceDataFrame or
DataFrame, depending on whether it is subset by a column or row, respectively.
The drop argument is ignored for column subsetting.
sdf[,1:2]## End5SequenceDataFrame with 1649 rows and 2 columns
## end5.treated.1 end5.treated.2
## <integer> <integer>
## 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
## <Rle> <IRanges> <Rle> | <integer> <character> <integer>
## [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...UGCGGUGGGCUUAUAAAUAUCUUAAAUAUUCUUACAsdf[1:3,]## DataFrame with 3 rows and 3 columns
## end5.treated.1 end5.treated.2 end5.treated.3
## <integer> <integer> <integer>
## 1 1 4 0
## 2 0 2 0
## 3 0 0 01.2 Modifier
Whereas, the SequenceData classes are used to hold the data, Modifier
classes are used to detect certain features within high throughput sequencing
data to assign the presence of specific modifications for an established
pattern. The Modifier class (and its nucleotide specific subclasses
RNAModifier and DNAModifier) is virtual and can be addapted to individual
methods. For example mapped reads can be analyzed using the ModInosine
class to reveal the presence of I by detecting a A to G conversion in normal
RNA-Seq data. Therefore, ModInosine inherits from RNAModifier.
To fix the data processing and detection strategy, for each type of sequencing
method a Modifier class can be developed alongside to detect modifications.
For more information on how to develop such a class and potentially a new
corresponding SequenceData class, please have a look at the vignette for
creating a new analysis.
For now three Modifier classes are available:
ModInosineModRiboMethSeqfrom theRNAmodR.RiboMethSeqpackageModAlkAnilineSeqfrom theRNAmodR.AlkAnilineSeqpackage
Modifier objects can use and wrap multiple SequenceData objects as elements
of a SequenceDataSet class. The elements of this class are different types of
SequenceData, which are required by the specific Modifier class. However,
they are required to contain data for the same annotation and sequence data.
Modifier objects are created with the same arguments as SequenceData objects
and will start loading the necessary SequenceData objects from these. In
addition they will automatically start to calculate any additional scores
(aggregation) and then start to search for modifications, if the optional
argument find.mod is not set to FALSE.
mi <- ModInosine(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
## Loading Pileup data from BAM files ... OK
## Aggregating data and calculating scores ... Starting to search for 'Inosine' ... done.(Hint: If you use an artificial genome, name the chromosomes chr1-chrN. It
will make some things easier for subsequent visualization, which relies on the
Gviz package)
Since the Modifier class wraps a SequenceData object the accessors to data
contained within work similarly to the SequenceData accessors described above.
What type of conditions the Modifier class expects/supports is usually
described in the man pages of the Modifier class.
names(mi) # matches the transcript names as returned by a TxDb object
bamfiles(mi)
ranges(mi) # generated from a TxDb object
sequences(mi)
seqinfo(mi)
sequenceData(mi) # returns the SequenceData 1.2.1 Settings
The behavior of a Modifier class can be fine tuned using settings. The
settings() function is a getter/setter for arguments used in the analysis and
my differ between different Modifier classes depending on the particular
strategy and whether they are implemented as flexible settings.
settings(mi)## $minCoverage
## [1] 10
##
## $minReplicate
## [1] 1
##
## $find.mod
## [1] TRUE
##
## $minScore
## [1] 0.4settings(mi,"minScore")## [1] 0.4settings(mi) <- list(minScore = 0.5)
settings(mi,"minScore")## [1] 0.51.3 ModifierSet
Each Modifier object is able to represent one sample set with multiple
replicates of data. To easily compare multiple sample sets the ModifierSet
class is implemented.
The ModifierSet object is created from a named list of named character vectors
or BamFileList objects. Each element in the list is a sample type with a
corresponding name. Each entry in the character vector/BamFileList is a
replicate (Alternatively a ModifierSet can also be created from a list of
Modifier objects, if they are of the same type).
sequences <- RNAmodR.Data.example.AAS.fasta()
annotation <- GFF3File(RNAmodR.Data.example.AAS.gff3())
files <- list("SampleSet1" = c(treated = RNAmodR.Data.example.wt.1(),
treated = RNAmodR.Data.example.wt.2(),
treated = RNAmodR.Data.example.wt.3()),
"SampleSet2" = c(treated = RNAmodR.Data.example.bud23.1(),
treated = RNAmodR.Data.example.bud23.2()),
"SampleSet3" = c(treated = RNAmodR.Data.example.trm8.1(),
treated = RNAmodR.Data.example.trm8.2()))msi <- ModSetInosine(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 ... OKThe creation of the ModifierSet will itself trigger the creation of a
Modifier object each containing data from one sample set. This step is
parallelized using the BiocParallel package. If a
Modifier class itself uses parallel computing for its analysis, it is switched
off unless internalBP = TRUE is set. In this case each Modifier object is
created in sequence, allowing parallel computing during the creation of each
object.
names(msi)## [1] "SampleSet1" "SampleSet2" "SampleSet3"msi[[1]]## A ModInosine object containing PileupSequenceData with 11 elements.
## | Input files:
## - treated: /home/biocbuild/.cache/R/ExperimentHub/979d453092f9b_2544
## - treated: /home/biocbuild/.cache/R/ExperimentHub/979d410801920_2546
## - treated: /home/biocbuild/.cache/R/ExperimentHub/979d49ef08d2_2548
## | Nucleotide - Modification type(s): RNA - I
## | Modifications found: yes (6)
## | Settings:
## minCoverage minReplicate find.mod minScore
## <integer> <integer> <logical> <numeric>
## 10 1 TRUE 0.4Again accessors remain mostly the same as described above for the Modifier
class returning a list of results, one element for each Modifier object.
bamfiles(msi)
ranges(msi) # generate from a TxDb object
sequences(msi)
seqinfo(msi)2 Analysis of detected modifications
Found modifications can be retrieved from a Modifier or ModifierSet object
via the modifications() function. The function returns a GRanges or
GRangesList object, respectively, which contains the coordinates of the
modifications with respect to the genome used. For example if a transcript
starts at position 100 and contains a modified nucleotide at position 50 of the transcript, the
returned coordinate will 150.
mod <- modifications(msi)
mod[[1]]## GRanges object with 6 ranges and 5 metadata columns:
## seqnames ranges strand | mod source type score
## <Rle> <IRanges> <Rle> | <character> <character> <character> <numeric>
## [1] chr2 34 + | I RNAmodR RNAMOD 0.900932
## [2] chr4 35 + | I RNAmodR RNAMOD 0.899622
## [3] chr6 34 + | I RNAmodR RNAMOD 0.984035
## [4] chr7 67 + | I RNAmodR RNAMOD 0.934553
## [5] chr9 7 + | I RNAmodR RNAMOD 0.709758
## [6] chr11 35 + | I RNAmodR RNAMOD 0.874027
## Parent
## <character>
## [1] 2
## [2] 4
## [3] 6
## [4] 7
## [5] 9
## [6] 11
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengthsTo retrieve the coordinates with respect to the transcript boundaries, use the
optional argument perTranscript = TRUE. In the example provided here, this
will yield the same coordinates, since a custom genome was used for mapping of
the example, which does not contain transcripts on the negative strand and per
transcript chromosomes.
mod <- modifications(msi, perTranscript = TRUE)
mod[[1]]## GRanges object with 6 ranges and 5 metadata columns:
## seqnames ranges strand | mod source type score
## <Rle> <IRanges> <Rle> | <character> <character> <character> <numeric>
## [1] chr2 34 * | I RNAmodR RNAMOD 0.900932
## [2] chr4 35 * | I RNAmodR RNAMOD 0.899622
## [3] chr6 34 * | I RNAmodR RNAMOD 0.984035
## [4] chr7 67 * | I RNAmodR RNAMOD 0.934553
## [5] chr9 7 * | I RNAmodR RNAMOD 0.709758
## [6] chr11 35 * | I RNAmodR RNAMOD 0.874027
## Parent
## <character>
## [1] 2
## [2] 4
## [3] 6
## [4] 7
## [5] 9
## [6] 11
## -------
## seqinfo: 11 sequences from an unspecified genome; no seqlengths2.1 Compairing results
To compare results between samples, a ModifierSet as well as a definition of
positions to compare are required. To construct a set of positions, we will use
the intersection of all modifications found as an example.
mod <- modifications(msi)
coord <- unique(unlist(mod))
coord$score <- NULL
coord$sd <- NULL
compareByCoord(msi,coord)## DataFrame with 6 rows and 6 columns
## SampleSet1 SampleSet2 SampleSet3 names positions mod
## <numeric> <numeric> <numeric> <factor> <factor> <character>
## 1 0.900932 0.998134 0.953651 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.943773 7 67 I
## 5 0.709758 0.766484 0.681451 9 7 I
## 6 0.874027 0.971474 0.954782 11 35 IThe result can also be plotted using plotCompareByCoord, which accepts an
optional argument alias to allow transcript ids to be converted to other
identifiers. For this step it is probably helpful to construct a TxDb object
right at the beginning and use it for constructing the Modifier/ModifierSet
object as the annotation argument.
txdb <- makeTxDbFromGFF(annotation)## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ... OKalias <- data.frame(tx_id = names(id2name(txdb)),
name = id2name(txdb))plotCompareByCoord(msi, coord, alias = alias)
Figure 1: Heatmap for identified Inosine positions
Additionally, the order of sample sets can be adjusted, normalized to any of the sample sets and the numbering of positions shown per transcript.
plotCompareByCoord(msi[c(3,1,2)], coord, alias = alias, normalize = "SampleSet3",
perTranscript = TRUE)