Method for detection and identification of RNA modifications

Abstract

Disclosed is a method of pre-processing and inferencing data for RNA chemical modification detection. The method employs direct RNA sequencing on synthetic RNA sequences, various outputs. Basecalled RNA reads are aligned to reference synthetic RNA sequences to obtain sequence alignment reads. Nucleotide assignments of the sequence alignment reads are then refined, data segmented and the segmented data aligned to the individual k-mers. Raw signal features corresponding to specific nucleotide positions in the reference synthetic RNA sequence are then identifiable. The features of the raw signals are subjected to a multiple-instance supervised learning model training framework for the detection of one RNA chemical modification or multiple RNA chemical modifications. The trained model infers from an unseen RNA sequence data set the likelihood of one RNA chemical modification or multiple RNA chemical modifications.

Description

[0001] Method for Detection and Identification of RNA Modifications

[0002] Technical Field

[0003] The present invention relates, in general terms, to the field of molecular biology and nucleic acid analysis. In particular, this disclosure provides methods for detection and identification of modifications on ribonucleic acid (RNA) molecules.

[0004] Background

[0005] This background is provided for generally presenting the context of the disclosure. Contents of this background section are neither expressly nor implied admitted as prior art against the present disclosure.

[0006] RNA modifications play crucial roles in regulating various aspects of RNA function, which ultimately impact cellular processes and organismal biology. RNA modifications can affect gene expression by regulating mRNA stability, splicing, localization, and translation efficiency. These modifications can finetune gene expression in response to developmental cues, environmental signals and cellular stress. Dysregulation of RNA modifications has been linked to various diseases including cancer, neurological disorders, and metabolic syndromes. For example, aberrant RNA modifications can disrupt normal gene expression patterns, leading to the development and progression of diseases.

[0007] RNA modifications play a critical role in the development and efficiency of RNA therapeutics, particularly mRNA vaccines. RNA modifications such as Nl- methyl-pseudouridine (mlMJ) are a vital ingredient in mRNA vaccines to reduce innate immune response to ensure vaccine efficacy and safety. Hence, efficient RNA modification detection is necessary for routine mRNA vaccine quality checks. Nanopore direct RNA sequencing, coupled with machine-learning approaches, has been shown to capture the presence of RNA modifications. Several computational methods have been developed that can specifically identify the type of modification from direct RNA-Seq data using supervised machine learning. These approaches require direct RNA-Seq data of a sample with known modifications that can be used to train the underlying machine learning models. These models incorporate features related to the sequence context of modified position and features from the nanopore signal data, which are then used to discriminate between modified and unmodified bases. To obtain data with known modified positions, experimental approaches such as m6ACE- Seq or GLORI-Seq can be used, or comparative approaches can identify modified positions from direct RNA-Seq data after knockout or inhibition of writer enzymes. As an alternative strategy, synthetic RNAs with specific modifications have been generated that enable the model training potentially even for modifications where otherwise transcriptome-wide maps of modifications are not available.

[0008] Using this approach, it has been demonstrated that m6A, m5C, mlA, MJ and other modifications can be detected. Furthermore, methods have been developed to jointly identify modifications such as m6A and J, enabling the profiling of different modifications within the same RNA molecule. While the identification of modified bases such as m6A has been reported with an accuracy that is comparable to alternative technologies, the discrimination of modifications such as MJ and mlMJ that can occur at the same base has been challenging for methods which are largely trained to detect single RNA modifications.

[0009] However, synthetic RNAs can be substantially different from cellular RNA due to a higher density of modified bases, the use of artificial base modifications, or an unusual sequence composition, which might not be reflected in such training data. Therefore, it would be desirable especially for the research and development community and pharmaceutical companies to have highly confident identification of multiple RIMA modifications to improve therapeutics.

[0010] Summary

[0011] Disclosed is a method for ribonucleic acid (RNA) modification detection, the method comprising: subjecting modified and unmodified synthetic RNA samples to direct RNA sequencing, thereby producing an RNA sequence data set comprising basecalled reads and raw electrical signal-level data; aligning the RNA sequence data set to reference synthetic RNA sequences to obtain sequence alignment reads indexed to specific sites in the reference synthetic RNA sequences; segmenting the raw electrical signal-level data into discrete events corresponding to nucleotide bases; creating an index that links each base in the sequence alignment reads with the segmented raw electrical signal-level data; aligning the segmented raw electrical signal-level data to the sequence alignment reads to produce features of the raw electrical signal-level data corresponding to specific nucleotide positions in the reference synthetic RNA sequence; labelling the modified and unmodified synthetic RNA samples using labels, based on the features at the specific nucleotide positions; and outputting labelled datasets, each comprising all RNA sequences having a particular said modification determined by the labels.

[0012] Further disclosed is a method for training a detection model for detecting RNA modifications, the method comprising: producing labelled datasets according to the method described above; for each labelled dataset, generating first summarised data for all sites of the respective base modification; for each unmodified RIMA sequence, generating second summarised data for all sites of the respective unmodified RNA sequence corresponding to the sites of the respective base modification, wherein the first summarised data and second summarised data comprise the features; and separately inputting the first summarised data and the second summarised data into a supervised machine learning based RNA modification detection model to train the detection model to predict presence of a modification at each site corresponding to the base.

[0013] Additionally, the detection model may include one or more single modification detection models, each trained using ones of said datasets corresponding to a respectively different single said base modification, and the unmodified RNA sequences.

[0014] More preferably, the detection model may include one or more multiple modification detection models, each trained using ones of said datasets corresponding to said modified sequences each having more than one said base modification, and the unmodified RNA sequences.

[0015] Further, the RNA modifications type may comprise one or more of Pseudouridine (pU or i ), 1-methylpseudouridine (mlpU or mlUJ), 5-methoxyuridine (mo5U), 2-thiouridine (s2U), 5-methylcytidine (m5C), 5-formylcytidine (f5C), N6- methyladenosine (m6A), Nl-methyladenosine (mlA), and thienoguanosine (thG) and unmodified nucleotides.

[0016] Furthermore, the features of the raw electrical signal-level data may include one or more normalized mean ( ), standard deviation (o), and dwell time (I) for the current and adjacent k-mers. Preferably, the direct RNA sequencing method is nanopore sequencing.

[0017] Further described is a computer-implemented method for RNA chemical modification detection, comprising: training one or more first detection model using unmodified and modified synthetic sequences to detect a single modification, using labelled datasets outputted from the method described above; training one or more second detection model using unmodified and modified synthetic sequences to detect multiple RNA modifications, using labelled datasets outputted from the method described above; subjecting an unseen RNA sequence data set to the first and second detection models, to detect one or more single RNA modifications or multiple RNA modifications in the unseen RNA sequence data set (detected RNA modifications); categorizing the detected RNA modifications into at least one specific RNA modification type; and generating, for each detected RNA modification, confidence probability that a base at a specific nucleotide position corresponding to the respective detected RNA modification in each read in the unseen RNA sequence set has been modified.

[0018] Additionally disclosed is a system for RNA modification detection, comprising: memory; at least one processor(s); and a pre-processing module, wherein the memory stores instructions that, when executed by the processor(s), cause the pre-processing module to perform any one of the method described above.

[0019] Also disclosed is system for RNA chemical modification detection, comprising: memory; at least one processor(s); and a detection module comprising a plurality of RNA modification detection models, at least one said model being trained to detect a single RNA modification and at least one said model being trained to detect multiple RNA modifications, wherein the memory stores instructions that, when executed by the processor(s), cause the processor(s) to: subjecting an unseen RNA sequence data set to the plurality of RNA modification detection models, to predict presence of a single modification or multiple modifications in the unseen RNA sequence data set (detected RNA modifications); categorizing the detected RNA modifications into at least one specific RNA modification type; and generating, for each detected RNA modification, a confidence probability that a base at a specific nucleotide position corresponding to the respective detected RNA modification in each read in the unseen RNA sequence set has been modified.

[0020] Brief description of the drawings

[0021] Some embodiments of systems and methods for detection and identification of RNA modifications, in accordance with the present disclosure, will now be described, by way of non-limiting example, with reference to the drawings in which:

[0022] Figure 1 illustrates an overview workflow diagram, according to various embodiments of the invention;

[0023] Figure 2 is a flow diagram of a method for ribonucleic acid (RNA) modification detection, according to various embodiments of the invention;

[0024] Figure 3 depicts a diagram representing the architecture of the invention for the prediction of RNA modifications, according to various embodiments of the invention; Figure 4a depicts components of synLibvl library containing 150- nucleotide synthetic sequences representing fragments of RIMA oligos with all k- mers (left) and fragments of the human transcriptome, including 5'UTR, CDS, and 3'UTR, according to various embodiments of the invention;

[0025] Figure 4b depicts nine RNA modifications incorporated into the synLibvl synthetic sequences through in vitro transcription (IVT), which then subjected to direct RNA-sequencing, according to various embodiments of the invention;

[0026] Figure 4c depicts the total number of sequenced reads and aligned reads across the IVT samples (nine modifications and one unmodified) from the RNA004 chemistry, according to various embodiments of the invention;

[0027] Figure 4d illustrates the percentage of synLibvl sequences that are covered by a minimum of twenty (20) sequencing reads across various in IVT samples generated using RNA004, according to various embodiments of the invention;

[0028] Figure 4e illustrates the number of nucleotide sites derived from k-mers containing a single modified base, wherein each site is supported by more than twenty (20) sequencing reads in RNA004, in accordance with various embodiments of the present invention.

[0029] Figure 4f illustrates the percentage of sequence identity between the reference nucleotides (adenine [A], cytosine [C], guanine [G], and uracil [U]) and the corresponding basecalled nucleotides that are mapped to the reference nucleotides in each sample, in accordance with various embodiments of the present invention.

[0030] Figure 5a illustrates the receiver operating characteristic (ROC) curves and precision-recall curves for the modNet single-modification models corresponding to nine RNA modifications: MJ, mlUJ, mo5U, s2U, m5C, f5C, m6A, mlA, and thG, in accordance with various embodiments of the present invention;

[0031] Figure 5b illustrates the ROC and area under precision-recall (PR AUC) curves for the modNet single-modification models, evaluated on held-out test data obtained from a predefined train-test split, as well as on an independent validation batch comprising newly sequenced samples of i , ml|P, mo5U, and s2U, in accordance with various embodiments of the present invention;

[0032] Figure 5c illustrates the receiver operating characteristic area under the curve (ROC AUC) values for the modNet multi-modification models, showing performance in distinguishing between modified and unmodified bases, as well as in differentiating among multiple distinct modification types, in accordance with various embodiments of the present invention;

[0033] Figure 5d illustrates ROC AUCs of the single-modification modNet models for discriminating unmodified and modified bases, in accordance with various embodiments of the present invention;

[0034] Figure 5e and 5f illustrate ROC precision-recall curves of the multimodification A-modification, C-modification, and U-modification modNet models, dorado's m6A, m5C, and i models, and the single-modification modNet-m6A, modNet-m5C, and modNet-’P models when multiple base modifications are present in the data, in accordance with various embodiments of the present invention;

[0035] Figure 5g illustrates ROC curves of m6anet, the single-modification m6A modNet model, and the multi-modification A modification modNet model applied to direct RNA-Seq data from HEK293T cells (DRACH sites), using m6ACE-Seq to define m6A sites. The modNet models were trained on the full IVT data set as they were applied to an unseen sample, in accordance with various embodiments of the present invention;

[0036] Figure 5h illustrates illustrates ROC AUCs of each DRACH kmer and all DRACH k-mers combined, in accordance with various embodiments of the present invention;

[0037] Figure 6a illustrates comparison of the number of all sequenced reads and the number of aligned reads in RNA002, in accordance with various embodiments of the present invention;

[0038] Figure 6b illustrates comparison of the number of alignment reads in RNA004 using the default alignment parameters and the sensitive alignment parameters, in accordance with various embodiments of the present invention;

[0039] Figure 6c illustrates percentage of synLibvl sequences in each IVT sample in RNA002, in accordance with various embodiments of the present invention; Figure 6d illustrates number of single-modified-base-containing sites supported by more than 20 reads in RNA002, in accordance with various embodiments of the present invention;

[0040] Figure 6e discloses percentage of single-modified-base kmer represented in each IVT R.NA004 sample, in accordance with various embodiments of the present invention;

[0041] Figure 6f discloses percentage of multi-modified-base kmer represented in each IVT RNA004 sample, in accordance with various embodiments of the present invention;

[0042] Figure 6g illustrates sequence identity of the basecalled mlUJ sites from Gunter et al. to the reference sites, in accordance with various embodiments of the present invention;

[0043] Figure 7a illustrates ROC and precision-recall curves of the RNA002-based single-modification modNet models of nine modifications (>P, mUP, mo5U, and s2U, m5C, f5C, m6A, mlA, and thG), in accordance with various embodiments of the present invention;

[0044] Figure 7b illustrates ROC and precision-recall curves of the singlemodification IP, mllJ, mo5U, and s2U modNet models in the new batch of UJ, mllJ, mo5U, and s2U samples for validation, in accordance with various embodiments of the present invention;

[0045] Figure 7c illustrates ROC AUCs of the RNA002-based multi-modification modNet models across unmodified and different modifications, in accordance with various embodiments of the present invention;

[0046] Figure 7d illustrates ROC AUCs of the RNA002-based single-modification modNet models between unmodified and modified, in accordance with various embodiments of the present invention;

[0047] Figure 7e illustrates ROC curves of the RNA002-based single-modification m6A modNet model and CHEUI m6A model, ROC curves of the RNA002-based single-modification m5C modNet model and CHEUI m5C model, ROC curves of the RNA004-based single-modification m6A modNet model and dorado's m6A model, and ROC curves of the RNA004-based single-modification ip modNet model and dorado's >P model, in accordance with various embodiments of the present invention;

[0048] Figure 7f illustrates ROC curves of dorado's m6A model applied on m6A and mlA samples. ROC curves of dorado's J model applied on MJ, mlUJ, mo5U, and s2U samples, in accordance with various embodiments of the present invention;

[0049] Figure 7g illustrates ROC curves of the single-modification m6A modNet model applied on m6A and mlA samples, and ROC curves of the singlemodification mlA modNet model applied on mlA and m6A samples, in accordance with various embodiments of the present invention;

[0050] Figure 7h illustrates ROC curves of m6anet, the single-modification m6A modNet model, and the multi-modification A-modification modNet model applied to HEK293T, in accordance with various embodiments of the present invention;

[0051] Figure 8a illustrates illustrates ROC AUCs of each DRACH kmer and all DRACH k-mers combined for the RNA002 chemistry, in accordance with various embodiments of the present invention;

[0052] Figure 8b illustrates illustrates ROC curves of each DRACH kmer and all DRACH k-mers combined for the RNA002 chemistry, in accordance with various embodiments of the present invention.

[0053] Detailed description

[0054] The present invention relates to systems and methods for the computational detection and classification of post-transcriptional ribonucleic acid (RNA) modifications from sequencing data. In particular, the present invention provides a machine learning framework, hereinafter referred to as modNet, that infers multiple RNA modifications from direct RNA sequencing (RNA-Seq), and which can be optimised for a specific sequence of interest. Throughout the specification, the terms "modNet” and "the present invention" may be used interchangeably to refer to the described system, its components, and methods of use, except where the context clearly indicates otherwise. Hereinafter, the system and method according to the present disclosure will be described in detail with reference to Figure 1 to Figure 8b according to the preferred embodiments. It is to be understood that limiting the description to the preferred embodiments of the invention is merely to facilitate discussion of the present invention. Various modifications may be made, to the preferred embodiments, without departing from the scope of the appended claims.

[0055] The modNet model is trained for predictions on the same sequence context while still maintaining the ability to generalise to unseen sequences to mitigate the impact of low coverage region in training data. The architecture of modNet consist of a raw signal deaggregation layer, a k-mer-embedding layer with 592 input channels, a read-level feature encoder layer, a read-level probability encoder layer outputting the read-level probabilities. Each described layer delivers its output to the next described layer. The single-modification modNet models use a noisy-OR pooling layer while the multi-modification modNet models use a softmax layer to summarise the read-probabilities to a site probability.

[0056] As used herein, the term "sequence context" refers to the specific nucleotide (A, U, C, or G) composition and structure surrounding a particular sequence or motif within the RNA molecule. Each model can either be trained to infer modification probabilities for a single modification or for multiple modifications. modNet provides modules for training, inference, and arguments to use pretrained models for 4J, mlUJ, mo5U, 2-thiouridine (s2U), m5C, 5-formylcytidine (f5C), m6A, Nl-methyladenosine (mlA), and thienoguanosine (thG) modification detection. An embodiment of the present invention enables the identification of nine RNA modifications (i.e. N6-methyladenosine (m6A), 1- methyladenosine (mlA), 5-formylcytidine (f5C), 5-methylcytidine (m5C), Pseudouridine (pU or *+>), 1-methylpseudouridine (mlpU or mlUJ), 5- methoxyuridine (mo5U), 2-thiouridine (s2U) and thienoguanosine (thG)) from direct RNA-sequencing data. RNA modifications are discriminated at the nucleotide resolution for single RNA molecules, enabling new avenues in scientific investigation, manufacturing quality control, and therapeutic development. When applied to the raw electrical signals generated from nanopore sequencing, embodiments of the present invention are able to convert the electrical signals to RNA sequences and base modifications.

[0057] Design and generation of synthetic sequences for machine learning (ML) training

[0058] To facilitate the training and evaluation of modNet, a synthetic library designated as synLibVl was constructed. This library comprises two primary components: (1) 150 base pair (bp) fragments derived from human complementary DNA (cDNA) sequences, and (2) 150 bp synthetic sequences engineered to comprehensively represent all possible 5-mer nucleotide combinations while being explicitly designed to exclude any sequences present within the known human transcriptome. Reference is made to Figure 4a for an illustration of the library design.

[0059] The synLibVl library served as a substrate for the generation of RNA molecules through in vitro transcription (IVT). Two categories of RNA products were synthesized: (i) unmodified RNA transcripts that retained standard nucleotides, and (ii) modified RNA transcripts wherein one or more nucleotide positions, and preferably each nucleotide position, were / was systematically substituted with a chemically modified base during the IVT process. These RNA products were subsequently analyzed via direct RNA sequencing - e.g., utilizing the RNA002 and / or RNA004 chemistry as shown in Figure 4b.

[0060] To further validate the performance of modNet and assess the impact of specific RNA modifications, independent validation datasets were generated. These datasets were prepared for a subset of selected chemical modifications, providing additional control data for comprehensive performance benchmarking. Sub-library 1 : Including all possible k-mers with curlcake

[0061] In one embodiment, the method utilizes publicly available platforms - e.g., curlcake, and RNAshapes version 2.1.6 - to generate synthetic modified and unmodified sequences. Curlcake is a synthetic sequence designing software that includes all possible kmers (k=5) and RNAshapes minimizes secondary RIMA structure. Curlcake was run with the following parameters: k=5, length=50, multiplicity=14 to 20, and attempts=100. A total of seven curlcake constructs were produced, one for each multiplicity parameter specified from 14 to 20. The lengths of the seven constructs are 29,070, 30,630, 33,390, 34,950, 37,230, 37,590, and 43,110 bp, respectively. The curlcake constructs were split into shorter sequences of 150-bp length, each overlapping by 30 nucleotides with the previous 150-bp sequence, resulting in 2,048 individual sequences.

[0062] Sub-library 2: Using human transcriptome fragments

[0063] This sub-library comprises synthetic nucleotide sequences resembling the human transcriptome. The generate these sequences, a selection process was conducted utilizing both cell-type-specific and constitutively expressed transcripts. These transcripts were identified using transcriptomic data - for experimentation purposes, this data was obtained from the Singapore IManopore Expression Project (SG-NEx), which comprises direct RNA-sequencing data from the following eight human cell lines: A549, MCF7, HCT116, HepG2, K562, H9, HEYA8, and HEK293T.

[0064] A transcript was considered to be expressed in a given cell line if it was supported by a minimum of n > 1, presently n = five (5), aligned read counts in each biological replicate of the direct RNA-seq datasets corresponding to that cell line. Read counts were determined based on the number of aligned sequencing reads per transcript, as obtained from transcriptome alignment BAM files downloaded from the SG-NEx project. Transcripts were categorized as "constitutively expressed" if they were present in n or more of the aforementioned eight (8) cell lines. Additionally, transcripts exclusively expressed in either the HEK293T or H9 cell line were designated as "HEK293T- specific" and "H9-specific," respectively. Transcripts possessing annotated features including a 5' untranslated region (5'UTR), 3' untranslated region (3'UTR), translational start codon, and translational stop codon were also filtered for further processing.

[0065] For each gene, the transcript exhibiting the highest read count was selected. The selected transcripts were computationally segmented into m-base pair (bp, where 150-bp is used for present purposes) sequences, with each sequence overlapping the preceding sequence by 30 nucleotides. A final filtering step was applied to eliminate any m-bp sequences containing the same nucleotide in a predetermined number of consecutive positions, such as 15 consecutive positions. Each remaining n-bp sequence was annotated according to its relative position within the original transcript, including but not limited to: transcription start site (TSS), 5'UTR, coding sequence (CDS), exon-exon junction, and 3'UTR.

[0066] Final library design

[0067] To enable in vitro transcription of the synthetic nucleotide sequences, a promoter sequence (e.g., T7 promoter sequence) was added to the 5' end of each n-bp synthetic sequence. Furthermore, to facilitate polymerase chain reaction (PCR) amplification, the complete construct— comprising the T7 promoter and the 150-bp synthetic sequence— was flanked at its 5' and 3' termini with primer binding sites. Specifically, the sequences 5'- TTGGACCCTCGTACAGAAGC-3' and 5'-TGGTCTTTGAATAAAGCCTGAGTAGGAAG- 3' were incorporated at the 5' and 3' ends of the construct, respectively. The resulting designed DNA sequences were synthesized and obtained as two distinct oligonucleotide pools from Twist Bioscience (South San Francisco, CA, USA).

[0068] PCR, in vitro transcription and sequencing of the synthetic sequences Both sub-libraries were subjected to PCR amplification (KAPA Hi Fi HotStart PCR, #07958897001) and DNA purification (Zymo Research DNA Clean &. Concentrator, D4013), before being pooled in equimolar concentrations.

[0069] The total number of classifications in the training dataset includes the total number of RNA modifications classifications and one unmodified classification. The training dataset may be a single dataset comprising all sequences with relevant labels (e.g., for training single modification detection models (i.e., for detecting single modifications of a particular base in an RNA sequence), the labels will identify the modification for those sequences including a single modification, and for training multiple modification detection models, the labels will identify the multiple (generally two) modifications in each sequence) or may include multiple datasets. Where multiple datasets are provided, each dataset for training a single modification detection model may include all sequences comprising a modification of a particular base (i.e., a modification of one of A, U, G and C) and the unmodified RNA sequences. Similarly, each dataset for training a multiple modification detection model may include all sequences comprising multiple modifications either of a particular base (i.e., two modifications of one of A, U, G and C), or modifications of two, or at least two, bases, and the unmodified RNA sequences. The pooled synLibVl DNA templates were transcribed in vitro (Thermo Fisher Scientific MEGAscript™ T7 Transcription Kit, AM1334) with unmodified and modified nucleotides. For the unmodified RNA, only unmodified nucleotides were provided during in vitro transcription. For the m6A and mlA RNA samples, the ATP was substituted with either m6A or mlA (TriLink Biotechnologies) while the other 3 nucleotides remained unmodified. The transcribed RNAs were then purified (New England Biolabs Monarch RNA Cleanup Kit, T2040), poly(A)-tailed (New England Biolabs E. coli Poly(A) Polymerase, M0276) and purified again, before they were subjected to nanopore direct RNA sequencing (Oxford Nanopore Technologies, SQK-RNA002 or SQK-RNA004). For the 4>, mlUJ, mo5U, s2U, m5C, and f5C modified RNA samples, a 3' primer carrying a 120-nt poly(T) sequence was used for the PCR amplification of the IVT templates, allowing bypassing of subsequent polyA-tailing reaction and second RNA purification before direct RNA sequencing. The modified nucleotides were also obtained from TriLink Biotechnologies.

[0070] On average, each RNA004 sample has a total of 3,258,611 sequenced reads while each RNA002 sample has a total of 1,000,935 sequenced reads, shown in Figure 4c and Figure 6a. Incorporation of modified bases during IVT was shown to introduce higher basecalling error rates and lower number of aligned reads 49,51,52. To address this, a highly sensitive minimap2 alignment option was used for the IVT sequences, leading to an average number of aligned reads of 1,610,594 per RNA004 sample while 486,981 per RNA002 sample when using modified RNA, as shown in Figure 4c and Figure 6b. Among the 6,756 unique 150 nt sequences in synLibVl, on average, 65% are sequenced with at least 20 reads with RNA004, with thG showing the lowest coverage (4.5%) as shown in Figure 4d and Figure 6c.

[0071] Preprocessing of the direct RNA-seq samples

[0072] Figure 2 discloses the method for ribonucleic acid (RNA) modification detection. The method comprises steps of subjecting modified and unmodified synthetic RNA samples to direct RNA sequencing (step 202), thereby generating an RNA sequence dataset comprising basecalled reads and corresponding raw electrical signal-level data. The method further includes aligning the RNA sequence dataset to reference synthetic RNA sequences to produce sequence alignment reads indexed to specific sites within the reference sequences (step 204). The raw electrical signal-level data is then segmented into discrete events corresponding to individual nucleotide bases (step 206). An index is created that links each base in the sequence alignment reads with the segmented raw electrical signal-level data (step 208). The segmented electrical signal-level data is aligned to the sequence alignment reads to extract features corresponding to specific nucleotide positions in the reference synthetic RNA sequences (step 210). Using these features, the synthetic RNA samples are labeled as modified or unmodified at particular nucleotide positions (step 212). Finally, the method includes outputting labeled datasets (step 214), each comprising all RNA sequences associated with a particular type of RNA modification, as determined by the applied labels.

[0073] Embodiments of the method further warrant subjecting the aligned data to data preparation. For modified RNA sequence, this involves generating summarised data for all sites of the base of modification. For unmodified RNA sequence, this involves generating summarised data for all sites of different bases separately.

[0074] Alignment

[0075] Figure 1 illustrates a nextflow workflow illustrating the pipeline of the present invention according to a preferred embodiment. In some embodiments, the basecalled reads from the sequencing run were aligned to the synLibVl reference with minimap2 version 2.22 with the options -ax map-ont — secondary=no -k5 -w5 -g2k -Al -B2 -02,32 -El,0 -z200 — cap-sw-mem -s 20. The alignment files were converted from sam file format to bam file format, which were then indexed, with samtools version 1.18.

[0076] Segmentation

[0077] As shown in Figure 1, kmer models for data segmentation were specifically trained for the modified samples. In some embodiments, this was achieved using nanopolish train (available in the rlO branch on github) with a lowering of the minimum required read length to work with shorter sequences - this lowering was achieved through modification to the code for nanopolish train. Prior to training the modification-specific kmer models, nanopolish index was used for indexing the fastq reads to the fasts raw signals. Then, nanopolish train was run with the following options: — train-kmer=all —rounds 5 --input-file- name r9.4_70bps.u_to_t_rna.5mer.template. model (the RNA template kmer model available in the nanopolish github repository's etc / r9-model directory). In a preferred embodiment, the fast5 files were converted to blows files with slow5tools version vl.1.0 with the f2s command, and the pod5 files were converted to blow5 files with blue-crab version 0.3.0 with the p2s command. The blow5 files were merged into a blow5 file with the slowStools merge command. Segmentation was performed with f5c version 1.2. The reads from the fastq file and the raw signals from blow5 file were indexed with the f5c index command. Then, f5c eventalign was run with the options — rna --min-mapq 0 - -min-recalib-events 100 -signal-index -scale-event. The — kmer_model argument was used to specify the modification-specific k-mer model for the RNA004 samples.

[0078] Data preparation

[0079] The output from f5c eventalign f5c eventalign (f5c call-methylation detects the methylated cytosine and f5c eventalign aligns raw nanopore signals (events) to the reference k-mers) was summarized for each read and / or each site. This was performed using modnet dataprep, which outputs a json file containing sites of the base of the modification (e.g. for m5C, only NNCNN sites were outputted). For the unmodified sample, sites of all bases were outputted separately into four json files (one for NNANN, NNCNN, NNUNN, and NNGNN) using modnet dataprep.

[0080] The f5c eventalign tool was utilized to align the current signals obtained from direct RNA sequencing to the corresponding reference sequence for all mapped reads. Following alignment, sequence positions supported by more than 20 reads were identified as suitable for model training. In the case of modified RNA, the number of such sites ranged from 2,345 for thG to 238,251 for m5C. For unmodified RNA, the average number of qualifying sites per base was approximately 281,213, as shown in Figure 4e and Figure 6d. These sites encompass all possible 5-mer sequences for each modification, with the exception of thG, which covered approximately 94% of the full 5-mer set as shown in Figures 6e and 6f. Although the modified RNA samples exhibited a lower number of aligned reads compared to unmodified samples, the datasets provide a sufficient number of informative sites for the effective training and evaluation of models for RIMA modification detection.

[0081] To evaluate whether the identified sites correspond to modified bases or unmodified bases that may have been erroneously incorporated during IVT by T7 RNA polymerase (RNAP), the sequence identity of all aligned reads was examined. The unmodified RNA sample demonstrated an average sequence identity of 98.7% at each base, indicating a low base incorporation error rate by T7 RNAP. In modified RNA samples, a reduction in sequence identity was observed, with unmodified bases showing an average identity of 96.6%, and modified bases exhibiting a lower identity of 72.2% (as shown in Figure 4f). With the exception and ml4J, which have been previously reported to be miscalled as cytidine (C) during basecalling (refer to Figure 6g), all other modified nucleotides were correctly identified according to their original base identities (as shown in Figure 4f). These results indicate that the dataset primarily consists of accurately transcribed and basecalled modified RNAs, which are suitable for use in the training and evaluation of RNA modification detection models.

[0082] RNA modification models

[0083] Figure 3 illustrates a diagram representing the modNet framework, disclosing amethod for ribonucleic acid (RNA) modification detection, and inferencing data for RNA chemical modification detection. The framework comprises a module modNet) designed for the prediction or detection (used synonymously herein, unless context indicates otherwise) of RNA modifications. modNet is optimized for processing new sequencing data derived from either known or unknown sequence contexts. Illustrated in Figure 3 are the A-modification models. Singlemodification and multi-modification modNet models use either a noisy-OR pooling layer or a softmax layer, respectively, to summarise the readprobabilities to a site probability. In another aspect, the present disclosure provides a method for training a detection model for identifying RIMA modifications. The method includes producing labelled datasets in accordance with the method described above for RNA modification detection. For each labelled dataset corresponding to a particular RNA base modification, the method further comprises generating first summarized data for all nucleotide sites associated with the respective base modification. Additionally, for each unmodified RNA sequence, second summarized data is generated for all sites corresponding to the same nucleotide positions as those associated with the respective base modification. Both the first summarized data and the second summarized data include features derived from raw electrical signal-level data aligned to specific nucleotide positions. The method further comprises separately inputting the first and second summarized data into a supervised machine learning-based RNA modification detection model, thereby training the model to predict the presence or absence of a modification at each corresponding nucleotide site.

[0084] The detection model of the present invention may include one or more single modification detection models, each trained using ones of said datasets corresponding to a respectively different single said base modification, and the unmodified RNA sequences.

[0085] The detection model of the present invention may include one or more multiple modification detection models, each trained using ones of said datasets corresponding to said modified sequences each having more than one said base modification, and the unmodified RNA sequences.

[0086] Single-modification detection models

[0087] Single-modification modNet models were trained using direct RNA sequencing data for nine modification types with RNA004 chemistry and eight modification types with RNA002 chemistry, as illustrated in Figure 7a. The sequencing data were partitioned into independent training and test datasets, wherein 50% of the reads from all sequences were allocated for model training and the remaining 50% were reserved for model evaluation. Model performance was assessed using the Area Under the Curve (AUC) metrics derived from both Precision-Recall (PR) curves and Receiver Operating Characteristic (ROC) curves.

[0088] Across all modification types, the modNet models demonstrated an average PR AUC of 0.88 and an average ROC AUC of 0.923 when evaluated on the test dataset (refer to Figure 5a). To assess the generalizability of the models to independent sequencing experiments, the single-modification modNet models corresponding to ip, mlUJ, mo5U, and s2U were applied to an external validation dataset. On this dataset, the models achieved an average PR AUC of 0.924 and an average ROC AUC of 0.969 (shown in Figure 5b and Figure 7b). These results demonstrate that the single-modification models developed using the modNet framework are capable of accurately detecting multiple RNA modifications across all four nucleotide bases and across distinct experimental batches.

[0089] (1 ) Model architecture

[0090] In some embodiments, modNet expands the model architecture in m6Anet, which was optimised for m6A, to enable the detection of any RNA modification, The original mSAnet model architecture contains a raw signal deaggregation layer, a k-mer-embedding layer with 66 input channels, a read-level feature encoder layer, a read-level probability encoder layer outputting the read-level probabilities, and a noisy-OR pooling layer summarizing the read-probabilities to a site probability. To facilitate the identification of any modification, the original k-mer embedding layer was adopted to allow training and inference on any k-mer. Specifically, the k-mer-embedding layer in modNet was specified to have a number of input channels corresponding to the number of possible rimers - e.g., 592 input channels to support all 592 possible 5-mers from the NNNXNNN sequence context k-mers (X can be A, C, U, or G).

[0091] (2) Estimation of model parameters (training) The model parameters are updated subject to minimizing the binary cross entropy loss as described for m6Anet.

[0092] Subsequently, multi-modification modNet models were developed to distinguish among RNA modifications that occur at the same nucleotide base. Specifically, the following models were trained: modNet-A to classify Nl-methyladenosine (mlA) and N6-methyladenosine (m6A); modNet-C to classify 5-methylcytidine (m5C) and 5-formylcytidine (fSC); and modNet-U to classify pseudouridine (ip), 1-methylpseudouridine (mlV), 5-methoxyuridine (mo5U), and 2-thiouridine (s2U). Each multi-modification model produces class probabilities for each site, such that the sum of all class probabilities (including the unmodified class) equals one.

[0093] Consistent with the training procedure used for single-modification models, multi-modification models were trained using 50% of the available reads and evaluated using the remaining 50%. Performance evaluation was conducted using a one-vs-one classification approach, wherein each modification class was individually compared against all other classes. Evaluation of the multimodification models demonstrated that the cytidine- and uridine-targeting models (modNet-C and modNet-U, respectively) achieved clear separation among the modification classes, with average ROC AUC values of 0.99 and 0.91, respectively (see Figure 5c; Figure 7c). By contrast, the adenosine-targeting model (modNet-A) achieved an ROC AUC of 0.63 when distinguishing m6A from mlA. This reduced performance may be attributed to the potential rearrangement of mlA into m6A during RIMA preparation and sequencing (i.e., Dimroth rearrangement), or may reflect inherent challenges in distinguishing specific modifications using direct RIMA sequencing data (see Figure 5c; Figure 7C).

[0094] Further evaluation was conducted to assess the capability of the multimodification models to differentiate between modified and unmodified nucleotide bases. In this context, the models achieved an average ROC AUC of 0.977, surpassing the average ROC AUC of 0.923 observed for the singlemodification modNet models (see Figures 5c, 5d; Figure 7c, 7d). This improvement in performance may be attributed to the increased volume of training data available due to the inclusion of multiple modification types per base. Collectively, these results demonstrate that the multi-modification modNet models enable accurate, simultaneous prediction and discrimination.

[0095] (1) Model architecture

[0096] In certain embodiments, the modNet muldmod architecture extends the original mSAnet multiple instance learning (MIL) framework to facilitate the concurrent detection of multiple RNA modification types using a single, unified multi-ciass classification model. To achieve this, the read-level probability encoder layer is modified to output a set of probability scores corresponding to the total number of classification categories present in the training data, which includes both the unmodified class and all specific modification types of interest.

[0097] Furthermore, rather than utilizing a noisy-or pooling layer, commonly employed to model the assumption that at least one read within a site exhibits a modification, the system Instead incorporates a mean pooling layer. This pooling mechanism aggregates the individual read-level probability scores into a composite site-level modification probability, thereby enabling a more generalized representation of site modification status across multiple classes.

[0098] Following the definition from Hendra et al, each position i corresponds to a reference k-mer sequence Stof length k - 5 with:

[0099] With st e {A, C, G, I / } corresponding to the unmodified nucleotide at position i. For each position i, the modification status is defined as with the index 0 representing unmodified bases yi 0= {lifpositioniisunmodified, Ootherwise} and with M representing all possible M base modifications included in the model: yt.m={Ufpositlonihasmodificationm, Ootherwise}

[0100] Here it is assumed that each position is either unmodified or has a single modification.

[0101] Similar to the single modification model in mSAnet, it is assumed that each read j at position t has an (unknown) modification status Ytj- {yij,0,yi,j,i, --->yij,M} with

[0102] Vi,j,o = {lifpositioniisunmodified, Ootherwise} and yi.i.m=[lifpositionihasmodificationm, Ootherwise} for m = 1, ... , M , with m representing all possible Mbase modifications included in the model, Here it is assumed that yt,j,m represents a single base from a single molecule.

[0103] Similar to the original mSAnet model, each read is described by a feature vector Xi j 6 7?15:

[0104] With p, a, I representing the original features obtained from modNet dat&prep corresponding to the normalized mean, standard deviation, and dwell time of read j at position i.

[0105] To encode all possible 592 5-mers Sj_1?Stand Sj+1with a fixed nucleotide at position i (i.e. st= A In the case of A modifications, st= C in the case of C modifications, s(= G in the case of G modifications, and s; = t / in the case of U modifications) into a two-dimensional vector, Ns- 592 in a mapping function f: Nsi?2. / (Sf)1and f(St)2which correspond to the first and second dimensions of the embedded vector of the 5-mer motif S) .

[0106] Each position is then described by Ntreads, with

[0107] Here, Nt- 20, where 20 reads are randomly sampled from all the reads. During inference, reads are sampled multiple times (default = 5), where the final site level modification probability is calculated as the average from the probabilities of each sampling iterations.

[0108] Firstly, modNet multimod estimates the modification probabilities for each read j at position i to have the modification m e {0,1,..., M} with m = 0 representing an unmodified base: where F-. R15RM+1Is parameterized by a neural network with two hidden layers of dimension 150 and 32 respectively, as described for mSAnet, with F(xi>y-)mrepresenting the mthdimension. The normalized modification probabilities were then calculated in accordance with equation Pitj = {Pi,j,o>Pij,i>- ■ ■ >Pt,j,M] each read using the softmax activation function to ensure that P£,7,mrepresents a probability distribution (i.e. = 1)- modNet multimod then summarizes the read-level probabilities using a mean pooling layer to estimate the site-level modification probabilities P£= {p£,0, Pi,i> ■ ■ ■ < PI,M)w>th m = 0 representing the unmodified base, and m = 1...M representing all modifications: far unmodified representing all possible M base modifications included in the model

[0109] (2) Estimation of model parameters (training)

[0110] Similar to m6Anetfwe train the network by minimizing the average cross entropy loss L between pi,mand yi,mfor all sites j“l,..JV and all modifications ents the unmodified base:

[0111] Here, / ‘and Fare updated end-to-end to minimize the cross-entropy loss L osing Adam optimizer.

[0112] Training and evaluation of &e / ?$M m©d / f / ©st / ©© detect / ©© m©de / s

[0113] (1) The generic madNet models

[0114] In accordance with the disclosed embodiments, the modNet models are specifically trained on data derived from known sequence contexts and are further evaluated for their capacity to generalize to previously unseen sequence contexts, including those arising from regions with low sequencing coverage.

[0115] To generate the training and evaluation datasets, synthetic sequences were obtained from two sources: (I) the curlcake library and (ii) a human transcriptome fragment library. For each synthetic sequence in these libraries, the dataset at each nucleotide position was divided such that approximately seventy-five percent (75%) of the reads were allocated to a training dataset, and the remaining twenty-five percent (25%) of the reads were allocated to a testing dataset. This allocation strategy ensures data exclusivity between training and testing phases, such that no Individual read is present in both datasets. With respect to the single-modification implementations of the modNet models, the labeling process was conducted as follows: all nucleotide sites derived from unmodified RNA samples were assigned a binary label of y, = 0 (indicating an unmodified state), and all nucleotide sites derived from modified R.NA samples that contained the targeted modified base were assigned a binary label of y^ - ldicating a modified state). The labeled data from training sequences containing the base of interest were compiled to form the training dataset, while the labeled data from test sequences containing the same base were complied to form the testing dataset.

[0116] To ensure data quality and consistency across model training and evaluation, specific filtering criteria were applied, hlucleotide sites that were supported by fewer than twenty (20) sequencing reads were excluded from both the training and testing datasets. Furthermore, any k-mer subsequences that contained more than one modified base were also excluded from the final datasets. This filtering process was employed to mitigate noise and prevent model bias introduced by low-read or ambiguously modified sequence regions.

[0117] The resulting curated datasets, comprising exclusively high-confidence, single- modification-labeled sites, were utilized for the training and evaluation of the modNet models, thereby enabling the development of accurate and generalizable RNA modification detection algorithms.

[0118] For the modNet multi-modification models, all sites from the unmodified sample were labeled as unmodified (y ) with 0 corresponding to the unmodified base and m representing all possible M base modifications included in the model. All sites from the modified samples with the modified bases were labeled as modified (yi / m= then combined all unmodified and modified sites from the training sequences with the base of interest into the training data set, and all unmodified and modified sites from the test sequences with the base of interest into the test data set, We used only sites with a coverage of at least 20 reads in all modified samples and the unmodified sample for training and testing.

[0119] In one aspect, the present invention discloses a computer-implemented method for detecting RNA modifications in nucleotide sequence data. The method comprises training one or more first modification detection models using synthetic RNA sequence datasets that include both unmodified sequences and sequences. These datasets are used to generate supervised labels enabling the models to learn distinguishing features associated with individual modification types. The method further comprises training one or more second modification detection models using synthetic RNA sequence datasets. Once trained, the first and second detection models are applied to an unseen RNA sequence dataset to identify potential modification events. The method includes classifying each detected modification according to its specific RNA modification type (e.g., m6A, ip, m5C) and determining, for each putative modification, a confidence probability that a nucleotide at a given position in each read of the sequence data is chemically modified.

[0120] Model training

[0121] (1) The generic modNet models

[0122] A total of nine (9) generic single-modification models were trained using the modNet training framework, each corresponding to a distinct RNA modification type selected from the group consisting of U*, mlM^, mo5U, s2U, m5C, f5C, m6AfmlA, and thG. In addition, three (3) generic multi-modification models were trained using the modNet training framework. These included: (i) a uridine- modification model capable of distinguishing between unmodified uridine and modified forms thereof, including tp, mW, mo5U, and s2U; (II) a cytidine- modification model configured to classify unmodified cytidine, m5C, and f5C; and (Hi) an adenosine-modification model configured to classify unmodified adenosine, m6A, and mlA. For the generic modNet models, reads from each site in the training set were split into training and validation (for model parameter optimization), and testing (for model selection) groups at a ratio of 1 : 1: 1, corresponding to the ratios 1 : 1 (training:validation) and 2: 1 (training+validation testing). For sites supported by more than 80 reads in the training data, we equally split the reads from each site in the training, validation, and testing groups. For sites supported by less than 80 reads, we selected 20 reads for the validation group, using all remaining reads for the training group, excluding the site from the testing group (as the minimum threshold of 20 reads is not met). This ensures that each site in the testing and validation groups is seen during training to optimise the models for the sequence context, where each site has been trained on a different set of reads from the same site during training.

[0123] Each model is trained independently using a fixed learning rate of 0.0001 over thirty (30) training epochs. The training procedure is initialized with a seed set to a value of twenty-five (25) to ensure reproducibility of results. Each training instance is further configured to execute over five (5) rounds of iterative training. The training process is carried out using a command-line interface provided by the modNet framework, executed with the following parameters: - ~lr 0.0001, —seed 25, —epochs 30, — save__perjepoch 1, and — numjterations 5.

[0124] (1) The generic modNet models

[0125] All modNet models were run on the test set using the rnodnet inference command, with the — modeLstate__dict argument used to specify each model respectively. For the single modification models, ail sites in the test data that contain the modification of interest were labeled as 1, and sites from the unmodified sample were labeled as 0 (corresponding to Yi). The multimodification models were evaluated separately for each modification, where only sites with the specific modification were labeled as 1, whereas sites with any other modification and unmodified sites were labeled as 0 (corresponding The site level probability was then used to rank all sites, and the ROC curve, Precision Recall (PR) curve, and the Area under the ROC and PR curve were calculated using the python packages scikit-learn and matplotlib.

[0126] Methods for specific analysis

[0127] Comparison to existing methods

[0128] CHEUI 4 was executed on the R.NA002 sequencing datasets corresponding to m6A, m5C, and unmodified RIMA samples in accordance with its GitHub repository. Specifically, the cheui__predict„modell.py script was employed, with the -m parameter set to the appropriate file path of either the m6A or m5C detection model, depending on the target RNA modification. The modification sites identified by CHEUI were then mapped to the corresponding modification sites predicted by the modNet models within their respective test datasets for m6A and m5C, To assess the comparative performance of the two methods, the modification probability scores generated by the modNet models were compared against the site-level probability scores reported by CHEUI.

[0129] The dorado basecaller was executed on RNA004 sequencing datasets corresponding to m6A, MJ, and unmodified RNA samples using the rna0G4...130bps.. S!jp@v5.0.0 model. The basecalling was performed with alignment parameters set to -k 5, — mm2-preset map-ont, and —modified -bases specified as either m6A or pseU, depending on the modification type being analyzed (m6A or MJ, respectively). Following basecalling, modification detection outputs were obtained using ModKit's pileup function, which processes the ModBAM files produced by dorado. The resulting bedMethyl files contain a "percent modified" value for each identified modification site, representing the proportion of reads supporting a given modification at that position. Modification sites reported by Dorado were then matched to the corresponding sites identified by the modNet models for m6A, MJ, adenosine (A), and uridine (U) in their respective test datasets. A comparative analysis was conducted between the modification probabilities predicted by the modNet models and the percent modified values reported by dorado.

[0130] When data with just a single base modification is used for evaluation, the modNet m6A and m5C models outperform CHEUI's m6A and m5C models (ROC AUC: m6A: 0.968, mSC: 0,951, vs. ROC AUC: m6A: 0.855, mSC: 0.866; Figure 7e), although the dorado m6A, m5Cfand models show higher accuracy than the modNet models (ROC AUC: m6A: 0.996, mSC: 0.987, IP: 0.999, Figure 7e). However, when other base modifications were present in the test data, they were frequently misclassified by dorado and the single modification modNet models, leading to a high number of false positives and reduced accuracy (Figures 7f, 7g, Figures 5e, 5f). In contrast, the multi-modification modNet models better maintain their accuracy in the presence of other base modifications, outperforming dorado in this scenario (Figures 5e, 5f).

[0131] Evaluation on human data

[0132] A single-modification modNet model specific to m6A was trained using all sequencing reads obtained from the mSA-modified and unmodified RNA samples. Additionally, a multi-modification modNet model targeting adenosine modifications was trained using all reads derived from mlA, m6A, and unmodified RIMA samples. It is noted that these models were trained on the combined training and validation datasets originating from a predefined training-validation-test data split, whereas models used solely for evaluation purposes were restricted to training on the training-validation subset. Inference was subsequently conducted on a HEK293T cell line RNA004 sample obtained from the Singapore Nanopore-Expression Project 60 (SG-NEx), using three models: (I) m6Anet, (ii) the single-modification m6A modNet model, and (Hi) the multi-modification adenoslne-modificatlon modNet model. Model performance was assessed through the generation of Receiver Operating Characteristic (ROC) curves. These evaluations were performed across all DRACH sequence motifs, as well as on a per-kmer basis within the DRACH motif class. The ROC curves and corresponding Area Under the Curve (AUC) metrics were computed using the Python libraries scikit-leam and matplotlib.

[0133] Read threshold optimisation

[0134] Similar to m6Anet, modNet's modification-rate estimation is based on the modification ratio, the estimated percentage of modified reads at a given site. Reads are defined as modified at any site if the read-level modification probability is above a specified threshold. As read probabilities are not calibrated, this threshold differs depending on the training data and model. To facilitate the estimation of modification rates, an optimised threshold is therefore estimated for each model. For the single-modification modNet models, we defined the 90th percentile of all read-level probabilities across all reads as the read-level probability threshold of modNet-^, -ml% and -mo5U, and the 80th percentile for modNet-s2U.

[0135] In contrast to single-modification models, which generate a read-level probability for a single predefined modification class, multi-modification models output read-level probabilities across all modification classes encompassed by the model architecture. Accordingly, for each modification class, an independent optimization of the read-level probability threshold was performed to enable accurate detection and quantification. For the uridine (U)-modificatlon modNet model, the following percentile values of read-level probabilities were defined as thresholds for their respective modification types: the 75th percentile for MJ, the 90th percentile for mlNP, the 60th percentile for mo5U, and the 25th percentile for s2U. These thresholds were derived from the empirical distribution of read-level probabilities for each modification class.

[0136] In one embodiment, the present invention provides a system for RNA chemical modification detection, comprising at least one processor, a memory communicatively coupled to the processor, and a detection module comprising a plurality of RNA modification detection models. The system includes at least one model trained to detect a single RNA modification and at least one model trained to detect multiple RIMA modifications, wherein each model is trained using labelled synthetic RNA sequence datasets comprising unmodified and chemically modified sequences with known modification types and positions. The memory stores computer-executable instructions which, when executed by the processor, cause the system to apply the plurality of RNA modification detection models to an unseen RNA sequence dataset to detect the presence of one or more RNA modifications, including individual and / or multiple modification types. The system is further configured to categorize each detected modification into a specific RNA modification type and to compute, for each detected modification event, a confidence probability that a nucleotide at a defined position within each read of the unseen dataset is chemically modified. This enables robust, probabilistic identification and classification of RNA modifications in high-throughput sequencing data.

[0137] It will be appreciated that many further modifications and permutations of various aspects of the described embodiments are possible. Accordingly, the described aspects are intended to embrace all such alterations, modifications, and variations that fall within the spirit and scope of the appended claims.

[0138] Throughout this specification and the claims which follow, unless the context requires otherwise, the word "comprise", and variations such as "comprises" and "comprising", will be understood to imply the inclusion of a stated integer or step or group of integers or steps but not the exclusion of any other integer or step or group of integers or steps.

[0139] The reference in this specification to any prior publication (or information derived from it), or to any matter which is known, is not, and should not be taken as an acknowledgment or admission or any form of suggestion that that prior publication (or information derived from it) or known matter forms part of the common general knowledge in the field of endeavour to which this specification relates.

Claims

Claims1. A method for ribonucleic acid (RIMA) modification detection, the method comprising: subjecting modified and unmodified synthetic RNA samples to direct RNA sequencing, thereby producing an RNA sequence data set comprising basecalled reads and raw electrical signal-level data; aligning the RNA sequence data set to reference synthetic RNA sequences to obtain sequence alignment reads indexed to specific sites in the reference synthetic RNA sequences; segmenting the raw electrical signal-level data into discrete events corresponding to nucleotide bases; creating an index that links each base in the sequence alignment reads with the segmented raw electrical signal-level data; aligning the segmented raw electrical signal-level data to the sequence alignment reads to produce features of the raw electrical signal-level data corresponding to specific nucleotide positions in the reference synthetic RNA sequence; labelling the modified and unmodified synthetic RNA samples using labels, based on the features at the specific nucleotide positions; and outputting labelled datasets, each comprising all RNA sequences having a particular said modification determined by the labels.

2. A method for training a detection model for detecting RNA modifications, the method comprising: producing labelled datasets according to the method of claim 1; for each labelled dataset, generating first summarised data for all sites of the respective base modification; for each unmodified RNA sequence, generating second summarised data for all sites of the respective unmodified RNAsequence corresponding to the sites of the respective base modification, wherein the first summarised data and second summarised data comprise the features; and separately inputting the first summarised data and the second summarised data into a supervised machine learning based RIMA modification detection model to train the detection model to predict presence of a modification at each site corresponding to the base.

3. The method according to claim 2, wherein the detection model includes one or more single modification detection models, each trained using ones of said datasets corresponding to a respectively different single said base modification, and the unmodified RNA sequences.

4. The method according to claim 2 or 3, wherein detection model includes one or more multiple modification detection models, each trained using ones of said datasets corresponding to said modified sequences each having more than one said base modification, and the unmodified RNA sequences.

5. The method according to claim 1, wherein the RNA modifications type comprises one or more of Pseudouridine (pll or MJ), 1-methylpseudouridine (mlpll or mli ), 5-methoxyuridine (mo5U), 2-thiouridine (s2U), 5- methylcytidine (m5C), 5-formylcytidine (f5C), N6-methyladenosine (m6A), l-methyladenosine (mlA), and thienoguanosine (thG) and unmodified nucleotides.

6. The method according to claims 1 or 2, wherein the features of the raw electrical signal-level data include one or more normalized mean ( ), standard deviation (a), and dwell time (I) for the current and adjacent k- mers.

7. The method according to claim 1, wherein the direct RNA sequencing method is nanopore sequencing.

8. A computer-implemented method for RINA modification detection, comprising: training one or more first detection model using unmodified and modified synthetic sequences to detect a single modification, using labelled datasets outputted from the method described above; training one or more second detection model using unmodified and modified synthetic sequences to detect multiple RNA modifications, using labelled datasets outputted from the method described above; subjecting an unseen RNA sequence data set to the first and second detection models, to detect one or more single RNA modifications or multiple RNA modifications in the unseen RNA sequence data set (detected RNA modifications); categorizing the detected RNA modifications into at least one specific RNA modification type; and generating, for each detected RNA modification, confidence probability that a base at a specific nucleotide position corresponding to the respective detected RNA modification in each read in the unseen RNA sequence set has been modified.

9. A system for RNA modification detection, comprising: memory; at least one processor(s); and a pre-processing module, wherein the memory stores instructions that, when executed by the processor(s), cause the pre-processing module to perform the method of any one of claims 1 - 7.

10. A system for RNA modification detection, comprising: memory; at least one processor(s); anda detection module comprising a plurality of RINA modification detection models, at least one said model being trained to detect a single RNA modification and at least one said model being trained to detect multiple RNA modifications, wherein the memory stores instructions that, when executed by the processor(s), cause the processor(s) to: subject an unseen RNA sequence data set to the plurality of RNA modification detection models, to detect presence of a single modification or multiple modifications in the unseen RNA sequence data set (detected RNA modifications); categorize the detected RNA modifications into at least one specific RNA modification type; and generate, for each detected RNA modification, a confidence probability that a base at a specific nucleotide position corresponding to the respective detected RNA modification in each read in the unseen RNA sequence set has been modified.

Images