Method for precise analysis of differential RNA editing sites based on signal-enhanced pretreatment

Abstract

Provided is a method for detecting differential RNA editing sites. The method includes the following steps: performing high-throughput sequencing and alignment on a sample set including two or more treatment conditions to obtain DNA and RNA sequencing data and alignment data; performing preliminary joint mutation detection on the DNA and RNA alignment data to obtain candidate RNA editing sites; with the candidate RNA editing sites acting as references, performing formal joint mutation detection on the DNA and RNA alignment data to obtain RNA editing sites; performing depth counting on the RNA editing sites to form a depth table, and evaluating the degree of differential editing of each editing site across different treatment conditions by using a statistical method, thereby selecting statistically significant RNA editing sites. The present method enables the high-precision and high-specificity identification of differential RNA editing sites across multiple types of RNA editing, which facilitates the elucidation of RNA editing mechanisms and the discovery of potential biomarkers, and assists in the formulation of diagnostic and therapeutic strategies for related diseases.

Description

A Precise Analysis Method for Differential RNA Editing Sites Based on Signal Enhancement Preprocessing Technical Field

[0001] This disclosure pertains to the fields of bioinformatics and biotechnology, and specifically relates to a method comprising a bioinformatics analysis process for identifying differential RNA editing events, which can identify C>U or A>I type RNA editing with high precision and specificity. Background Technology

[0002] RNA editing is an important post-transcriptional modification mechanism that regulates gene expression and function by altering the base sequence of RNA molecules. This phenomenon is particularly prominent in A-to-I (adenosine to inosine) editing catalyzed by adenosine deaminases (ADAR family) and C-to-U (cytidine to uridine) editing catalyzed by cytidine deaminases (APOBEC family). However, RNA editing events are typically low in abundance and scattered in distribution, and due to factors such as single nucleotide variants (SNVs) at the DNA level, accurate detection remains challenging despite the availability of high-throughput RNA sequencing (RNA-seq) technology.

[0003] Differentially variable RNA (DVR) editing sites are RNA single nucleotide sites that exhibit significant differences in the degree of editing between different experimental conditions or sample groups. DVR sites must meet the requirement of statistically significant differences in the degree of editing under different conditions to ensure high reliability and accuracy in detection, thus being considered high-confidence RNA editing sites.

[0004] However, the detection of differential RNA editing sites using existing technologies still faces the following major challenges: First, sequencing noise interference and false positives are problems. Existing models (such as rMATS-DVR or JACUSA2) have insufficient noise filtering capabilities, leading to high false positive rates. Second, joint DNA and RNA mutation detection has not been effectively implemented. Existing detection tools (such as VaDiR) attempt to filter DNA mutation interference using whole-genome sequencing (WGS) data from the same sample, but these methods often struggle to retain the true RNA editing signal in subsequent correction steps. Furthermore, the detection capability for complex editing types (such as C>U) is limited. Compared to the more common A-to-I editing, low-abundance editing types such as C-to-U often have weaker signals in sequencing background noise, and existing technologies have significant limitations in accurately identifying these types of edits.

[0005] In addition, base quality correction (BQSR) is typically performed before mutation detection. This step is to correct systematic biases that may be introduced during sequencing, thereby improving the accuracy of variant detection. However, traditional BQSR methods usually rely on known reference variant sites for correction. For unknown editing sites, especially RNA editing sites with limited understanding, conventional BQSR methods may not be able to effectively identify and correct these newly emerging variants. If BQSR is performed solely based on known reference sites, these unknown editing sites may be misclassified as sequencing noise, thus missing genuine editing events.

[0006] Therefore, there is an urgent need in this field for a method to detect RNA editing sites that is highly sensitive, highly specific, and covers a wider range of editing types. Summary of the Invention

[0007] The purpose of this disclosure is to provide a method for detecting differentially edited RNA sites based on signal optimization.

[0008] In one aspect, a method for detecting differentially edited RNA sites is provided, comprising the steps of:

[0009] A) Provide an independent sample set containing N samples of differentially edited RNA sites to be detected; where N is a positive integer ≥2;

[0010] Each of the independent sample sets includes: (i) RNA alignment data for each sample, which is obtained by aligning the RNA sequencing data of a single sample with a reference genome, denoted as the first dataset; and (ii) DNA alignment data for each sample, which is obtained by aligning the DNA sequencing data of a single sample with a reference genome, and after necessary sequence information preprocessing and base quality scoring correction, denoted as the second dataset.

[0011] B) Joint mutation detection was performed on the first and second datasets for each sample containing a differentially edited RNA site to be detected, and this dataset was designated as the third dataset.

[0012] The third dataset of each of the samples is merged to obtain a dataset of known RNA editing sites, denoted as the fourth dataset;

[0013] C) Using the fourth dataset as a reference input, perform base quality score correction on the first dataset of all samples; collect the first dataset after base quality score correction from all samples, retain the independent sample information, and denote it as the fifth dataset;

[0014] Collect the second dataset from all sample sets, retain the independent sample information, and denote it as the sixth dataset;

[0015] D) Perform joint mutation detection on the fifth dataset and the sixth dataset to obtain a candidate RNA editing site dataset, denoted as the seventh dataset;

[0016] E) From the first dataset of each sample containing a differential RNA editing site to be detected, obtain the allele depth of each candidate RNA editing site in the seventh dataset, thereby obtaining an allele depth table for each sample containing candidate RNA editing site information and editing site depth data, denoted as the eighth dataset;

[0017] F) Merge the eighth dataset of each sample to obtain the total allele depth table, denoted as the ninth dataset; perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset to determine the differential RNA editing sites.

[0018] In some implementations, the independent sample set includes: two or more samples with different experimental treatments, samples of normal tissue and tumor tissue, and two or more samples of time series.

[0019] In some implementations, the independent sample set contains two samples.

[0020] In some implementations, the independent sample set includes samples of normal tissue and tumor tissue.

[0021] In some implementations, the samples are duplicated.

[0022] In some implementations, the replication is selected from the group consisting of biological replication, technical replication, or a combination thereof.

[0023] In some implementations, the replication is a biological replication.

[0024] In some implementations, the number of biological replicates is 4.

[0025] In some implementations, the step of obtaining the sequencing data includes:

[0026] (s1) Extract nucleic acids, which include DNA and RNA;

[0027] (s2) Construct sequencing libraries, which include DNA sequencing libraries and RNA sequencing libraries;

[0028] (s3) Sequencing is performed to obtain the sequencing data.

[0029] In some implementation schemes, in step (s1), nucleic acids are extracted using a Roche MagNA Pure 96 instrument and its accompanying reagents.

[0030] In some implementations, the A260 / A280 of the DNA is 1.8-2.0.

[0031] In some implementations, the RIN value of the RNA is >8.

[0032] In some implementations, the DNA sequencing library is a whole-genome DNA sequencing library.

[0033] In some implementations, the sequencing depth of the DNA sequencing library is ≥30×.

[0034] In some implementations, the RNA sequencing library is constructed using operations selected from the group consisting of:

[0035] (1) Construct a chain-specific library;

[0036] (2) Reduce the proportion of rRNA;

[0037] Or a combination thereof.

[0038] In some implementations, a chain-specific library is constructed using a first conjugation method.

[0039] In some implementations, the proportion of rRNA is reduced by methods selected from the group consisting of: ribosomal RNA removal, Poly(A) enrichment, or a combination thereof.

[0040] In some implementations, the RNA sequencing library has a read length of 100 bp at both ends.

[0041] In some implementations, the sequencing depth of the RNA sequencing library is 60 million reads per sample.

[0042] In some implementations, the sequencing data is filtered.

[0043] In some implementations, the filtering operation is selected from the group consisting of:

[0044] (I) Remove the connector sequence;

[0045] (II) Removal of low-quality bases;

[0046] Or a combination thereof.

[0047] In some implementations, the software that performs the filtering includes cutadapt.

[0048] In some implementations, step (A) specifically includes:

[0049] (A1) Provides RNA sequencing data from N samples and DNA sequencing data from said N samples, where N is a positive integer ≥2;

[0050] (A2) Align the RNA sequencing data from the N samples with the reference genome to obtain N RNA alignment data and correct them; and

[0051] (A3) The DNA sequencing data from the N samples are compared with the reference genome to obtain N DNA alignment data and then corrected.

[0052] Steps (A2) and (A3) can be interchanged, performed sequentially, or performed simultaneously.

[0053] In some implementations, the software used in step (A2) to compare the RNA sequencing data with a reference genome includes STAR and HISAT2, preferably STAR.

[0054] In some implementations, in step (A2), the software used to align the RNA sequencing data with a reference genome is STAR.

[0055] In some implementations, in step (A2), the RNA sequencing data is aligned with a reference genome using a two-pass alignment method in STAR.

[0056] In some implementations, the software used in step (A3) to align the DNA sequencing data with a reference genome includes BWA.

[0057] In some implementations, in step (A3), the DNA sequencing data is compared with a reference genome using the MEM algorithm in BWA.

[0058] In some implementations, the DNA sequencing data is whole-genome sequencing data.

[0059] In some implementations, the tools used to perform the correction include Picard Tools.

[0060] In some implementations, step (A2) includes the following steps:

[0061] (A2a) Reorder the sequences according to the reference genome;

[0062] (A2b) Add sample set labels to distinguish different sample sources;

[0063] (A2c) labeling of repetitive sequences generated by PCR;

[0064] (A2d) segmentation of sequencing fragments containing variable splicing sites (including N fragments);

[0065] (A2e) Locating the target region that requires local sequence re-alignment;

[0066] (A2f) performs local re-sequence alignment for insertion or deletion (indel) mutations.

[0067] In some implementations, step (A2) further includes the step of chain-specific splitting.

[0068] In some implementations, step (A3) includes the following steps:

[0069] (A3a) Reorder the sequences according to the reference genome;

[0070] (A3b) Add sample set labels to distinguish different sample sources;

[0071] (A3c) labeling of repetitive sequences generated by PCR;

[0072] (A3d) Locating the target region that requires local sequence re-alignment;

[0073] (A3e) performs local re-sequence alignment for insertion or deletion (indel) mutations.

[0074] In some implementations, the reference genome is a human reference genome.

[0075] In some implementations, the human reference genome is GRCh38.

[0076] In some implementations, in step (B), joint mutation detection is performed using the first dataset of each sample as the experimental group and the second dataset of each sample as the control group.

[0077] In some implementations, in step (B), the software that performs the joint mutation detection includes GATK.

[0078] In some implementations, in step (B), the combined mutation detection is performed using Mutect2 from GATK.

[0079] In some implementations, step (B) further includes the step of:

[0080] (B1) Quality control, obtaining data that passes quality control;

[0081] (B2) Screen for valid variants to obtain the dataset of the known RNA editing sites;

[0082] (B3) Data indexing: Obtain an indexed dataset of the known RNA editing sites.

[0083] In some implementations, in step (B1), the software performing the quality control includes GATK.

[0084] In some implementations, the quality control is performed in step (B1) using FilterMutectCalls from GATK.

[0085] In some implementations, in step (B1), the maximum number of mutation events allowed in each assembly region (max-events-in-region) is 4.

[0086] In some implementations, in step (B1), the data passed through quality control is labeled.

[0087] In some implementations, the software used to perform the screening in step (B2) includes bcftools.

[0088] In some implementations, valid variants are screened using the labels.

[0089] In some implementations, the known RNA editing site dataset is an RNA mutation site dataset Z.

[0090] In some implementations, the software that performs the data indexing in step (B3) includes GATK.

[0091] In some implementations, in step (B3), the data indexing is performed using IndexFeatureFile from GATK.

[0092] In some implementations, the software used to perform the base quality correction in step (C) includes GATK.

[0093] In some implementations, in step (C), the base quality correction is performed using the BaseRecalibrator and ApplyBQSR from GATK.

[0094] In some implementations, in step (D), joint mutation detection is performed using the fifth dataset as the experimental group and the sixth dataset as the control group.

[0095] In some implementations, in step (D), the software that performs the joint mutation detection includes GATK.

[0096] In some implementations, in step (D), the joint mutation detection is performed using Mutect2 from GATK.

[0097] In some implementations, step (D) further includes the step:

[0098] (D1) Mutation site filtering;

[0099] (D2) Screening for single nucleotide variant sites;

[0100] (D3) Data processing.

[0101] In some implementations, the software used to perform the mutation site filtering in step (D1) includes GATK.

[0102] In some implementations, in step (D1), the mutation site filtering is performed using FilterMutectCalls from GATK.

[0103] In some implementations, in step (D1), the maximum number of mutation events allowed in each assembly region is 4.

[0104] In some implementations, the software used to perform the single nucleotide variant site screening in step (D2) includes GATK.

[0105] In some implementations, in step (D2), the single nucleotide variant sites are screened using SelectVariants in GATK.

[0106] In some implementations, step (D3) further includes the following step:

[0107] (D3a) Homopolymer nucleotide sequence filtering;

[0108] (D3b) Filtering of genomic repetitive sequences.

[0109] In some implementations, the software used to perform the homopolymer nucleotide sequence filtering in step (D3a) includes SNPIR.

[0110] In some implementations, in step (D3a), the homopolymer nucleotide sequence is filtered using filter_homopolymer_nucleotides.pl in SNPIR.

[0111] In some implementations, the software that performs the genomic repetitive sequence filtering in step (D3b) includes: SNPIR.

[0112] In some implementations, in step (D3b), the genomic repetitive sequence filtering is performed using pblat_candidates_ln.pl in SNPIR.

[0113] In some implementations, step (E) specifically includes:

[0114] (E1) Split the first dataset from N samples into positive and negative chains to obtain N positive chain subsets and N negative chain subsets;

[0115] (E2) Calculate the editing site depth data for each candidate RNA editing site in the seventh dataset in two subsets.

[0116] In some implementations, the software that performs step (E2) includes SAMtools.

[0117] In some implementations, step (E2) is performed using mpileup from SAMtools.

[0118] In some implementations, the candidate RNA editing site information includes: the chromosome where the candidate RNA editing site is located, the base position of the candidate RNA editing site, the strand where the candidate RNA editing site is located, the reference genome base type (REF base), and the edited base type (ALT base).

[0119] In some implementations, the edit site depth data includes: reference base depth, edited base depth, and base editing ratio.

[0120] In some implementations, the statistical analysis is selected from the group consisting of: model-based differential analysis, multiple hypothesis testing, or a combination thereof.

[0121] In some implementations, the model is selected from the group consisting of: generalized linear mixture models, Dirichlet-multinomial distribution models, or combinations thereof.

[0122] In some implementations, the model is a combination of a generalized linear mixture model and a Dirichlet-multinomial distribution model.

[0123] In some implementations, the multiple hypothesis testing is selected from the group consisting of: Bonferroni correction or Benjamini-Hochberg FDR correction.

[0124] In some implementations, the software used to perform the statistical analysis in step (F) includes rMATS-DVR.

[0125] In some implementations, the threshold for FDR correction is ≤0.05.

[0126] In some implementations, step (F) further includes the step of annotating the differentially edited RNA sites.

[0127] In another aspect, an apparatus or system is provided for detecting differentially edited RNA sites, the apparatus or system comprising:

[0128] (M1) Input module, which is configured to input RNA alignment data and DNA alignment data of N samples to be tested;

[0129] (M2) Scanning module, the scanning module is configured to perform the following operations: perform joint mutation detection on the RNA alignment data and DNA alignment data of each of the test samples to obtain a known RNA editing site dataset; take the known RNA editing site dataset as input, and perform base quality scoring correction on the RNA alignment data of each of the test samples to obtain corrected RNA alignment data of each of the test samples;

[0130] (M3) Detection module, configured to perform the following operations: collect corrected RNA alignment data for each test sample to obtain a fifth dataset; collect DNA alignment data for each test sample to obtain a sixth dataset; perform joint mutation detection on the fifth and sixth datasets to obtain a seventh dataset; obtain the allele depth of each candidate RNA editing site in the seventh dataset from the RNA alignment data of each test sample to obtain an eighth dataset; merge the eighth datasets of each test sample to obtain a ninth dataset; and perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset to determine differentially expressed RNA editing sites.

[0131] (M4) Output module, which is configured to output information about the differentially edited RNA sites.

[0132] In some implementations, N is a positive integer ≥2.

[0133] In some implementations, the scanning module (M2) includes the following sub-modules:

[0134] (M2.1) Preliminary detection submodule, which is configured to perform the following operations: perform joint mutation detection on the RNA alignment data and DNA alignment data of the sample for each differential RNA editing site to be detected, thereby obtaining a preliminary RNA editing site dataset;

[0135] (M2.2) Quality control submodule, which is configured to perform the following operations: perform quality control on the preliminary RNA editing site dataset to obtain a quality-controlled RNA editing site dataset;

[0136] (M2.3) Effective variant screening submodule, which is configured to perform the following operation: screen by labels in a quality-controlled RNA editing site dataset to obtain a known RNA editing site dataset;

[0137] (M2.4) Data Indexing Submodule, the submodule being configured to perform the following operations: index the known RNA editing site dataset to obtain the indexed known RNA editing site dataset;

[0138] (M2.5) Quality Correction Submodule, which is configured to perform the following operations: take the indexed known RNA editing site dataset as input, perform base quality scoring correction on the RNA sequencing data, thereby obtaining corrected RNA sequencing data.

[0139] In some implementations, the detection module (M3) includes the following sub-modules:

[0140] (M3.1) Formal detection submodule, which is configured to perform the following operations: perform joint mutation detection on the fifth dataset and the sixth dataset to obtain a formal RNA editing site dataset;

[0141] (M3.2) Mutation site filtering submodule, the submodule being configured to perform the following operation: perform quality filtering on the formal RNA editing site dataset based on mutation sites, thereby obtaining a filtered RNA editing site dataset;

[0142] (M3.3) Screening submodule, the submodule being configured to perform the following operation: extract single nucleotide variant sites from the filtered RNA editing site dataset to obtain the filtered RNA editing site dataset;

[0143] (M3.4) Data processing submodule, the submodule is configured to perform the following operations: filter the RNA editing site dataset that has passed the screening to obtain a candidate RNA editing site dataset, denoted as the seventh dataset;

[0144] (M3.5) Depth Sampling Submodule, configured to perform the following operations: obtain the allele depth of each candidate RNA editing site in the seventh dataset from the RNA alignment data of each of the test samples, thereby obtaining an allele depth table containing candidate RNA editing site information and editing site depth data, denoted as the eighth dataset; merge the eighth dataset of each of the test samples to obtain the ninth dataset.

[0145] (M3.6) Statistical Analysis Submodule, which is configured to perform the following operations: perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset to determine differential RNA editing sites.

[0146] In some implementations, within the data processing submodule (M3.4), the filtering includes:

[0147] (m1) Homopolymer nucleotide sequence filtering;

[0148] (m2) Filtering of genomic repetitive sequences.

[0149] In some implementations, the software that performs the homopolymer nucleotide sequence filtering in step (m1) includes SNPIR.

[0150] In some implementations, in step (m1), the homopolymer nucleotide sequence is filtered using filter_homopolymer_nucleotides.pl in SNPIR.

[0151] In some implementations, the software that performs the genomic repetitive sequence filtering in step (m2) includes: SNPIR.

[0152] In some implementations, in step (m2), the genomic repetitive sequence filtering is performed using pblat_candidates_ln.pl in SNPIR.

[0153] It should be understood that, within the scope of this disclosure, the above-described technical features and the technical features specifically described below (such as in the embodiments) can be combined with each other to form new or preferred technical solutions. Due to space limitations, they will not be described in detail here. Attached Figure Description

[0154] Figure 1 shows a schematic diagram of the specific steps in the analysis process of this disclosure.

[0155] Figure 2 shows the baseline results of differential RNA editing site (DVR) detection based on simulated data.

[0156] Figure 3 shows the expression levels of different genes induced by doxycycline treatment in cells for 72 hours, relative to TPB.

[0157] Figure 4 shows a sequence logo plot of nucleotide frequencies around the C>U differential RNA editing site (DVR) in A3B-mediated T-47D cells detected using different methods, using A3B-mediated DNA editing site data obtained from the GSE193225 dataset.

[0158] Figure 5 shows the sequence logo diagram of nucleotide frequencies around the C>U differential RNA editing site (DVR) in A3B-mediated SK-OV-3 cells, detected using different methods.

[0159] Figure 6 shows a sequence logo plot of nucleotide frequencies around the A>G(I) differential RNA editing site (DVR) in A3B-mediated T-47D cells, detected using different methods, using data from 10,000 A>I site editing sites sampled from REDIportal.

[0160] Figure 7 shows the sequence logo diagram of nucleotide frequencies around the A>G(I) differential RNA editing site (DVR) in A3B-mediated SK-OV-3 cells, detected using different methods.

[0161] Figure 8 shows the minimum folding energy (MFE) density plot of C>U and A>G(I) differential RNA editing sites (DVR) in A3B-mediated T-47D cells detected using different methods, based on RNAFold analysis.

[0162] Figure 9 shows the minimum folding energy (MFE) density plot of C>U and A>G(I) differential RNA editing sites (DVR) in A3B-mediated SK-OV-3 cells detected using different methods, based on RNAFold analysis.

[0163] Figure 10 shows the enhanced crosslinking and immunoprecipitation sequencing (eCLIP-seq) signal enrichment of A3B-mediated C>U differential RNA editing sites (DVRs) in T-47D cells detected using different methods, indicating the binding of A3B around the identified DVRs. Each row in the figure represents a specific C>U DVR, ordered by standard.

[0164] Figure 11 shows the enhanced crosslinking and immunoprecipitation sequencing (eCLIP-seq) signal profiles of the C>U differential RNA editing site (DVR) in A3B-mediated T-47D cells detected using different methods. The vertical dashed line represents the mean, and the area inside the horizontal dashed line represents the 95% confidence interval of the signal.

[0165] Figure 12 shows the coverage of A3B-mediated C>U differential RNA editing sites (DVR) in T-47D cells, the sequence logo diagram of surrounding nucleotide frequencies, and the signal profile of enhanced crosslinking and immunoprecipitation sequencing (eCLIP-seq) detected using the method described in this invention and the JACUSA2 RDD-RRD combined detection method. The vertical dashed line represents the mean, and the area inside the horizontal dashed line represents the 95% confidence interval of the signal. Detailed Implementation

[0166] Through extensive and in-depth research, the inventors have proposed a novel method for detecting RNA editing sites that is highly sensitive, highly specific, and covers a wide range of editing types. Specifically, the inventors have innovatively introduced a "signal optimization pre-scanning" step into the DNA / RNA joint analysis workflow, optimizing the signal-to-noise ratio of the entire analysis process. This step identifies and labels potential RNA editing sites before base quality correction, ensuring that these sites are not misjudged as sequencing noise during the correction process. These editing sites are then used as references in subsequent detection steps, thereby preserving and enhancing the true RNA editing signal. Experimental results show that, compared with existing detection methods, this method significantly improves the sensitivity and specificity for detecting editing events. Based on this, the present invention was completed.

[0167] The specific implementation steps of this disclosure are as follows:

[0168] 1. Sample preparation, nucleic acid extraction and sequencing library construction.

[0169] The experimental design disclosed herein includes two or more sample sets, which can be divided according to different experimental treatments, normal / tumor tissue controls, and time series data. Each sample set should contain at least two replicate types, which can be biological replicates or technical replicates. In a preferred embodiment, the number of replicates is four, all of which are biological replicates.

[0170] Regarding the selection of nucleic acid extraction methods, conventional genomic DNA or whole RNA extraction methods can be used, ensuring that the extracted products meet certain purity requirements. In a preferred embodiment, nucleic acid extraction is performed using a Roche MagnaPure96 instrument and its matching reagents. The A260 / A280 value of the obtained genomic DNA sample should be maintained between 1.8 and 2.0, and the RIN value of the RNA sample should be greater than 8.

[0171] For the construction of whole-genome sequencing libraries, standard whole-genome DNA library construction methods compatible with the sequencing instrument can be selected. For RNA sequencing library construction, conventional construction methods matching the sequencing instrument used can also be chosen. This disclosure particularly recommends using strand-specific RNA sequencing libraries, as this type of library helps to clarify the directionality of RNA sequencing. The preferred construction method is the first-strand method, i.e., first-strand synthesis. During the construction of RNA sequencing libraries, to reduce the proportion of rRNA in the library, ribosome removal or poly-A enrichment methods can be selected. In the preferred embodiment of this disclosure, ribosome removal is used. After library construction, single-end sequencing or paired-end sequencing can be performed, selected according to experimental requirements. In the preferred embodiment, paired-end sequencing is used to improve the accuracy of sequence alignment and the depth of downstream analysis.

[0172] 2. Sequencing

[0173] In this disclosure, sequencing is performed using a sequencing platform compatible with the sequencing library. Preferred sequencing instruments include BGISeq, MGISeq, and Illumina HiSeq to ensure data quality and reliability. The sequencing depth should be selected within a reasonable range to meet the needs of downstream analysis. In a preferred embodiment, for human cell samples, the whole-genome sequencing depth is 30X or higher, and the RNA sequencing depth is 60 million reads per sample to ensure sufficient coverage and statistical significance. The sequencing read length used should be no less than 100 bp. In a preferred embodiment, if paired-end sequencing is used, the read length is 100 bp.

[0174] 3. Comparison and preprocessing of RNA / DNA combined mutation detection data

[0175] The main objective is to process the raw sequencing files to generate sequence alignment files. In a preferred embodiment of this disclosure, the raw data includes high-throughput sequencing data of control group RNA, control group DNA, treated group RNA, and treated group DNA, in FASTQ format. The sequence filtering step involves removing adapter sequences and low-quality bases from the sequencing data. In a preferred embodiment of this disclosure, the software used is cutadapt, which is used to efficiently filter low-quality data and ensure that the data quality meets the requirements of subsequent analysis. Next, the RNA sequence is aligned to a reference genome using STAR software, preferably a two-pass alignment method (STAR ​​2-pass) to improve alignment accuracy and the ability to detect alternative splicing events. For DNA sequence alignment, Burrows-Wheeler Aligner (BWA) software is used, specifically the BWA MEM algorithm, to ensure high-quality DNA sequence alignment.

[0176] RNA data preprocessing primarily aims to correct RNA sequence alignment results to meet the needs of subsequent mutation detection. This work can be accomplished using the Picard Tools suite, including: reordering sequences according to a reference genome; adding sample set markers to distinguish different sample sources; labeling PCR-generated repetitive sequences to reduce false positive rates; segmenting sequencing fragments containing alternative splicing sites (including N-fragments) to improve alignment accuracy; locating target regions requiring local sequence realignment; and performing local realignment of insertion or deletion mutations to ensure accurate mutation location. Base quality correction is not involved in RNA data preprocessing.

[0177] DNA data preprocessing: The main purpose is to perform necessary corrections to the DNA sequence alignment results to meet the needs of subsequent mutation detection. This part of the work can also be done using Picard Tools, including: reordering sequences according to the reference genome; adding sample set markers to distinguish the origin of different samples; marking repetitive sequences generated by PCR to reduce false positives caused by PCR amplification; locating all target regions that need to be re-aligned locally; and performing local re-alignment of insertion or deletion mutations to ensure the accuracy of mutation sites.

[0178] Preliminary base quality correction of DNA data: Base quality score correction (BQSR) is performed on the sequencing data processed above to further improve data quality.

[0179] Sequencing alignment files split by strand: This step aims to split RNA sequencing alignment files strand-specifically, dividing each sample into positive strand (relative to the reference genome) and negative strand alignment files for more accurate downstream analysis.

[0180] 4. Signal optimization pre-scan and base quality correction

[0181] The primary purpose of signal optimization pre-scanning is to accurately identify and label potential RNA editing sites before base quality correction, thereby reducing the impact of sequencing noise and non-editing mutations on subsequent analyses. This step ensures that these potential editing sites are not erroneously degraded during BQSR, thus improving the accuracy of RNA editing events.

[0182] In the signal optimization pre-scan, this disclosure preferably uses the Mutect2 function in the GATK (version 4 and above) tool for joint DNA / RNA mutation detection. During this process, RNA sequencing data is treated as a "tumor" sample, while genomic DNA sequencing data serves as a "normal" control. This joint mutation detection strategy effectively distinguishes RNA editing events from DNA mutations, thereby optimizing the identification and detection results of RNA editing sites.

[0183] The specific process is as follows:

[0184] (1) Mutation detection: The GATK Mutect2 tool was run on each sample individually to identify preliminary mutation sites in the RNA data.

[0185] (2) Quality Control: The detected mutation sites were quality controlled using GATK's FilterMutectCalls tool. In this step, the maximum number of allowed mutation events per assembly region was set to 4 to prevent the false filtering of true mutation sites.

[0186] (3) Screening for valid variants: Use bcftools to screen and retain single nucleotide variant sites that pass the quality control "PASS" tag, generating the RNA mutation site dataset Z.

[0187] (4) Data indexing: Use GATK's IndexFeatureFile tool to create an index for dataset Z to support subsequent BQSR steps.

[0188] After completing the signal optimization pre-scan, this disclosure uses the BaseRecalibrator and ApplyBQSR modules in the GATK toolchain to perform base quality correction (BQSR) on the RNA sequencing data. BQSR is an important step in ensuring the accuracy of RNA editing detection, eliminating sequencing errors, and reducing pseudo-mutations.

[0189] In this step, this disclosure uses the list of known variant sites generated in the "signal optimization pre-scan" stage as one of the inputs to BQSR. These known sites are specially processed during the BQSR process to reduce the impact of sequencing noise on RNA editing sites, thereby improving the reliability of editing sites.

[0190] 5. Formal DNA / RNA Combined Mutation Detection

[0191] Following the signal optimization pre-scan and base quality correction (BQSR) steps, this disclosure further employs a DNA / RNA combined mutation detection step to more accurately identify and differentiate RNA editing sites from DNA mutations. This step, by combining RNA sequencing data with genomic DNA sequencing data, further improves the detection specificity and accuracy of RNA editing events and reduces false positives due to DNA mutations.

[0192] The specific process is as follows:

[0193] (1) Combined mutation detection: GATK’s Mutect2 tool was used to perform combined mutation detection on RNA sequencing data and DNA sequencing data. RNA sequencing data was treated as “tumor” samples and DNA sequencing data was used as “normal” controls to identify the difference between RNA editing events and DNA mutations.

[0194] (2) Mutation Site Filtering: For mutation sites identified by Mutect2, the FilterMutectCalls tool is used for quality filtering. In this disclosure, the maximum number of mutation events in a region (max-events-in-region) is set to 4, thereby allowing the detection of RNA editing site clusters and preventing genuine RNA editing sites from being incorrectly filtered out. During this process, low-quality mutation data is removed, ensuring that only mutation sites with high confidence are retained. This filtering step helps reduce false positives caused by sequencing errors and low-quality data.

[0195] (3) Screening for single nucleotide variant sites

[0196] The SelectVariants tool is used to filter the variant data, extract single nucleotide variant sites, and output new, high-quality variant data files. This step ensures that the output dataset focuses on sites with high mutation frequency and high reliability, further improving the accuracy of RNA editing event detection.

[0197] (4) Subsequent data processing

[0198] The screened data undergoes further cleaning, including homopolymer nucleotide sequence filtering and genomic repetitive sequence filtering, to remove interfering information such as homologous polynucleotides that could lead to misclassification. After these refined processes, a more accurate set of mutation sites is obtained. Finally, through intersection analysis, combining RNA and DNA mutation site data, the final mutation results are generated, confirming the mutations that truly belong to RNA editing events.

[0199] 6. Editing depth sampling and differential statistical analysis

[0200] Editing Depth Detection: Sequence counting and editing depth calculation are performed on candidate RNA editing sites to obtain the counts of mutant and normal sequences, and the editing depth is calculated. In a preferred embodiment of this disclosure, the mpileup command in the samtools software is used for sequence counting to ensure accurate calculation of editing depth. Further, information from each sample is collected to compile a list containing detailed information and editing depth of candidate RNA editing sites. This list should at least include: the chromosome and base position of the candidate RNA editing site; the reference genome base type (REF base) and the edited base type (ALT base); the strand on which the RNA editing site is located; and the editing status of each sample in each sample set at that editing site, including the reference base depth, the edited base depth, and the proportion of edited bases to all sequences at that site (base editing ratio).

[0201] Statistical analysis of differences: The RNA editing site information and editing depth list were imported into the rMATS-DVR software package. The GLMM model was used to compare whether there were statistical differences in the editing ratio of candidate RNA editing sites in the comparison of different sample sets, and the false discovery rate (FDR) of each site was calculated to assess the significance of editing events under different conditions.

[0202] Final selection of differentially edited RNA sites (DVRs): Based on the statistical analysis results, RNA editing sites with a sequencing depth greater than 10 and an FDR of less than 5% were selected as the final confirmed DVRs.

[0203] 7. Examples of computer instructions and parameters, etc.

[0204] (1) Remove connectors and filter low-quality bases. Taking the cutadapt software as an example, the command format is as follows:

[0205] cutadapt-g ADAPTER-O 5-e 0-o sample.trimmed.fastq sample.fastq--minimum-length35--discard-untrimmed--info-file=reads.adapter.txt

[0206] cutadapt-q 10-o output.fastq input.fastq

[0207] (2) RNA sequencing alignment, using STAR software as an example, the command format is as follows:

[0208] STAR--runThreadN 20--runMode genomeGenerate--genomeDir$GENOME_DIR--genomeFastaFiles$REFERENCE.FA--sjdbGTFfile$GENCODE.GTF--sjdbOverhang 99

[0209] STAR--runThreadN 20--genomeDir$GENOME_DIR--outFileNamePrefix$PREFIX_RNA --readFilesIn 1.fq.gz 2.fq.gz--readFilesCommand zcat--outSJfilterReads Unique--outFilterMultimapNmax 1

[0210] (3) DNA sequencing and alignment, using BWA and SAMtools software as examples, the command format is as follows:

[0211] bwa mem-t 20$REFERENCE.FA read_1.fq read_2.fq|samtools view-b>GENOMIC_DNA.bam

[0212] (4) DNA / RNA sequencing BAM files are used to perform sequencing group division. Taking Picard as an example, the command format is as follows:

[0213] picard AddOrReplaceReadGroups INPUT=$label_reordered.bam OUTPUT=$label_addrg.bam RGID=$label RGLB=$label RGPL=COMPLETE RGPU=lane1RGSM=$label

[0214] (5) Reorder the input RNA / RNA sequencing BAM files. Taking the Picard tool as an example, the command format is as follows:

[0215] picard ReorderSam INPUT=$bam OUTPUT=${label}_reordered.bam S=true R=$REFERENCE.FA

[0216] (6) DNA / RNA sequence BAM files are used to mark repetitive sequences. Taking the Picard tool as an example, the command syntax is as follows:

[0217] picard MarkDuplicates INPUT=${label}_addrg.bam OUTPUT=${label}_dedup.bam CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT READ_NAME_REGEX=null METRICS_FILE=$label_metrics.txt

[0218] (7) Processing RNA sequence BAM files with long insertions or deletions, and segmenting CIGAR strings. Taking the GATK tool as an example, the command paradigm is as follows:

[0219] gatk SplitNCigarReads-R$REFERENCE.FA-I${label}_dedup.bam-O${label}_split.bam

[0220] (8) BQSR preprocessing of DNA sequencing BAM files

[0221] gatk BaseRecalibrator-I${label}_dedup.bam-R$REFERENCE.FA-O${label}_recalibration_report.grp-known-sites dbsnp.vcf

[0222] gatk ApplyBQSR-R$REFERENCE.FA-I${label}_dedup.bam-bqsr-recal-file${label}_recalibration_report.grp-O${label}_recalibration.bam

[0223] (9) Preliminary DNA / RNA combined mutation detection, using GATK software as an example, the command format is as follows:

[0224] gatk Mutect2-R$REFERENCE.FA-I$RNA_BAM-I$DNA_BAM-normal$DNA_LABEL-O output_mutect.vcf

[0225] The command paradigm for filtering preliminary DNA / RNA joint mutation detection results is as follows:

[0226] gatk FilterMutectCalls-V output.mutect.vcf-R$REFERENCE.FA-O output.filter.vcf

[0227] bcftools view-f PASS output.filter.vcf>known_variants.vcf

[0228] (10) Perform BQSR processing on the RNA sequence BAM file. The command format is as follows:

[0229] gatk BaseRecalibrator-I${label}_split.bam-R$REFERENCE.FA-O${label}recalibration_report.grp --known-sites known_variants.vcf

[0230] gatk ApplyBQSR-R$REFERENCE.FA-I${label}_split.bam--bqsr-recal-file${label}recalibration_report.grp-O${label}RNA_recalibrated.bam

[0231] (11) Perform DNA / RNA combined mutation detection again. Taking GATK software as an example, the command format is as follows:

[0232] gatk Mutect2-R$REFERENCE.FA -I$RNA_BAM-I$DNA_BAM-normal$DNA_LABEL-O output_mutect.vcf

[0233] Filtering of DNA / RNA combined mutation detection results

[0234] gatk FilterMutectCalls-V output.mutect.vcf-R$REFERENCE.FA--min-median-base-quality 12--max-events-in-region 4-O output.filter.vcf

[0235] bcftools view-f PASS output.filter.vcf>{label}.vcf

[0236] Extracting single-base variations

[0237] gatk SelectVariants-R$REFERENCE.FA -V${label}.vcf --select-type-to-include SNP-O${label}.SNV.vcf

[0238] (12) Filtering of homologous polynucleotides and repetitive sequence regions, using the Perl script in the SNPIR package as an example, the command paradigm is as follows:

[0239] perl filter_homopolymer_nucleotides.pl-infile${label}.SNV.vcf-outfile{label}_homo.vcf-refgenome$REFERENCE.FA

[0240] (13) Prepare the BAM file for pblat analysis. Generate a reference file by merging all RNA-seq BAM files, and sort and index them. The command paradigm is as follows:

[0241] samtools merge-f${label}all.bam${RNA_bam[@]}&&samtools sort-@$threads-o${label}all_sorted.bam${label}all.bam&&samtools index${label}all_sorted.bam

[0242] (14) Use pblat for alignment and filter duplicate alignment sequences. Use the Perl script in the SNPIR package to process candidate variants and further generate high-confidence variant files:

[0243] perl${script_dir}pblat_candidates_ln.pl-infile${label}.homo.vcf-outfile${label}.pblat.vcf-bamfile${label}all_sorted.bam-refgenome$REFERENCE.FA-minbasequal 5

[0244] awk'{OFS="\t";$2=$2-1"\t"$2;print$0}'${label}.pblat.vcf>${label}.knownedit.vcf&&bedtools intersect-a${label}.SNV.vcf-b${label}.knownedit.vcf-wa-header>${label}.final.vcf

[0245] (15) Edit depth sampling to obtain an allele depth table. Taking samtools as an example, the command paradigm is as follows:

[0246] samtools mpileup-Bd 100000-f$REFERENCE.FA -l final.vcf-q 30-Q 17-a-oresults.pileup${RNA_bam[@]}

[0247] (16) Statistical analysis of DVR: The rMATS-DVR package script was used to analyze the mutation depth, count the RNA editing events with significant differences, and the FDR.py script was used to calculate the false detection rate (FDR). The command paradigm is as follows:

[0248] python vcf_to_mats_input_For_Mutect2.py$label.final.vcf${label}.inc.txt$RNA_bam_sample1$RNA_bam_sample2 20 5 T${label}.pileup T

[0249] python MATS_LRT.py${label}.inc.txt${label}_rMATS-DVR_results 4 0.0001

[0250] python FDR.py${label}_rMATS-DVR_results_rMATS_Result_P.txt${label}_rMATS-DVR_results_rMATS_Result_FDR.txt

[0251] (17) Annotate and summarize the DVR:

[0252] python snv_annotation.py --input${label}_rMATS-DVR_results_rMATS_Result_FDR.txt --output${label}_rMATS-DVR_Result.txt --summary${label}rMATS-DVR_Result_summary.txt --label1 DMSO --label2 DOX--snp$knownSNV--repeat NA --editing$knownediting--gene$geneanno

[0253] 8. Implementation of differential RNA editing site detection methods described in previous studies: Data preprocessing

[0254] Sequencing data from the same sources as those described in the preceding embodiments of this disclosure were selected, including: whole-genome sequencing data; DNA information for providing a reference genome; and RNA sequencing data under appropriate conditions, including an induced group and a control group, each with at least three biological replicates. Routine quality control and adapter removal were performed on the raw WGS and RNA-seq data (Fastq format) to ensure data quality met the requirements for subsequent analysis. Sequence alignment was then performed using the same version of the human genome (e.g., GRCh38) and corresponding gene annotation information (e.g., GENCODE GRCh38.p13).

[0255] 9. Development of previously described methods for detecting differentially expressed RNA editing sites: JACUSA method

[0256] Data input and alignment file preparation: For RNA-seq alignment files, use the Picard tool to mark duplicates and retain sequencing / alignment quality information. Following the JACUSA developer's recommendations, prepare the RNA sequencing alignment file (BAM format) to ensure that samples can be identified as different conditions or groups (i.e., induced group vs. control group).

[0257] RDD and RRD pattern detection: RDD (RNA-DNA Difference) mode: The RNA alignment file is compared with WGS or a reference genome to detect differences between RNA and DNA sequences. RRD (RNA-RNA Difference) mode: Only the RNA sequencing results of the induced group and the control group are compared to detect sites that differ under different treatment conditions. In this embodiment, to ensure consistency with the "double comparison" procedure of this disclosure, "joint analysis of RDD and RRD modes" is also performed, that is, the results of the two modes are merged or their intersection is taken to obtain a candidate set of differentially edited sites.

[0258] Filtering and Result Output: JACUSA's built-in variant filtering criteria (such as sequencing coverage, quality score, and homology region determination) are used to obtain preliminary candidates. Further analysis of allele frequency differences among candidate variants in different groups is conducted, and the sites with significant differences are output, i.e., the DVRs identified by JACUSA.

[0259] 10. Development of previously described differential RNA editing site detection methods: rMATS-DVR method

[0260] Alignment file preparation: Generate an RNA alignment file (BAM) using an RNA-seq alignment method consistent with this disclosure (such as STAR alignment). As recommended by the developers, no additional DNA correction is introduced; if WGS information is used, it is typically only used for routine BQSR or filtering.

[0261] RNA variant detection: As recommended by the developers, HaplotypeCaller in the GATK software package was used to perform variant detection on the RNA data to obtain an initial set of candidate RNA variants. Because rMATS-DVR does not have the DNA-RNA joint mutation detection function of the method described in this disclosure, some DNA SNVs may be mistakenly identified as RNA variants.

[0262] Differential analysis: Based on the obtained candidate sites, the depth of variant alleles in RNA sequencing of the induced group and the control group was compared, and the original GLMM statistical model of rMATS was used for differential detection. If a variant site showed a significant difference in editing frequency between the two treatment conditions and passed the internal filtering threshold, it was identified as a DVR and output.

[0263] 11. Development of previously described methods for detecting differentially expressed RNA editing sites: VaDiR method

[0264] Single-condition RNA variation detection: VaDiR is designed with RNA-DNA variation comparison in mind, but it often only supports single-sample comparisons or is not specifically designed for multi-repetitive omics data, making it unable to directly perform joint statistical modeling of multiple biological replicates. In this embodiment, to be as consistent as possible with this disclosure, RNA variation detection was first performed on the induced group and the control group, and obviously low-quality or suspected repetitive regions were excluded.

[0265] Differential selection and determination: Among the sites detected in the induction group and the control group, variant sites that appear frequently and meet the internal filtering criteria of VaDiR (such as sequencing coverage ≥10, sequencing quality score ≥20, etc.) are screened. The allele frequencies between the two groups are compared. If the frequency difference reaches the specified threshold, it is determined to be a differentially edited site.

[0266] It should be understood that the specific methods and experimental conditions described below with varying degrees of detail are intended to provide a substantive understanding of this disclosure. Definitions of certain terms used in this specification are provided below. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this disclosure pertains.

[0267] the term

[0268] As used herein, the terms “containing” or “including (comprise)” can be open-ended, semi-closed, or closed-ended. In other words, the terms also include “consistently made of” or “made of”.

[0269] As used herein, the term “and / or” refers to and covers any and all possible combinations of one or more of the related listed items.

[0270] Cytidine deaminase

[0271] As used in this article, the term "cytidine deaminase" refers to a class of enzymes that can deaminoknose cytosine (or cytidine) to uracil (or uridine). These enzymes not only play a crucial role in regulating genomic mutations and diversity at the DNA level, but also mediate C>U nucleotide base changes during RNA editing, affecting gene transcription and translation, thus playing an important role in various physiological and pathological processes (such as tumorigenesis).

[0272] RDD

[0273] As used in this article, the term "RDD" (RNA-DNA Difference) refers to the identification of sites at the RNA level where the genomic DNA (gDNA) and RNA-complementary DNA (cDNA) sequences differ at the DNA level. This stage of analysis aims to eliminate false positives caused by single nucleotide variants (SNVs) and extract variations truly resulting from RNA editing.

[0274] RRD

[0275] As used in this article, the term "RRD" (RNA-RNA Difference) refers to the identification of differentially edited RNA sites that exhibit significant changes in editing levels between two conditions by statistically analyzing the differences in editing depth at variant sites through comparison of RNA sequencing (RNA-seq) data from different experimental conditions or biological states. The analysis in the RRD phase emphasizes the impact of variability in experimental conditions on RNA editing events.

[0276] TPR

[0277] As used in this article, the term "TPR" (True Positive Rate) refers to the ability of a model or method to detect real RNA editing events, defined as the proportion of detected true positive events to the total number of actual positive events.

[0278] Wherein: TP (True Positives) represents the number of detected true positive events; FN (False Negatives) represents the number of missed true positive events. TPR reflects the sensitivity of the detection method and is an important indicator of the ability to detect edit sites.

[0279] Precision

[0280] As used in this article, the term "precision" refers to the proportion of true positive events among all RNA editing events detected by a model or method.

[0281] Where: TP (True Positives) represents the number of true positive events detected; FP (False Positives) represents the number of false positive events detected. Precision is an important parameter for evaluating the reliability of a detection method; higher precision indicates fewer false positive events and more reliable results.

[0282] Accuracy

[0283] As used in this article, the term "accuracy" refers to the overall accuracy of a model or method in classifying all edit events (including positive and negative), defined as the proportion of correctly classified events to all classified events.

[0284] Wherein: TP (True Positives) represents the number of true positive events detected; TN (True Negatives) represents the number of true negative events detected; FP (False Positives) represents the number of events falsely detected as positive; and FN (False Negatives) represents the number of true positive events missed. Accuracy comprehensively measures the sensitivity and specificity of the method.

[0285] F-score

[0286] As used in this article, the term "F-score" is an indicator that combines Precision and TPR, defined as the harmonic mean of Precision and TPR.

[0287] Among them, Precision and TPR are as described above. The F-score is used to balance the sensitivity and reliability of the detection method and is an important evaluation of the overall performance of the model.

[0288] APOBEC3B(A3B)

[0289] As used herein, the term "APOBEC3B (A3B)" is a member of the cytidine deaminase family, whose functions include deamination of cytosine residues in DNA or RNA molecules. APOBEC3B plays an important role in viral defense and tumorigenesis; its aberrant expression can lead to genomic instability and the accumulation of numerous mutations, thereby promoting tumorigenesis and development. Therefore, in the RNA editing research disclosed herein, APOBEC3B is used as a typical model enzyme for C>U editing events at the RNA level.

[0290] DNA / RNA Joint Variant Calling

[0291] In this disclosure, "DNA / RNA combined mutation detection" refers to the simultaneous use of DNA sequencing data and RNA sequencing data in the same sample to identify and distinguish RNA-level base variations (e.g., RNA editing) from single nucleotide variants (SNVs) or insertion / deletion variations carried by the genome itself. This process can accurately eliminate false positives caused by genomic variations, improving the specificity and accuracy of RNA editing site detection.

[0292] Single nucleotide variant (SNV)

[0293] As used in this article, the term "SNV" refers to a single nucleotide change in a genome sequence that differs from a reference sequence by only a single base. SNVs can occur in coding, non-coding, or regulatory regions and may cause changes in protein coding (such as missense mutations) or functional effects.

[0294] RNA editing site

[0295] As used herein, the term "RNA editing site" refers to a base sequence change that is inconsistent with the template DNA sequence but occurs in the RNA transcript. The main types include adenosine (A) to inosine (I) conversions and cytosine (C) to uracil (U) conversions. This disclosure focuses on the high-precision identification and validation of these editing sites using a combined DNA / RNA strategy and statistical analysis.

[0296] Differentially Edited RNA Site (DVR)

[0297] As used herein, the term "differentiated RNA editing site" refers to an editing site where the level of RNA editing (base substitution ratio) differs significantly between two or more groups, under different experimental conditions, or at different time points. This disclosure quantifies and compares the editing depth of each sample using generalized linear models or other robust statistical methods to screen for RNA editing sites that are significantly affected by experimental or environmental conditions, thereby establishing the "differentiated RNA editing sites (DVR)" described in this disclosure.

[0298] Signal optimization pre-scan

[0299] This disclosure introduces a signal optimization pre-scan in the DNA / RNA joint analysis process to improve the signal-to-noise ratio of the process and enhance the accuracy and specificity of detection.

[0300] In signal optimization pre-scanning, multi-level joint analysis and pre-labeling strategies significantly improved signal preservation and false positive filtering in RNA editing site detection.

[0301] First, it uses combined DNA / RNA mutation detection to initially label candidate RNA editing sites during the initial screening stage of RNA sequencing data. This preprocessing strategy maximizes the identification of potential real RNA editing signals, preventing them from being downweighted or lost in subsequent BQSR steps.

[0302] Secondly, by jointly analyzing mutations at the DNA and RNA levels, signal optimization pre-scanning selects high-confidence RNA editing events instead of DNA editing events, reducing false positive sites caused by background noise or artifacts. This is particularly effective in precisely filtering low-complexity sequences and genomic repetitive regions, minimizing false positive interference. For complex editing types such as C-to-U, signal optimization pre-scanning, combining DNA and RNA data, effectively improves the detection efficiency of low-abundance editing signals. Highlighting these complex editing sites during the initial screening stage provides a strong foundation for subsequent differential analysis.

[0303] Base Quality Correction (BQSR)

[0304] As used in this article, the terms "base quality correction" and "BASR" are used interchangeably and are common steps before mutation detection. Base quality correction aims to correct systematic biases that may be introduced during sequencing, thereby improving the accuracy of mutation detection. When sequencers read DNA sequences, factors such as chemical reactions and instrument performance may lead to low quality scores for certain bases. These low-quality bases may be misidentified as mutations, affecting the reliability of subsequent analyses. By correcting these low-quality bases in the BQSR step, the interference of sequencing errors on mutation detection results can be effectively reduced, improving the accuracy and reliability of mutation detection. Currently, BQSR is mainly conducted based on known reference sequences, adjusting (generally lowering) the quality scores of mutated bases to reduce the interference of sequencing noise on mutation detection.

[0305] The main advantages of this disclosure include:

[0306] (1) This disclosure develops a method for detecting RNA editing sites, and for the first time introduces a signal optimization pre-scanning step in the method. This step reduces the number of bases that are misjudged as noise by recalibrating the base mass fraction, thereby significantly improving the retention rate of editing signals, and thus significantly increasing the number of RNA editing events detected, enhancing the identification of real RNA editing signals, greatly reducing the interference of background noise, and having the characteristics of high sensitivity.

[0307] (2) Compared with existing RNA editing site detection methods and combined detection methods, the method disclosed herein, by combining DNA and RNA data and using signal optimization pre-scanning, can identify and filter false positive sites caused by single nucleotide variations, thereby effectively reducing the false positive rate.

[0308] (3) This disclosure uses a generalized linear mixed statistical analysis method to perform differential analysis on RNA editing sites. Compared with traditional methods, it can significantly improve the detection rate and reliability of differential RNA editing sites and has the characteristics of high accuracy.

[0309] (4) This disclosure provides a method for studying RNA editing site detection and / or RNA editing mechanisms under different biological conditions, different experimental treatments or disease states at the whole transcriptome and whole genome levels.

[0310] The present disclosure is further illustrated below with reference to specific embodiments. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the disclosure. Experimental methods in the following embodiments, unless otherwise specified, are generally performed under conventional conditions, such as those described in Sambrook et al., Molecular Cloning: A Laboratory Manual (New York: Cold Spring Harbor Laboratory Press, 1989), or as recommended by the manufacturer. Unless otherwise stated, percentages and parts are weight percentages and parts by weight.

[0311] Materials and methods

[0312] 1. Data Acquisition and Sequence Alignment

[0313] 1.1 Sample Preparation

[0314] Sample source: RNA sequencing (RNA-seq) and whole genome sequencing (WGS) data were collected from multiple independent sample sets (N samples, N≥2). The sample sets can cover different biological conditions, experimental treatments, or disease states.

[0315] Library Construction:

[0316] RNA sequencing libraries: Strand-specific library construction methods (such as the First-Strand method) are preferred to ensure accurate acquisition of positive and negative strand information. Ribosomal RNA Depletion or Poly(A) Enrichment is used to reduce the proportion of rRNA in the library and improve mRNA coverage.

[0317] DNA sequencing library: Standard whole-genome DNA library construction methods are used to ensure sequencing depth of 30× or higher to provide high-quality genomic data.

[0318] Sequencing platform: Select a high-throughput sequencing platform (such as BGISeq-500, MGISeq or Illumina HiSeq), with the preferred RNA sequencing read length being 100bp at both ends and the sequencing depth being 60 million reads per sample.

[0319] 1.2 Sequence Alignment and Preprocessing Steps

[0320] RNA sequencing data alignment:

[0321] The RNA-seq data were aligned to a reference genome (e.g., GRCh38) and a reference transcriptome (e.g., Gencode) using STAR software (version 2-pass mode).

[0322] Preliminary adapter removal, low-quality sequence filtering, and repetitive sequence labeling are performed to generate a corrected RNA alignment file (.bam format).

[0323] DNA sequencing data alignment:

[0324] The BWA-MEM algorithm was used to align WGS data to the same reference genome.

[0325] Adapter removal, low-quality filtering, and repetitive sequence labeling are performed to generate a corrected DNA alignment file (.bam format).

[0326] 2. Signal optimization pre-scanning steps

[0327] Tool Selection: This disclosure utilizes the Mutect2 function in GATK (version 4 and above) for combined DNA / RNA mutation detection, where RNA sequencing data is treated as a "tumor" sample, while DNA sequencing data serves as a "normal" control. This combined mutation detection strategy effectively distinguishes between RNA editing sites and genomic DNA mutations, optimizing the detection results of RNA editing events.

[0328] Mutation detection process:

[0329] Gatk Mutect2 mutation detection was performed on each sample individually to generate preliminary RNA mutation sites.

[0330] The gatk FilterMutectCalls tool was used for quality control of mutation sites, with a maximum allowed number of mutation events of 4 per assembly region to prevent genuine mutations from being falsely filtered out.

[0331] Using bcftools, we retain those single nucleotide variants that pass the quality control "PASS" tag, forming the RNA mutation site dataset Z.

[0332] Use gatk IndexFeatureFile to create an index for dataset Z, preparing for the subsequent BQSR steps.

[0333] Objective: The primary purpose of signal optimization pre-scanning is to accurately identify potential RNA editing sites before BQSR correction and label these sites as known variants, ensuring that their quality fraction is not mistakenly reduced during subsequent base quality fraction correction. This step significantly reduces misjudgments caused by sequencing noise or non-editing mutations, improving the reliability and accuracy of RNA editing sites.

[0334] 3. Formal DNA / RNA Combined Mutation Detection

[0335] 3.1 BQSR Correction

[0336] Tool Selection: In RNA editing detection, BQSR correction is a crucial step in ensuring high accuracy of mutation sites and reducing sequencing noise. This disclosure uses the BaseRecalibrator and ApplyBQSR modules from the GATK toolchain to correct the quality of read data during RNA-seq data analysis, thereby eliminating false mutations caused by errors in the sequencing process.

[0337] Correction Procedure: Use the list of known variant sites generated in the "Signal Optimization Pre-scan" step as the "Known Sites" section of the BQSR reference input. Perform BQSR correction.

[0338] 3.2 DNA / RNA Combined Mutation Detection

[0339] Tool Selection: This disclosure uses the GATK Mutect2 tool for combined DNA / RNA mutation detection. By combining data from RNA sequencing and genomic DNA sequencing, it can effectively distinguish between RNA editing and DNA mutation.

[0340] process:

[0341] Combined mutation detection: Using the Mutect2 tool, RNA sequencing data is used as a "tumor" sample and DNA sequencing data is used as a "normal" control to identify RNA editing events and DNA mutations.

[0342] Filtering mutation sites: Use the FilterMutectCalls tool to filter the mutation data generated by Mutect2, retaining only mutation sites that meet the standards and excluding low-quality data.

[0343] Selecting single nucleotide variant sites: Use the SelectVariants tool to filter out single nucleotide variant sites and output new high-quality variant data files.

[0344] Subsequent processing: The screened single nucleotide variant data undergoes further filtering to remove interfering information such as homologous polynucleotides, ensuring the accuracy of mutation sites. Specifically, this includes homopolymer sequence filtering: removing mutations located in homopolymer repetitive sequences to prevent false positives due to sequencing errors; and genomic repetitive sequence filtering: removing mutations located in highly repetitive regions of the genome to reduce variations introduced by multiple alignment errors. Finally, intersection analysis of RNA and DNA mutation sites is performed to generate the final variant results.

[0345] Objective: By using combined DNA / RNA mutation detection, this disclosure effectively distinguishes between RNA editing events and DNA mutations, significantly improves the detection specificity of RNA editing sites, avoids misjudgments caused by DNA mutations, and ensures more accurate RNA editing detection.

[0346] 4. Editing depth sampling and differential statistical analysis

[0347] 4.1 Editing Depth Sampling

[0348] Tool selection: Use the SAMtools mpileup command to count the strand-specific editing depth of the selected candidate RNA editing sites.

[0349] Implementation steps:

[0350] The RNA-seq alignment file is split into positive and negative strands, generating n positive strand subsets and n negative strand subsets respectively.

[0351] For each candidate RNA editing site, the depth of the reference base and the edited base were calculated in each sample.

[0352] Calculate the editing ratio (edited base depth / total depth) for each site in each sample.

[0353] 4.2 Statistical Analysis of Differences

[0354] Statistical methods: The Generalized Linear Mixed Model (GLMM) or Dirichlet-multinomial distribution model was used to analyze the differences in RNA editing ratios among multiple biological conditions or experimental groups.

[0355] Tool selection: Use efficient statistical analysis packages such as rMATS and edgeR, combined with the R language for data processing and analysis.

[0356] Implementation steps:

[0357] The editing depth data of each candidate RNA editing site in each sample were compiled into a matrix file.

[0358] Statistical models were applied to the matrix data to assess whether there were significant differences in the editing ratio of each site under different conditions.

[0359] Perform multiple hypothesis testing (such as Benjamini-Hochberg FDR correction) to screen out DVRs with statistical significance (such as FDR≤0.05).

[0360] 4.3. Final DVR Confirmation and Output

[0361] Confirmation criteria: RNA editing sites with sequencing depth ≥10 and differential statistical analysis result FDR ≤0.05 were selected as the final differential RNA editing sites (DVR).

[0362] Output: Generate a DVR list containing the following information: chromosome location (chromosome number, base position); reference base type (REF base); edited base type (ALT base); strand (positive or negative); edit depth and edit ratio in each sample; statistical analysis results (such as T-value, FDR, etc.).

[0363] Example 1: Construction of virtual data required for evaluating the detection efficiency of differentially expressed RNA editing sites

[0364] This description evaluates the effectiveness of the disclosed method in detecting differentially edited RNA sites (DVRs) through virtual data construction and simulation analysis, and compares it with other existing DVR detection technologies. To achieve this goal, specific construction and simulation analysis were carried out based on virtual data from whole-genome sequencing (WGS) and transcriptome sequencing (RNA-seq).

[0365] In terms of data construction, this embodiment uses two data types for simulation: WGS and RNA-seq. WGS data uses a 2×150nt sequencing strategy with an average sequencing depth of 33× and two biological replicates. RNA-seq data uses a 2×100nt strand-specific sequencing strategy with a total read length of approximately 16 million and three to five biological replicates. In both datasets, 50,000 identical single nucleotide polymorphisms (SNVs) located in DNA and RNA are inserted to simulate the presence of single nucleotide polymorphisms (SNPs) in the actual genome. Furthermore, an additional 50,000 RNA editing sites existing only at the RNA level are inserted into the RNA-seq data to further simulate the real RNA editing event environment, thus providing basic data support for subsequent differential RNA editing site detection.

[0366] To evaluate the detection effectiveness of differentially modified RNA editing sites, this embodiment sets differential characteristics for the inserted RNA editing sites. By adjusting the variation frequency of RNA editing sites between the "treatment group" and "control group" experimental conditions, approximately 6,000 sites were made to meet the DVR criteria. These criteria are based on thresholds defined by methods such as rMATS GLMM or JACUSA. Sites that meet the above conditions are defined as "successfully detected" true positives, used to evaluate the sensitivity and accuracy of the method described in this disclosure in DVR detection.

[0367] To further simulate real-world applications, this embodiment also inserted approximately 2,000 SNVs that meet the requirements of the RNA-RNA detection method into the WGS and RNA-seq data to simulate false positive results. This operation aims to examine the differences in theoretical precision and accuracy between the pure RNA-RNA detection mode (without reference DNA sequencing) and the DNA-RNA combined detection mode. This design further limits the upper limits of the theoretical precision and accuracy of different detection methods. For example, the upper limits of precision and accuracy for the pure RNA-RNA detection mode are limited to approximately 0.75 and 0.12, respectively, while the upper limits for precision and accuracy for the DNA-RNA combined detection mode are approximately 0.98 and 0.56, respectively (see Table 2 for details).

[0368] Based on virtual data, this embodiment comprehensively evaluates the performance of the method described in this disclosure through a series of indicators. In the simulated data, if the detection method can detect an RNA editing site consistent with the set differential RNA editing site, it is considered positive; false detections outside the DVR definition range are considered false positives. Based on the aforementioned number of true positives and false positives, key indicators such as precision, accuracy, and sensitivity of each detection method are calculated to quantify its detection efficiency and reliability.

[0369] The virtual dataset and evaluation system established through the above steps allow for a comprehensive comparison of the differences between the detection efficiency and accuracy of the tools described in this disclosure and other mainstream DVRs. Tables 1 and 2 detail the sequencing parameters used and statistical information on the inserted RNA editing sites, providing data support for subsequent result analysis and comparison.

[0370] Table 1. Overview of Simulation Data Used in the Evaluation

[0371] Table 2. Intrinsic parameters involved in the simulation data used for evaluation

[0372] Example 2: Comparison of Methodologies Involved in Cross-sectional Evaluation

[0373] To further verify the impact of the signal optimization pre-scanning step (hereinafter referred to as "the correction step") introduced in the method described in this disclosure on DVR detection, this embodiment compares the detection effects of retaining or omitting the correction step in the CADRES workflow. Simultaneously, based on the same set of simulation data, the existing JACUSA2 was used for analysis in both RNA-RNA difference (RRD) and RNA-DNA difference (RDD) modes, and the intersection of the variants obtained from RRD and RDD was subsequently calculated (i.e., the JACUSA2 RDD-RRD joint detection method) to maximize the potential of DVR detection. Furthermore, existing methodologies such as VaDiR and rMATS-DVR were evaluated to examine the detection efficiency and applicability of each method from multiple perspectives. Table 3 below provides a comparison of the methods in this embodiment and an overview of their characteristics.

[0374] Table 3 lists the functional characteristics of each methodology used in the comparative evaluation in this embodiment, including: whether base quality fraction recorrection (BQSR) is performed, the method of comparing RNA-RNA differences (RRD) and RNA-DNA differences (RDD), the degree of support for multiple conditions or multiple samples, the biological replicate requirement, the variant detection engine used, the statistical analysis method for variant depth, and additional alignment artifact filtering measures. This table provides a clear understanding of the differences in settings at different detection stages of each analytical procedure, facilitating subsequent comparative analysis of the detection results.

[0375] Table 3. Comparison of methodologies involved in cross-sectional evaluation of simulation data in the examples

[0376] Example 3: Evaluating the effectiveness of differential locus detection methods using virtual data

[0377] This embodiment evaluates the sensitivity and specificity of different detection methods in DVR detection using a virtual data simulation experimental system, and focuses on the impact of the number of repeated experiments on detection performance.

[0378] The experimental results are shown in Figure 2. The same virtual dataset was used in the experimental design, applied to the methods described in this disclosure (with or without signal optimization pre-scanning steps), the JACUSA2 method (RDD-RRD joint detection mode), and traditional methods that only consider the differences between RNA and DNA (such as VaDiR). These methods were compared under the same conditions, and their performance was quantitatively evaluated using metrics such as True Positive Rate (TPR), Precision, Accuracy, and F1 Score. The True Positive Rate (TPR) measures the method's ability to correctly identify DVR events; Precision assesses the proportion of sites in the detection results that are true DVR events; Recall reflects the proportion of true DVR sites that the method can detect; and the F1 Score is a performance metric that comprehensively balances precision and recall.

[0379] The experimental results (Figure 2) show that all detection methods performed similarly in true positive rate (TPR), but exhibited significant differences in precision and accuracy. In particular, the method described in this disclosure, combined with the JACUSA2 RDD-RRD joint detection method, effectively reduced the number of false positive sites by integrating the analysis of RDDs and RRDs, demonstrating significantly better precision and detection reliability than traditional methods that only consider RDDs (such as VaDiR). Specifically, the precision of the method described in this disclosure is 10-15% higher than that of the JACUSA2 RDD-RRD joint detection method, indicating a significant advantage in reducing false positives. However, the recall rate of the method described in this disclosure is slightly lower than that of the JACUSA2 RDD-RRD joint detection method by 3-5%, suggesting that the method may have certain limitations in detecting certain edge editing events.

[0380] By introducing a signal optimization pre-scanning step, the method described in this disclosure significantly improves performance in terms of true positive rate (TPR) and accuracy. Signal optimization pre-scanning effectively identifies and preserves genuine RNA editing signals while reducing false positives, thus enhancing the overall sensitivity of the detection. This result fully demonstrates the crucial role of signal optimization pre-scanning in the method described in this disclosure, giving it a significant advantage in distinguishing genuine editing events from background noise. Further experiments analyzed the impact of the number of experimental repetitions on detection performance. The results showed that when the number of repetitions exceeded four, the improvement in evaluation metrics tended to plateau. This indicates that increasing the number of repetitions within a certain range can improve the stability and reliability of the detection, but the optimization effect of increasing the number of repetitions beyond a certain threshold is limited.

[0381] This embodiment verifies the high efficiency of the method described in this disclosure in the detection of differential RNA editing sites. Especially after incorporating optimized pre-scanning steps, the method outperforms other mainstream methods in terms of accuracy and true positive rate (TPR). After optimizing the number of experimental replicates, the method demonstrates high accuracy and sensitivity in RNA editing event identification. By effectively identifying statistically significant RNA editing sites, the method not only improves the reliability of DVR detection but also demonstrates its potential as an important tool in RNA editing research. These experimental results provide strong support for the scientific validity and practical application value of the method disclosed in this disclosure.

[0382] Example 4: Construction of a cell line expressing cytidine deaminase

[0383] This embodiment describes the construction of a lentivirus-mediated cell line inducibly expressing the cytidine deaminase APOBEC3B, used to study cytidine deaminase-mediated RNA editing function and to verify its feasibility as an evaluation tool for differential RNA editing site detection. The T-47D and SK-OV-3 cell lines were selected as model systems, with T-47D cells exhibiting moderate levels of endogenous A3B expression, while SK-OV-3 cells express almost no endogenous A3B. These two cell lines were chosen to compare A3B-mediated DVR events, providing a contrast between endogenous and exogenous A3B expression.

[0384] First, the gene coding sequence of APOBEC3B (NM_004900.5) was synthesized by a contractor (Shanghai Sangon Biotech Co., Ltd.), and its human cell codons were optimized to improve expression efficiency. The coding sequence was cloned into the pTRIPZ lentiviral vector (GE Healthcare) via AgeI and ClaI (BspDI) restriction sites, respectively, to obtain lentiviral expression plasmids for inducible expression of APOBEC3B and its inactive mutant control. To prepare lentiviral particles, the above expression plasmids were co-transfected with helper plasmids pMD2.G and psPAX2 into 293TN cells (Clonetech). The culture medium was changed 24 hours after transfection, and the supernatant was collected at 48 and 72 hours post-transfection. The collected supernatant was filtered through a 0.44 μm microfiltration filter, and the lentiviral particles were concentrated according to the standard procedure provided by the Peg-IT kit (System Biosciences).

[0385] T-47D and SK-OV-3 cells were transduced using the concentrated lentiviral particles described above. The culture medium was changed 48 hours after transduction to remove residual virus from uninfected cells, and 1 μg / ml puromycin was added to the medium 72 hours post-transduction for cell selection. Under continuous puromycin selection pressure, doxycycline-induced APOBEC3B-expressing T-47D and SK-OV-3 cell lines were successfully established.

[0386] Example 5: Cell sample collection, nucleic acid extraction and high-throughput sequencing

[0387] This embodiment describes the steps for collecting samples and extracting nucleic acids from constructed T-47D and SK-OV-3 cells to ensure data suitable for subsequent high-throughput sequencing. Cells were cultured under standard conditions, with passaged in RPMI-1640 medium, 10% fetal bovine serum, and 1% pen / strep, for both T-47D and SK-OV-3 cells. In the induced group cells, doxycycline at a concentration of 1 μg / mL was added for 72 hours to ensure adequate APOBEC3B expression. The uninduced control cells did not receive doxycycline.

[0388] After induction, cells were isolated from the culture plate using trypsin digestion and washed with PBS buffer to remove residual substances from the culture medium. The collected cell samples were used for RNA and genomic DNA (gDNA) extraction. RNA extraction was performed using Trizol reagent or other similar reagents, strictly following the manufacturer's instructions to ensure the integrity and purity of the RNA samples. Genomic DNA extraction was performed using the QIAamp DNA Mini Kit (Qiagen), following the manufacturer's standard operating procedures to ensure high-quality DNA samples. After extraction, quality control was performed on the RNA and DNA samples. The concentration and purity of RNA and DNA were measured using a NanoDrop instrument, ensuring that their OD260 / 280 ratio was between 1.8 and 2.0. Further electrophoretic analysis of the RNA was performed using an Agilent Bioanalyzer 2100 to ensure that the RIN value of the extracted RNA was greater than 9.0. RNA samples were then enriched for mRNA and used for library construction using ribosome removal methods. High-throughput RNA sequencing was subsequently performed on the samples using the BGI-SEQ500 platform (BGI Genomics, Shenzhen). The sequencing depth was optimized based on the experimental design (greater than 60 million reads) to ensure coverage and data accuracy.

[0389] RNA-seq results were counted using the Subread software package, and the results were normalized using the TPM method and further corrected using TBP gene expression levels. The results are shown in Figure 3. After 72 hours of doxycycline treatment, the level of exogenously induced APOBEC3B expression was significantly increased in T-47D and SK-OV-3 cells, indicating successful construction of the induced expression cells.

[0390] Example 6: Detection of RNA differential editing sites using the method described in this disclosure

[0391] This embodiment describes the detection of differentially edited RNA sites in T47D and SK-OV-3 cells using the methods described in the "Specific Implementation Methods" section of this disclosure. Statistical analysis results (Tables 4 and 5) show that after 72 hours, under 1 μg / mL doxycycline-induced APOBEC3B expression in cells, the C>U editing frequency changed significantly, and RNA editing sites were successfully detected and screened.

[0392] Table 4. Comparison of the number of differentially expressed RNA editing sites (DVRs) detected by different methods in APOBEC3B-induced T-47D cells.

[0393] Table 5. Comparison of the number of differentially expressed RNA editing sites (DVRs) detected by different methods in APOBEC3B-induced SK-OV-3 cells.

[0394] Example 7: Detection of differentially edited RNA sites using control methods (JACUSA, rMATS-DVR, VaDiR)

[0395] This embodiment aims to illustrate the specific procedure for detecting DVR based on sequencing data from T-47D and SK-OV-3 induced expression cell lines, under experimental conditions identical or similar to those described in this disclosure, to verify the improvements and advantages of this disclosure in detection effectiveness and accuracy. RNA differential editing site results were collected from sequencing data of the induced expression cell treatment groups / control groups described in Examples 4-6 using JACUSA, rMATS-DVR, and VaDiR methods, respectively. The number of different types of editing events, such as C>U and A>G(I), detected by each method is shown in Tables 4 and 5.

[0396] The data in Table 4 show that, with the introduction of the signal optimization pre-scanning step, the method disclosed herein significantly outperforms the version without signal optimization pre-scanning in terms of RNA editing event detection capability. Specifically, the signal optimization pre-scanning increased the number of detected C>U and A>G(I) type RNA editing events, with the number of C>U DVR sites increasing from 739 to 816 and the number of A>G(I) DVR sites increasing from 1451 to 3764. This improvement indicates that the signal optimization pre-scanning step can more effectively identify low-abundance RNA editing events, especially complex editing types.

[0397] Experimental results showed significant overlap between the DVRs detected by rMATS-DVR and SNVs in T-47D cells (see Table 4). This result is expected because the rMATS-DVR method does not employ a mechanism to filter SNVs. In the rMATS-DVR analysis including SNVs, 33% of C>U DVRs showed a decrease in candidate allele frequency after APOBEC3B induction, a pattern not observed in DVRs obtained using the method described in this disclosure. Further testing was conducted using the SK-OV-3 cell model, where the exogenous expression level of APOBEC3B was low but biologically relevant. The results showed that the detection rate of the method described in this disclosure was slightly lower than that of the JACUSA2 RDD-RRD combined detection method (see Tables 4 and 5). However, the overall difference was not significant. This result indicates that the method described in this disclosure significantly improves the false positive problem caused by SNVs compared to methods that only consider RNA mutations.

[0398] Example 8: Comparison of sequence features of C>U RNA differentially edited sites detected by different methods

[0399] This embodiment compares the sequence characteristics of the disclosed method and a control method in detecting C>U RNA differential editing sites, including upstream and downstream base bias and enzyme-specific motif enrichment, to verify the ability of different methods to detect differential editing sites accurately. The study of RNA editing site sequence characteristics is particularly important, as the RNA editing mechanism of APOBEC3B is highly dependent on the surrounding base sequence background. APOBEC3B has known RNA editing target sites that favor a specific sequence feature (motif), namely 5'-UUCM (where M represents A or C). By detecting whether C>U events are enriched in the 5'-UUCM motif and its upstream and downstream specific sequence environment, it is possible to more accurately verify whether the editing site is mediated by APOBEC3B, thereby eliminating false positive signals.

[0400] In the experiment, sequences within a certain base range (e.g., ±5 or ±10 nt) upstream and downstream of the DVR sites detected by different methods were first extracted, and uniform filtering and quality control were performed to ensure the consistency and reliability of all sequence data. Nucleotide spectrum analysis and Logo visualization were performed on these sequence data using tools such as WebLogo to clarify the base enrichment patterns of the DVRs detected by each method around the editing sites. The analysis results (Figures 4 and 5) show that the C>U sites detected by the method described in this disclosure significantly exhibit the enrichment characteristics of the APOBEC3B specific target site motif (5'-UUCM) in the editing core region (e.g., near the C site). In contrast, the control methods (including the JACUSA2 RDD-RRD combined detection method, JACUSA2 RDD, JACUSA2 RRD, and VaDiR method) failed to adequately filter DNA SNVs or sequencing noise, resulting in weaker sequence bias and a certain degree of deviation in the detected C>U sites. Furthermore, the sequence characteristics of the APOBEC3B differential RNA editing sites detected by the method described in this disclosure are significantly different from the sequence characteristics of DNA editing mediated by APOBEC3B (obtained from the dataset GSE193225) (Figure 4). The above sequence characteristic analysis demonstrates that the method described in this disclosure exhibits higher accuracy and effectiveness in the specific identification of RNA editing events compared to existing methods (including JACUSA2 RDD-RRD combined detection, JACUSA2 RDD, JACUSA2 RRD, rMATS-DVR, and VaDiR methods).

[0401] Example 9: Comparison of sequence features of differentially edited A>I RNA sites detected by different methods

[0402] This embodiment compares the performance of the method disclosed herein with other detection methods (JACUSA2, rMATS-DVR, VaDiR) in terms of specificity and accuracy by detecting and analyzing the A>I RNA editing site, and further verifies the advantages of the method described herein in enriching A>I RNA editing site-specific motifs (5'-YAS, Y = U or C, S = G or C). A>I RNA editing is a post-transcriptional modification event catalyzed by ADAR (adenosine deaminase acting on RNA) family enzymes, characterized by the deamination of adenosine (A) to inosine (I). This editing event is usually manifested as A>G substitution in high-throughput RNA sequencing. ADAR enzymes are highly dependent on target sequences, and their editing preferences are significantly influenced by the adjacent base environment. In particular, the 5′-YAS motif (where Y represents U or C, and S represents G or C) located in the sequences before and after the editing site is considered an important marker of ADAR activity. Detecting whether this specific motif can be enriched is key to verifying the biological authenticity and specificity of the RNA editing event.

[0403] The experiment used two independent cell model datasets, including the T-47D induced expression cell line and the SK-OV-3 cell line. The experimental results (Figures 6 and 7) show that the method described in this disclosure exhibits significant detection capability for A>I editing sites in both T-47D and SK-OV-3 cells. The A>I sites detected by the method described in this disclosure are enriched with 5′-YAS motifs in both upstream and downstream sequences, a feature highly consistent with known ADAR targeting sites in the REDIportal database. In contrast, although the rMATS-DVR and JACUSA methods can detect a large number of A>I events, their failure to effectively exclude false positive signals (such as SNVs or sequencing noise) results in significantly weaker 5′-YAS motif enrichment capabilities for the detected A>I sites compared to the method described in this disclosure. Furthermore, the A>I sites detected by these control methods contain some signals interfered with by SNVs, further weakening the biological relevance of the editing sites. In further comparisons, although the sensitivity of the method described in this disclosure in the SK-OV-3 cell model was slightly lower than that of the JACUSA combined with RDD-RRD combined detection, its significant advantage in enriching A>I site motifs indicates that it is more suitable for specific detection. In particular, through the "signal optimization pre-scan" step, the method described in this disclosure can separate the real RNA editing signal from the noise, significantly reduce the interference of false positive signals, and ensure a high degree of consistency between the detection results of A>I sites and ADAR activity characteristics. Through the above analysis, this embodiment verifies the significant advantages of the method described in this disclosure in terms of specificity and accuracy in detecting A>I RNA editing events, especially in the accuracy and effectiveness of enriching ADAR-preferred motifs (5'-YAS).

[0404] Example 10: Comparison of physicochemical characteristics of detected differentially edited RNA sites

[0405] This embodiment verifies the biological relevance and accuracy of the detected sites by comparing the physicochemical characteristics of the RNA differential editing sites (i.e., RNA secondary structure energy, Minimum Folding Energy, VaDiR) with control methods (JACUSA2, rMATS-DVR, VaDiR). The experiment used high-throughput sequencing data from the same induced APOBEC3B cells (T-47D and SK-OV-3) as in the previous embodiments, and ran the methods described in this disclosure and the control methods under consistent analytical conditions to obtain the output DVR sites. Using RNA secondary structure prediction tools such as RNAFold, MFE was calculated for sequences within a certain range (±100 nt) upstream and downstream of each DVR site, and the MFE distribution characteristics of the DVRs obtained by each method were compared.

[0406] The analysis results (Figures 8 and 9) show that the DVRs detected by both the method described in this disclosure and the control method exhibit significant RNA secondary structure folding characteristics, and their MFE distribution is consistent with the characteristics of known RNA editing sites. This indicates that the detected DVRs are more likely located in regions where RNA secondary structures can form, consistent with the biological characteristics of RNA editing. Overall, this embodiment verifies that both the method described in this disclosure and the control method can identify biologically characteristic DVR sites when detecting differentially edited RNA sites, which further confirms that the detected DVRs are real events at the RNA level, rather than DNA SNVs or interference from background noise.

[0407] Example 12: Comparison of detected C>U RNA differential editing sites with cytidine deaminase binding RNA sites

[0408] This embodiment utilizes enhanced crosslinking and immunoprecipitation sequencing (eCLIP-seq) technology to verify the accuracy and specificity of the disclosed method and control method in detecting differentially edited RNA sites. eCLIP-seq is an experimental technique for studying the interaction between RNA-binding proteins (RBPs) and RNA. It involves immobilizing RNA and its binding proteins into a complex through ultraviolet crosslinking, followed by immunoprecipitation to enrich specific protein-nucleic acid complexes. The RNA is then isolated, reverse transcribed, and analyzed using high-throughput sequencing to determine the binding sites of the RNA-binding proteins and the distribution of target RNA. Applied to this study, this technique specifically captures the APOBEC3B-RNA complex through immunoprecipitation, and further sequences the bound RNA fragments to obtain the APOBEC3B binding site on RNA.

[0409] In this study, publicly available eCLIP-seq data (GSE193225) was used to analyze the APOBEC3B protein binding signal around DVRs identified by different detection methods (including the method described in this disclosure, JACUSA2 RDD-RRD combined detection, rMATS-DVR, and VaDiR method). The results showed that the C>U DVR site detected by the method described in this disclosure was significantly better than other methods in terms of eCLIP-seq signal enrichment. Specifically, the DVR detected by the method described in this disclosure exhibited strong APOBEC3B binding signal enrichment in the upstream region of the editing site, which is consistent with the biological mechanism by which APOBEC3B functions by binding RNA during C>U RNA editing (Figure 10).

[0410] In contrast, while methods such as JACUSA2 RDD-RRD combined detection and rMATS-DVR also detected a certain number of C>U DVRs, their binding signal enrichment was low, and the consistency with the APOBEC3B binding pattern was poor. The statistically significant difference in signal enrichment compared to the method described in this disclosure (Figure 11) indicates that the method described in this disclosure can not only identify RNA editing events but also accurately capture signals of RNA-binding proteins associated with editing events. In particular, it has a significant advantage over other methods in capturing the binding signal of the APOBEC3B-RNA complex.

[0411] Example 13: Differences in RNA differential editing sites detected by the method described in this disclosure compared to the best previous methods.

[0412] This embodiment verifies the improvements and technical advantages of the disclosed method in detecting differentially edited RNA (DVR) by directly comparing it with the currently recognized superior method for detecting DVR (JACUSA2RDD-RRD combined detection). Experimental data were obtained from high-throughput sequencing of induced APOBEC3B cells (T-47D and SK-OV-3), and the same quality control and alignment strategies were used for analysis to ensure the fairness of the comparison. Intersection, union, and difference analyses were performed on the DVR sets detected by CADRES and JACUSA2RDD-RRD combined detection to further evaluate the specificity and biological significance of overlapping and unique sites of the two methods in detecting editing types such as C>U and A>G(I).

[0413] The results (Figure 12) show that, compared with the JACUSA2 RDD-RRD combined detection, the method described in this disclosure exhibits higher specificity and a lower false positive rate, especially in C>U type editing events, where it can more effectively exclude interference sites originating from sequencing noise. Analysis of overlapping sites revealed that the common DVRs detected by the method described in this disclosure and JACUSA2 were mostly sites with significant changes in editing rate and high coverage. These sites were highly consistent with the sequence characteristics of known APOBEC3B RNA editing sites, and the eCLIP-seq signal also showed significant enrichment. For the RNA differential editing sites unique to the method described in this disclosure, these sites showed stronger enrichment in A3B-specific sequence preferences (such as 5′-UUCM), and the protein signals near the binding sites were more significant, further validating that these sites are A3B-mediated genuine editing events. In contrast, the RNA differential editing sites unique to the JACUSA2 RDD-RRD combined detection no longer possessed the 5′-UUCM sequence characteristics, but instead exhibited the 5′-UCA sequence characteristics. This sequence signature is consistent with the DNA editing site sequence signature (5'-TCA) of APOBEC3B, and some sites are associated with genomic variation background or low coverage regions. Furthermore, the eCLIP-seq signal did not show significant enrichment, suggesting that the results may contain more false positive signals caused by SNVs.

[0414] By comparing with current best-in-class RNA editing detection protocols, this embodiment demonstrates that the method described in this disclosure exhibits higher specificity and accuracy in detecting real-world DVRs, particularly showing a significant advantage in identifying C>U type RNA editing events. This method not only provides more reliable basic data for research on RNA editing mechanisms but also demonstrates greater practical value for subsequent applications in biological function analysis and clinical biomarker screening.

[0415] All documents mentioned in this disclosure are incorporated herein by reference as if each document were individually incorporated herein by reference. Furthermore, it should be understood that after reading the foregoing teachings of this disclosure, those skilled in the art can make various alterations or modifications to this disclosure, and these equivalent forms also fall within the scope defined by the appended claims.

Claims

1. A method for detecting differentially edited RNA sites, characterized in that, Including the following steps: A) Provide an independent sample set containing N samples of differentially edited RNA sites to be detected; where N is a positive integer ≥2; Each of the independent sample sets includes: (i) RNA alignment data for each sample, which is obtained by aligning the RNA sequencing data of a single sample with a reference genome, denoted as the first dataset; and (ii) DNA alignment data for each sample, which is obtained by aligning the DNA sequencing data of a single sample with a reference genome, and after necessary sequence information preprocessing and base quality scoring correction, denoted as the second dataset. B) Joint mutation detection was performed on the first and second datasets for each sample containing a differentially edited RNA site to be detected, and this dataset was designated as the third dataset. The third dataset of each of the samples is merged to obtain a dataset of known RNA editing sites, denoted as the fourth dataset; C) Using the fourth dataset as a reference input, perform base quality score correction on the first dataset of all samples; collect the first dataset after base quality score correction from all samples, retain the independent sample information, and denote it as the fifth dataset; Collect the second dataset from all sample sets, retain the independent sample information, and denote it as the sixth dataset; D) Perform joint mutation detection on the fifth dataset and the sixth dataset to obtain a candidate RNA editing site dataset, denoted as the seventh dataset; E) From the first dataset of each sample containing a differential RNA editing site to be detected, obtain the allele depth of each candidate RNA editing site in the seventh dataset, thereby obtaining an allele depth table for each sample containing candidate RNA editing site information and editing site depth data, denoted as the eighth dataset; F) Merge the eighth dataset of each sample to obtain the total allele depth table, denoted as the ninth dataset; perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset to determine the differential RNA editing sites.

2. The method as described in claim 1, characterized in that, Step (A) specifically includes: (A1) Provides RNA sequencing data from N samples and DNA sequencing data from said N samples, where N is a positive integer ≥2; (A2) Align the RNA sequencing data from the N samples with the reference genome to obtain N RNA alignment data and correct them; and (A3) The DNA sequencing data from the N samples are compared with the reference genome to obtain N DNA alignment data and then corrected. Steps (A2) and (A3) can be interchanged, performed sequentially, or performed simultaneously.

3. The method as described in claim 1 or 2, characterized in that, Step (B) also includes the following steps: (B1) Quality control, obtaining data that passes quality control; (B2) Screen for valid variants to obtain the dataset of the known RNA editing sites; (B3) Data indexing: Obtain an indexed dataset of the known RNA editing sites.

4. The method as described in claim 3, characterized in that, In step (B1), the maximum number of mutation events allowed in each assembly region (max-events-in-region) is 4.

5. The method according to any one of claims 1-3, characterized in that, Step (D) also includes the following step: (D1) Mutation site filtering; (D2) Screening for single nucleotide variant sites; (D3) Data processing.

6. The method as described in claim 5, characterized in that, Step (D3) also includes the following steps: (D3a) Homopolymer nucleotide sequence filtering; (D3b) Filtering of genomic repetitive sequences.

7. The method according to any one of claims 1-3 or 5, characterized in that, Step (E) specifically includes: (E1) Split the first dataset from N samples into positive and negative chains to obtain N positive chain subsets and N negative chain subsets; (E2) Calculate the editing site depth data for each candidate RNA editing site in the seventh dataset in two subsets.

8. The method according to any one of claims 1-3 or 5, characterized in that, The candidate RNA editing site information includes: the chromosome where the candidate RNA editing site is located, the base position of the candidate RNA editing site, the strand where the candidate RNA editing site is located, the base type of the reference genome, and the base type after editing.

9. The method according to any one of claims 1-3 or 5, characterized in that, The edit site depth data includes: reference base depth, edited base depth, and base editing ratio.

10. An apparatus or system for detecting differentially edited RNA sites, characterized in that, The device or system includes: (M1) Input module, which is configured to input RNA alignment data and DNA alignment data of N samples to be tested; (M2) Scanning module, the scanning module is configured to perform the following operations: perform joint mutation detection on the RNA alignment data and DNA alignment data of each of the test samples to obtain a known RNA editing site dataset; take the known RNA editing site dataset as input, and perform base quality scoring correction on the RNA alignment data of each of the test samples to obtain corrected RNA alignment data of each of the test samples; (M3) Detection module, configured to perform the following operations: collect corrected RNA alignment data for each test sample to obtain a fifth dataset; collect DNA alignment data for each test sample to obtain a sixth dataset; perform joint mutation detection on the fifth and sixth datasets to obtain a seventh dataset; obtain the allele depth of each candidate RNA editing site in the seventh dataset from the RNA alignment data of each test sample to obtain an eighth dataset; merge the eighth datasets of each test sample to obtain a ninth dataset; and perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset to determine differentially expressed RNA editing sites. (M4) Output module, which is configured to output information about the differentially edited RNA sites.

11. The device or system as claimed in claim 10, characterized in that, The scanning module (M2) includes the following sub-modules: (M2.1) Preliminary detection submodule, which is configured to perform the following operations: perform joint mutation detection on the RNA alignment data and DNA alignment data of the sample for each differential RNA editing site to be detected, thereby obtaining a preliminary RNA editing site dataset; (M2.2) Quality control submodule, which is configured to perform the following operations: perform quality control on the preliminary RNA editing site dataset to obtain a quality-controlled RNA editing site dataset; (M2.3) Effective variant screening submodule, which is configured to perform the following operation: screen by labels in a quality-controlled RNA editing site dataset to obtain a known RNA editing site dataset; (M2.4) Data Indexing Submodule, the submodule being configured to perform the following operations: index the known RNA editing site dataset to obtain the indexed known RNA editing site dataset; (M2.5) Quality Correction Submodule, which is configured to perform the following operations: take the indexed known RNA editing site dataset as input, perform base quality scoring correction on the RNA alignment data, thereby obtaining corrected RNA alignment data.

12. The device or system as claimed in claim 10 or 11, characterized in that, The detection module (M3) includes the following sub-modules: (M3.1) Formal detection submodule, which is configured to perform the following operations: perform joint mutation detection on the fifth dataset and the sixth dataset to obtain a formal RNA editing site dataset; (M3.2) Mutation site filtering submodule, the submodule being configured to perform the following operation: perform quality filtering on the formal RNA editing site dataset based on mutation sites, thereby obtaining a filtered RNA editing site dataset; (M3.3) Screening submodule, the submodule being configured to perform the following operation: extract single nucleotide variant sites from the filtered RNA editing site dataset to obtain the filtered RNA editing site dataset; (M3.4) Data processing submodule, the submodule is configured to perform the following operations: filter the RNA editing site dataset that has passed the screening to obtain a candidate RNA editing site dataset, denoted as the seventh dataset; (M3.5) Depth Sampling Submodule, configured to perform the following operations: obtain the allele depth of each candidate RNA editing site in the seventh dataset from the RNA alignment data of each of the test samples, thereby obtaining an allele depth table containing candidate RNA editing site information and editing site depth data, denoted as the eighth dataset; merge the eighth dataset of each of the test samples to obtain the ninth dataset. (M3.6) Statistical Analysis Submodule, which is configured to perform the following operations: perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset to determine differential RNA editing sites.

Images