A method for precise analysis of differential RNA editing sites based on signal enhancement preprocessing

Abstract

The application discloses a method for detecting differential RNA editing sites. The method comprises the following steps: performing high-throughput sequencing and alignment on a sample set comprising 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; performing formal joint mutation detection on the DNA and RNA alignment data with the candidate RNA editing sites as references to obtain RNA editing sites; performing depth counting on the RNA editing sites and forming a depth table; and using a statistical method to evaluate the degree of differential editing of each editing site between different treatment conditions to screen out RNA editing sites with statistical significance. The method can be used for high-precision and high-specificity identification of differential RNA editing sites of various types of RNA editing, and is helpful for explaining RNA editing mechanisms, discovering potential biomarkers and assisting in diagnosis and treatment strategy making of related diseases.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of bioinformatics and biotechnology, and specifically relates to a method composed of a bioinformatics analysis process for identifying differential RNA editing events, which can recognize C>U type or A>I type RNA editing with high precision and high specificity. BACKGROUND

[0002] RNA editing is an important post-transcriptional modification mechanism that regulates gene expression and function by changing 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 usually low in abundance and distributed sparsely, and due to factors such as single nucleotide variations (SNVs) at the DNA level, accurate detection is still challenging despite the use of high-throughput RNA sequencing (RNA-seq) technology.

[0003] Differential variants on RNA (DVR) refers to RNA single nucleotide sites with significant differences in editing degree between different experimental conditions or sample groups. Differential RNA editing sites need to meet the statistical significant difference in editing degree under different conditions, making them have higher reliability and accuracy in detection, and thus are considered as high-confidence RNA editing sites.

[0004] However, when using existing technologies to detect differential RNA editing sites, the following main problems still exist: First, sequencing noise interference and false positive problems, the existing models (such as rMATS-DVR or JACUSA2) have insufficient noise filtering, which can lead to a high false positive rate; second, the inability to effectively carry out DNA and RNA combined mutation detection, existing detection tools (such as VaDiR) attempt to use whole genome sequencing (WGS) data of the same sample to filter DNA mutation interference, but these methods are usually difficult to retain true RNA editing signals in subsequent correction steps; in addition, the detection ability for complex editing types (such as C>U) is limited. Compared with the more common A-to-I editing, the signal of low-abundance editing types such as C-to-U is often weaker in the sequencing background noise, and existing technologies have significant limitations in the accurate identification of such editing.

[0005] In addition, base quality score recalibration (BQSR) is usually performed before mutation detection. This step is to correct the systematic bias 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 currently limited RNA editing sites, conventional BQSR methods may not effectively identify and correct these emerging variants. If BQSR is only performed relying on known reference sites, it may lead to misjudgment of these unknown editing sites as sequencing noise, thereby missing real editing events.

[0006] Therefore, there is an urgent need in the art for a high-sensitivity, high-specificity, and wide-coverage RNA editing site detection method. SUMMARY

[0007] The purpose of the present application is to provide a method for detecting differential RNA editing sites based on signal optimization.

[0008] In a first aspect, the present application provides a method for detecting differential RNA editing sites, comprising the steps of:

[0009] A) providing an independent sample set containing N samples to be detected for differential RNA editing sites; wherein N is a positive integer ≥ 2;

[0010] wherein each of the independent sample set comprises: (i) RNA alignment data of each of the samples, which is obtained by aligning the RNA sequencing data of a single sample with a reference genome, denoted as a first data set; (ii) DNA alignment data of each of the samples, which is obtained by aligning the DNA sequencing data of a single sample with a reference genome, and is subjected to necessary sequence information preprocessing and base quality score recalibration, denoted as a second data set;

[0011] B) performing joint mutation detection on the first data set and the second data set of each sample to be detected for differential RNA editing sites, denoted as a third data set;

[0012] merging the third data set of each sample, thereby obtaining a known RNA editing site data set, denoted as a fourth data set;

[0013] C) inputting the fourth data set as a reference to perform base quality score recalibration on the first data set of all samples; collecting the first data set after base quality score recalibration in all samples, retaining independent sample information, denoted as a fifth data set;

[0014] collecting the second data set in all sample sets, retaining independent sample information, denoted as a sixth data set;

[0015] D) performing joint mutation detection on the fifth data set and the sixth data set, thereby obtaining a candidate RNA editing site data set, denoted as a seventh data set;

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

[0017] F) merging the eighth data set of each sample to obtain a total allele depth table, denoted as a ninth data set; and performing statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth data set, thereby determining the differential RNA editing site.

[0018] In another preferred embodiment, the independent sample set comprises two or more samples of different experimental treatments, samples of normal tissue and tumor tissue, and two or more samples of time series.

[0019] In another preferred embodiment, the independent sample set contains 2 samples.

[0020] In another preferred embodiment, the independent sample set comprises samples of normal tissue and tumor tissue.

[0021] In another preferred embodiment, the samples are repeated.

[0022] In another preferred embodiment, the repetition is selected from the group consisting of biological repetition, technical repetition, or a combination thereof.

[0023] In another preferred embodiment, the repetition is biological repetition.

[0024] In another preferred embodiment, the number of biological repetitions is 4.

[0025] In another preferred embodiment, the step of obtaining the sequencing data comprises:

[0026] (s1) extracting nucleic acids, wherein the nucleic acids comprise DNA and RNA;

[0027] (s2) constructing sequencing libraries, wherein the sequencing libraries comprise DNA sequencing libraries and RNA sequencing libraries;

[0028] (s3) sequencing, thereby obtaining the sequencing data.

[0029] In another preferred embodiment, in step (s1), the nucleic acids are extracted using Roche MagNA Pure 96 instrument and its matching reagents.

[0030] In another preferred embodiment, the A260 / A280 of the DNA is 1.8-2.0.

[0031] In another preferred embodiment, the RIN value of the RNA is >8.

[0032] In another preferred embodiment, the DNA sequencing library is a whole genome DNA sequencing library.

[0033] In another preferred embodiment, the sequencing depth of the DNA sequencing library is >30x.

[0034] In another preferred embodiment, the RNA sequencing library is constructed by an operation selected from the group consisting of:

[0035] (1) constructing a strand-specific library;

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

[0037] or a combination thereof.

[0038] In another preferred embodiment, the strand-specific library is constructed by a first joint method.

[0039] In another preferred embodiment, the proportion of rRNA is reduced by a method selected from the group consisting of a ribosomal RNA removal method, a Poly(A) enrichment method, or a combination thereof.

[0040] In another preferred embodiment, the read length of the RNA sequencing library is dual-end 100bp.

[0041] In another preferred embodiment, the sequencing depth of the RNA sequencing library is 60 million reads per sample.

[0042] In another preferred embodiment, the sequencing data is filtered.

[0043] In another preferred embodiment, the operation of filtering is selected from the group consisting of:

[0044] (I) removing adapter sequences;

[0045] (II) removing low-quality bases;

[0046] or a combination thereof.

[0047] In another preferred embodiment, the software for performing the filtering comprises cutadapt.

[0048] In another preferred embodiment, step (A) specifically comprises:

[0049] (A1) providing RNA sequencing data from N samples, and DNA sequencing data from the N samples, wherein N is a positive integer >2;

[0050] (A2) aligning the RNA sequencing data from the N samples to a reference genome respectively to obtain N RNA aligned data and correct;

[0051] (A3) aligning the DNA sequencing data from the N samples to a reference genome respectively to obtain N DNA aligned data and correct;

[0052] wherein, steps (A2), (A3) can be interchanged, performed in sequence or simultaneously.

[0053] In another preferred embodiment, in step (A2), the software for aligning the RNA sequencing data to a reference genome comprises: STAR, HISAT2, preferably STAR.

[0054] In another preferred embodiment, in step (A2), the software for aligning the RNA sequencing data to a reference genome is STAR.

[0055] In another preferred embodiment, in step (A2), the RNA sequencing data is aligned to a reference genome using the 2-pass method in STAR.

[0056] In another preferred embodiment, in step (A3), the software for aligning the DNA sequencing data to a reference genome comprises: BWA.

[0057] In another preferred embodiment, in step (A3), the DNA sequencing data is aligned to a reference genome using the MEM algorithm in BWA.

[0058] In another preferred embodiment, the DNA sequencing data is whole genome sequencing data.

[0059] In another preferred embodiment, the tool for performing the correction comprises: Picard Tools.

[0060] In another preferred embodiment, in step (A2), the correction comprises steps:

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

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

[0063] (A2c) labeling PCR-generated duplicated sequences;

[0064] (A2d) splitting sequencing fragments (containing N fragments) containing variable splicing sites;

[0065] (A2e) locating target regions that need to be locally sequence realigned;

[0066] (A2f) local re-alignment of indel mutations.

[0067] In another preferred embodiment, in step (A2), further comprising the step of: chain-specific splitting.

[0068] In another preferred embodiment, in step (A3), the correction comprises the step of:

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

[0070] (A3b) adding sample set markers to distinguish different sample origins;

[0071] (A3c) marking PCR-generated duplicated sequences;

[0072] (A3d) locating target regions that need local sequence realignment;

[0073] (A3e) local re-alignment of indel mutations.

[0074] In another preferred embodiment, the reference genome is a human reference genome.

[0075] In another preferred embodiment, the human reference genome is GRCh38.

[0076] In another preferred embodiment, in step (B), the joint mutation calling is performed with the first dataset of each sample as the experimental group and the second dataset of each sample as the control group.

[0077] In another preferred embodiment, in step (B), the software performing the joint mutation calling comprises: GATK.

[0078] In another preferred embodiment, in step (B), the joint mutation calling is performed using Mutect2 in GATK.

[0079] In another preferred embodiment, in step (B), further comprising the step of:

[0080] (B1) quality control, obtaining data passing the quality control;

[0081] (B2) screening effective variations, obtaining the known RNA editing site dataset;

[0082] (B3) data indexing, obtaining the indexed known RNA editing site dataset.

[0083] In another preferred embodiment, in step (B1), the software performing the quality control comprises: GATK.

[0084] In another preferred embodiment, in step (B1), the quality control is performed using FilterMutectCalls in GATK.

[0085] In another preferred embodiment, in step (B1), the maximum number of variant events allowed in each assembled region (max-events-in-region) is 4.

[0086] In another preferred embodiment, in step (B1), the data passing the quality control is tagged.

[0087] In another preferred embodiment, in step (B2), the software performing the filtering includes bcftools.

[0088] In another preferred embodiment, the valid variants are filtered by the tag.

[0089] In another preferred embodiment, the known RNA editing site dataset is RNA mutation site dataset Z.

[0090] In another preferred embodiment, in step (B3), the software performing the data indexing includes GATK.

[0091] In another preferred embodiment, in step (B3), the data indexing is performed using IndexFeatureFile in GATK.

[0092] In another preferred embodiment, in step (C), the software performing the base quality recalibration includes GATK.

[0093] In another preferred embodiment, in step (C), the base quality recalibration is performed using BaseRecalibrator and ApplyBQSR in GATK.

[0094] In another preferred embodiment, in step (D), joint mutation calling is performed using the fifth dataset as the experimental group and the sixth dataset as the control group.

[0095] In another preferred embodiment, in step (D), the software performing the joint mutation calling includes GATK.

[0096] In another preferred embodiment, in step (D), the joint mutation calling is performed using Mutect2 in GATK.

[0097] In another preferred embodiment, in step (D), further comprising a step:

[0098] (D1) mutation site filtering;

[0099] (D2) Single nucleotide variant site filtering;

[0100] (D3) Data processing.

[0101] In another preferred embodiment, in step (D1), the software performing the mutation site filtering comprises: GATK.

[0102] In another preferred embodiment, in step (D1), the mutation site filtering is performed using FilterMutectCalls in GATK.

[0103] In another preferred embodiment, in step (D1), the maximum number of variant events allowed within each assembled region is 4.

[0104] In another preferred embodiment, in step (D2), the software performing the single nucleotide variant site filtering comprises: GATK.

[0105] In another preferred embodiment, in step (D2), the single nucleotide variant site filtering is performed using SelectVariants in GATK.

[0106] In another preferred embodiment, in step (D3), further comprises a step of:

[0107] (D3a) Homopolymer nucleotide sequence filtering;

[0108] (D3b) Genomic repeat sequence filtering.

[0109] In another preferred embodiment, in step (D3a), the software performing the homopolymer nucleotide sequence filtering comprises: SNPiR.

[0110] In another preferred embodiment, in step (D3a), the homopolymer nucleotide sequence filtering is performed using filter_homopolymer_nucleotides.pl in SNPiR.

[0111] In another preferred embodiment, in step (D3b), the software performing the genomic repeat sequence filtering comprises: SNPiR.

[0112] In another preferred embodiment, in step (D3b), the genomic repeat sequence filtering is performed using pblat_candidates_ln.pl in SNPiR.

[0113] In another preferred embodiment, step (E) specifically comprises:

[0114] (E1) Splitting the first data set from N samples by positive and negative strands to obtain N positive strand subsets and N negative strand subsets;

[0115] (E2) calculating the edit site depth data of each candidate RNA editing site in the seventh dataset in two subsets, respectively.

[0116] In another preferred embodiment, the software performing step (E2) comprises: SAMtools.

[0117] In another preferred embodiment, step (E2) is performed using mpileup in SAMtools.

[0118] In another preferred embodiment, the candidate RNA editing site information comprises: the chromosome where the candidate RNA editing site is located, the base position where the candidate RNA editing site is located, the strand where the candidate RNA editing site is located, the reference genome base type (REF base), the edited base type (ALT base).

[0119] In another preferred embodiment, the edit site depth data comprises: reference base depth, edited base depth, base editing ratio.

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

[0121] In another preferred embodiment, the model is selected from the group consisting of: generalized linear mixed model, Dirichlet-multinomial distribution model, or a combination thereof.

[0122] In another preferred embodiment, the model is a combination of generalized linear mixed model and Dirichlet-multinomial distribution model.

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

[0124] In another preferred embodiment, in step (F), the software performing the statistical analysis comprises: rMATS-DVR.

[0125] In another preferred embodiment, the threshold of the FDR correction is ≤0.05.

[0126] In another preferred embodiment, in step (F), it further comprises the step of: annotating the differential RNA editing sites.

[0127] In a second aspect of the present application, a device or system for detecting differential RNA editing sites is provided, which comprises:

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

[0129] (M2) a scanning module configured to perform the following operations: performing joint mutation calling on the RNA aligned data and DNA aligned data of each of the test samples to obtain a preliminary RNA editing site dataset;

[0130] (M3) a detecting module configured to perform the following operations: collecting the corrected RNA aligned data of each of the test samples to obtain a fifth dataset; collecting the DNA aligned data of each of the test samples to obtain a sixth dataset; performing joint mutation calling on the fifth dataset and the sixth dataset to obtain a seventh dataset; obtaining allele depth of each candidate RNA editing site in the seventh dataset from the RNA aligned data of each of the test samples to obtain an eighth dataset; merging the eighth dataset of each of the test samples to obtain a ninth dataset; performing statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset to determine differential RNA editing sites;

[0131] (M4) an output module configured to output information of the differential RNA editing sites.

[0132] In another preferred embodiment, the N is a positive integer ≥ 2.

[0133] In another preferred embodiment, the scanning module (M2) comprises the following sub-modules:

[0134] (M2.1) a preliminary detecting sub-module configured to perform the following operations: performing joint mutation calling on the RNA aligned data and DNA aligned data of the sample of each of the differential RNA editing sites to be detected to obtain a preliminary RNA editing site dataset;

[0135] (M2.2) a quality control sub-module configured to perform the following operations: performing quality control on the preliminary RNA editing site dataset to obtain a RNA editing site dataset passing quality control;

[0136] (M2.3) an effective variant screening sub-module configured to perform the following operations: screening by the label in the RNA editing site dataset passing quality control to obtain a known RNA editing site dataset;

[0137] (M2.4) a data indexing submodule, configured to index the known RNA editing site dataset, thereby obtaining an indexed known RNA editing site dataset;

[0138] (M2.5) a quality correction submodule, configured to perform base quality score correction on the RNA sequencing data, taking the indexed known RNA editing site dataset as input, thereby obtaining corrected RNA sequencing data.

[0139] In another preferred example, the detecting module (M3) comprises the following submodules:

[0140] (M3.1) an official detecting submodule, configured to perform joint mutation detection on the fifth dataset and the sixth dataset, thereby obtaining an official RNA editing site dataset;

[0141] (M3.2) a mutation site filtering submodule, configured to perform mutation site-based quality filtering on the official RNA editing site dataset, thereby obtaining a filtered RNA editing site dataset;

[0142] (M3.3) a screening submodule, configured to extract single nucleotide variant sites in the filtered RNA editing site dataset, thereby obtaining a screened RNA editing site dataset;

[0143] (M3.4) a data processing submodule, configured to perform filtering on the screened RNA editing site dataset, thereby obtaining a candidate RNA editing site dataset, denoted as a seventh dataset;

[0144] (M3.5) a depth sampling submodule, configured to obtain 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 an eighth dataset; and merging the eighth dataset of each of the test samples to obtain a ninth dataset

[0145] (M3.6) a statistical analysis submodule, configured to perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset, thereby determining differential RNA editing sites.

[0146] In another preferred example, in the data processing submodule (M3.4), the filtering comprises:

[0147] (m1) homopolymer nucleotide sequence filtering;

[0148] (m2) genomic repeat sequence filtering.

[0149] In another preferred embodiment, in step (m1), the software performing the homopolymer nucleotide sequence filtering comprises: SNPiR.

[0150] In another preferred embodiment, in step (m1), the homopolymer nucleotide sequence filtering is performed using filter_homopolymer_nucleotides.pl in SNPiR.

[0151] In another preferred embodiment, in step (m2), the software performing the genomic repeat sequence filtering comprises: SNPiR.

[0152] In another preferred embodiment, in step (m2), the genomic repeat sequence filtering is performed using pblat_candidates_ln.pl in SNPiR.

[0153] It should be understood that, within the scope of the present application, all the technical features described above and the technical features described in detail below (such as the examples) can be combined with each other to form new or preferred technical solutions. Due to the limited space, they will not be listed one by one here. BRIEF DESCRIPTION OF DRAWINGS

[0154] Figure 1 The figure shows the specific steps of the analysis flowchart of the present application.

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

[0156] Figure 3 The figure shows the expression levels of different genes induced for expression relative to TPB after 72 hours of doxycycline treatment.

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

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

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

[0160] Figure 7 Sequence logo of nucleotide frequencies around A>G (I) differential RNA editing sites (DVRs) in SK-OV-3 cells mediated by A3B detected using different methods.

[0161] Figure 8 Density plot of RNAFold-based minimum folding energy (MFE) of C>U vs. A>G (I) differential RNA editing sites (DVRs) in T-47D cells mediated by A3B.

[0162] Figure 9 Density plot of RNAFold-based minimum folding energy (MFE) of C>U vs. A>G (I) differential RNA editing sites (DVRs) in SK-OV-3 cells mediated by A3B.

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

[0164] Figure 11 Profile of enhanced cross-linking and immunoprecipitation sequencing (eCLIP-seq) signal of C>U differential RNA editing sites (DVRs) in T-47D cells mediated by A3B detected using different methods, vertical dashed line represents the mean, and the area within the horizontal dashed line represents the 95% confidence interval of the signal.

[0165] Figure 12 Coverage, sequence logo of nucleotide frequencies around, and profile of enhanced cross-linking and immunoprecipitation sequencing (eCLIP-seq) signal of C>U differential RNA editing sites (DVRs) in T-47D cells mediated by A3B detected using the method described in the present application and JACUSA2 RDD-RRD joint detection method, vertical dashed line represents the mean, and the area within the horizontal dashed line represents the 95% confidence interval of the signal. DETAILED DESCRIPTION

[0166] The inventors have proposed a novel RNA editing site detection method with high sensitivity, high specificity and wide coverage of editing types through extensive and in-depth research. Specifically, the inventors first innovatively introduced a "signal optimization pre-scanning" step in the DNA / RNA joint analysis process, which optimized 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, and using these editing sites as references in the subsequent detection process, thereby preserving and enhancing the true RNA editing signal. Experimental results show that compared with existing detection methods, the method significantly improves the detection sensitivity and specificity of editing events. On this basis, the present application is completed.

[0167] The specific implementation steps of the present application are as follows:

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

[0169] The experimental design of the present application includes two or more sample sets of groups, which can be divided according to different experimental treatment methods, normal / tumor tissue controls, and time series, etc. Each sample set should contain at least two types of repeats, which can be biological repeats or technical repeats. In a preferred embodiment, the number of repeats is four, and all are biological repeats.

[0170] In terms of the selection of nucleic acid extraction methods, conventional genomic DNA or total RNA extraction methods can be used, provided that the extracted product meets certain purity requirements. In a preferred embodiment, nucleic acid extraction is completed using Roche MagnaPure96 instruments and their matching reagents, and 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 sequencing instruments can be selected. For the construction of RNA sequencing libraries, conventional construction methods compatible with the used sequencing instruments can also be selected. In the present application, it is particularly recommended to use chain-specific RNA sequencing libraries, which help to determine the directionality of RNA sequencing. The preferred construction method is the first-strand method, i.e., the first-strand synthesis method. In the construction of RNA sequencing libraries, in order to reduce the proportion of rRNA in the library, ribosome removal or poly-A enrichment methods can be selected. In the preferred embodiment of the present application, the ribosome removal method is used. After library construction, single-end sequencing or double-end sequencing can be performed, according to the experimental requirements. In a preferred embodiment, double-end sequencing is used to improve the accuracy of sequence alignment and the depth of downstream analysis.

[0172] 2. Sequencing

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

[0174] 3. RNA / DNA joint mutation detection data alignment and preprocessing

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

[0176] RNA data preprocessing: The main purpose is to correct the RNA sequence alignment results as needed to meet the needs of subsequent mutation detection. This part of the work can be completed using the Picard Tools tool set, including: reordering sequences according to the reference genome; adding sample set markers to distinguish different sample sources; marking PCR-generated repeated sequences to reduce false positive rates; splitting sequencing fragments containing variable splicing sites (N fragments) to improve alignment accuracy; locating target regions that require local sequence realignment; local re-sequence alignment for insertion or deletion (indel) mutations to ensure the accuracy of the mutation position. In RNA data preprocessing, base quality correction is not involved.

[0177] DNA data preprocessing: The main purpose is to correct the DNA sequence alignment results as needed to adapt to the needs of subsequent mutation detection. This part of the work can also be completed using Picard Tools, including: reordering sequences according to the reference genome; adding sample set markers to distinguish the source of different samples; marking PCR-generated repeated sequences to reduce false positives caused by PCR amplification; locating all target regions that need local sequence realignment; local re-sequence alignment for insertion or deletion (indel) mutations to ensure the accuracy of the mutation site.

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

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

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

[0181] The main purpose of signal optimization pre-scan is to accurately identify and mark potential RNA editing sites before base quality correction, to reduce the impact of sequencing noise and non-editing mutations on subsequent analysis. Through this step, we can ensure that these potential editing sites are not incorrectly degraded during the BQSR process, thereby improving the accuracy of RNA editing events.

[0182] In signal optimization pre-scan, the present application preferably uses the Mutect2 function in GATK (version 4 and above) tools for DNA / RNA combined mutation detection. In this process, RNA sequencing data is treated as a "tumor" sample, while genomic DNA sequencing data is used as a "normal" control. Through this combined mutation detection strategy, RNA editing events can be effectively distinguished 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: Perform GATK Mutect2 tool on each sample separately to identify preliminary variant sites in RNA data.

[0185] (2) Quality control: Use the FilterMutectCalls tool of GATK to perform quality control on the detected mutation sites. In this step, the maximum number of allowed variation events in each assembly region is set to 4 to prevent false filtering of true variant sites.

[0186] (3) Screening effective variants: Using bcftools to screen and retain single nucleotide variant sites with "PASS" label passed quality control, generating RNA mutation site dataset Z.

[0187] (4) Data indexing: Using the IndexFeatureFile tool of GATK to index dataset Z, providing support for subsequent BQSR steps.

[0188] After completing the signal optimization pre-scan, the present application performs base quality score recalibration (BQSR) on the RNA sequencing data through the BaseRecalibrator and ApplyBQSR modules in the GATK tool chain. BQSR is an important step to ensure the accuracy of RNA editing detection, exclude sequencing errors, and reduce false mutations.

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

[0190] 5. Formal DNA / RNA joint mutation detection

[0191] After the signal optimization pre-scan and base quality score recalibration (BQSR) steps, the present application further adopts a DNA / RNA joint mutation detection step to further accurately identify and distinguish RNA editing sites and DNA mutations. This step further improves the detection specificity and accuracy of RNA editing events by combining RNA sequencing data with genomic DNA sequencing data, and reduces false positives caused by DNA mutations.

[0192] The specific process is as follows:

[0193] (1) Joint mutation detection: Using the Mutect2 tool of GATK to perform joint variant detection on RNA sequencing data and DNA sequencing data. RNA sequencing data is treated as a "tumor" sample, and DNA sequencing data is used as a "normal" control to distinguish between RNA editing events and DNA mutations.

[0194] (2) Mutation site filtering: For the mutation sites identified by Mutect2, the FilterMutectCalls tool is used for quality filtering. In the present application, the maximum number of variant events in the region (max-events-in-region) is set to 4, thereby allowing the detection of RNA editing site clusters and preventing the false filtering of true RNA editing sites. In this process, low-quality variant data is removed to ensure that only mutation sites with high reliability are retained. This filtering step helps to reduce false positives caused by sequencing errors and low-quality data.

[0195] (3) Single nucleotide variant site screening

[0196] The SelectVariants tool is used to screen the variant data, extract single nucleotide variant sites, and output a new high-quality variant data file. This step ensures that the output data set 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 filtered data is further cleaned, including homopolymer sequence filtering and genomic repeat sequence filtering, to remove interfering information such as homologous polynucleotides that may cause misjudgment. After these refined treatments, a more accurate set of mutation sites will be obtained. Finally, through intersection analysis, combining RNA and DNA mutation site data, the final variant result is generated, confirming the mutations that are truly RNA editing events.

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

[0200] Editing depth detection: Perform sequence counting and editing depth calculation of candidate RNA editing sites, obtain the count of mutant sequences and normal sequences, and calculate the editing depth. In the preferred embodiment of the present application, the mpileup command in the samtools software is used for sequence counting to ensure accurate calculation of editing depth. Further collect sample information and prepare a list containing detailed information of candidate RNA editing sites and editing depth. The list should at least contain: the chromosome where the candidate RNA editing site is located, the base position; reference genome base type (REF base) and edited base type (ALT base); the strand where the RNA editing site is located; the editing situation of each sample in the sample set at the editing site, including reference base depth, edited base depth, and the proportion of edited base to all sequences at the site (base editing proportion).

[0201] Differential statistical analysis: Import the RNA editing site information and editing depth list into the rMATS-DVR software package, use the GLMM model to compare whether the editing proportion of the candidate RNA editing site is statistically different in the comparison of each sample set, and calculate the False Discovery Rate (FDR) of each site to evaluate the significance of the editing event under different conditions.

[0202] Final selection of differential RNA editing sites (DVR): According to the results of statistical analysis, RNA editing sites with sequencing depth greater than 10 and FDR less than 5% are selected as the final confirmed DVR.

[0203] 7. Examples of computer instructions and parameters

[0204] (1) De-adaptor and low-quality base filtering, for example, using cutadapt software, the command paradigm is as follows:

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

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

[0207] (2) RNA sequencing alignment, for example, using STAR software, the command paradigm 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 --outSJfilterReadsUnique --outFilterMultimapNmax 1

[0210] (3) DNA sequencing alignment, for example, BWA and SAMtools software, command paradigm 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 file to carry out sequencing group division, for example, Picard tool, command paradigm as follows:

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

[0214] (5) to input RNA / RNA sequencing BAM file for reordering, for example, Picard tool, command paradigm as follows:

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

[0216] (6) DNA / RNA sequence BAM file mark repeat sequence, for example, Picard tool, command paradigm 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) to RNA sequence BAM file processing with long insertion or deletion read, split CIGAR string, for example, GATK tool, command paradigm 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 joint mutation detection, using GATK software as an example, the command paradigm is as follows:

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

[0225] Filtering of preliminary DNA / RNA joint mutation detection results, the command paradigm is as follows:

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

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

[0228] (10) BQSR processing of RNA sequence BAM files, the command paradigm 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) Again, DNA / RNA combined mutation detection, for example, GATK software, command paradigm 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 again

[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] Extract single base variation

[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 repeat sequence regions, for example, Perl scripts in SNPiR software package, command paradigm as follows:

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

[0240] (13) Prepare BAM files for pblat analysis, generate reference files by merging all RNA-seq BAM files, and sort and index, command paradigm 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 Perl scripts in SNPiR software 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 -w a -header >${label}.final.vcf

[0245] (15) Edit depth sampling to obtain allele depth table, for example, samtools, command paradigm as follows:

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

[0247] (16) Statistical analysis of DVR, using the script of rMATS-DVR software package to analyze the variant depth, count the significant RNA editing events, and calculate the false discovery rate (FDR) using the FDR.py script. 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 5T ${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) Annotation and summary of 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. Previous detection of differential RNA editing site method: data preprocessing

[0254] Select sequencing data from the same source as the examples described above, including: whole genome sequencing data; for providing reference genomic DNA information; RNA sequencing data under corresponding conditions, including induction group and control group, each group has at least three biological repeats. The WGS and RNA-seq raw data (fastq format) are subjected to conventional quality control and de-ligation treatment to ensure that the data quality meets the requirements of subsequent analysis. Then use the same version of human genome (such as GRCh38) and the corresponding gene annotation information (such as GENCODE GRCh38.p13) for sequence alignment.

[0255] 9. The development of the differential RNA editing site detection means described in the prior art: JACUSA method

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

[0257] RDD and RRD mode detection: RDD (RNA-DNA Difference) mode: align the RNA alignment file with the WGS or reference genome to detect the difference between RNA and DNA sequences. RRD (RNA-RNA Difference) mode: only compare the RNA sequencing results of the induction group and the control group to detect the sites that differ under different treatment conditions. In this embodiment, to be consistent with the "double comparison" of the process of the present application, "RDD and RRD mode joint analysis" is also performed, that is, the results of the two modes are combined or intersected respectively to obtain a candidate set of differential editing sites.

[0258] Filtering and result output: apply the JACUSA built-in variant filtering criteria (such as sequencing coverage, quality score, homologous region determination, etc.) to obtain the preliminary candidates. Further, the allele frequency difference of the candidate variants in different groups is calculated, and the sites with significant differences are output, that is, the DVR identified by JACUSA.

[0259] 10. The development of the differential RNA editing site detection means described in the prior art: rMATS-DVR method

[0260] Alignment file preparation: use the same RNA-seq alignment method (such as STAR alignment) as the present application to generate RNA alignment files (BAM). According to the developer's suggestion, no additional DNA correction is introduced; if WGS information is used, it is usually only used for conventional BQSR or filtering.

[0261] RNA variant calling: According to the developer's suggestion, use HaplotypeCaller in GATK package to call variants for RNA data, and obtain the initial candidate RNA variant set. Because rMATS-DVR does not have the DNA-RNA combined mutation detection function of the method described in the present application, part of the SNV of DNA may be misidentified as RNA variants.

[0262] Difference analysis: On the basis of obtaining candidate sites, compare the variant allele depth in the RNA sequencing of the induction group and the control group, and use the original GLMM statistical model of rMATS for difference detection. If the editing frequency of a variant site has a significant difference between the two treatment conditions and passes the internal filtering threshold, it is determined as DVR and output.

[0263] 11. The development of the differential RNA editing site detection method described by the previous person: VaDiR method

[0264] Single-condition RNA variant calling: VaDiR focuses on RNA-DNA variant comparison in development and design, but often only supports single sample comparison or is not specifically designed for multiple replicates of omics data, and cannot directly perform joint statistical modeling on multiple biological replicates. In order to be consistent with the present application as much as possible, the present embodiment first performs RNA variant calling on the induction group and the control group respectively, and excludes obvious low-quality or suspected repetitive regions.

[0265] Different screening and difference determination: Among the sites detected in the induction group and the control group, screen the variant sites with high frequency and meet the internal filtering standards of VaDiR (such as sequencing coverage ≥10, sequencing quality score ≥20, etc.), compare the allele frequencies between the two groups, and if the frequency difference reaches the specified threshold, it is determined as a differential editing site.

[0266] It should be understood that the following description of specific methods and experimental conditions of the present application in various degrees of detail is provided to provide a substantial understanding of the present application. The definitions of certain terms used in the present specification are provided below. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as generally understood by one of ordinary skill in the art to which the present application belongs.

[0267] Terms

[0268] Where a numerical range is provided, unless the context clearly indicates otherwise, it is intended to include every calculated intermediate value, to the tenth of the unit of the lower limit, to the higher limit and any other intermediate value within the stated range, even if the base unit is not explicitly stated. The upper and lower limits of these smaller ranges can be independently combined to form further ranges, and are also contemplated within the application, subject to any explicitly excluded limit in the stated range. For example, "1 to 50" includes "2 to 25," "5 to 20," "25 to 50," "1 to 10," and the like.

[0269] As used herein, the terms "containing" or "including" can be open, semi-closed and closed. In other words, the terms also include "consisting essentially of" or "consisting of."

[0270] As used herein, the term "and / or" relates to and encompasses any and all possible combinations of one or more of the associated listed items.

[0271] Cytidine Deaminase

[0272] As used herein, the term "cytidine deaminase" refers to a class of enzymes that can convert cytosine (or cytidine) to uracil (or uridine) by removing the amino group. Such enzymes not only play a key role in regulating genomic mutation and diversity at the DNA level, but also can mediate C>U nucleotide base changes in the RNA editing process, affecting gene transcription and translation, thereby playing an important function in various physiological and pathological processes (such as tumor occurrence).

[0273] RDD

[0274] As used herein, the term "RDD" (RNA-DNA Difference) refers to the identification of sites with base differences at the RNA level from DNA sequences by aligning genomic DNA (gDNA) and RNA complementary DNA (cDNA) sequences. The analysis at this stage aims to eliminate false positives due to genomic single nucleotide variations (SNVs) and extract variations caused by RNA editing.

[0275] RRD

[0276] As used herein, the term "RRD" (RNA-RNA Difference) refers to the identification of RNA differential editing sites with significant editing level changes between two conditions by comparing RNA sequencing (RNA-seq) data under different experimental conditions or biological states, and statistically analyzing the editing depth difference of the variation sites. The analysis at the RRD stage emphasizes the influence of experimental condition variability on RNA editing events.

[0277] TPR

[0278] As used herein, the term “TPR” (True Positive Rate) refers to the detection ability of a model or method for true RNA editing events, defined as the proportion of the number of detected true positive events to the number of all actual positive events.

[0279]

[0280] wherein: TP (True Positives) is the number of detected true positive events; FN (False Negatives) is the number of missed true positive events. TPR reflects the sensitivity of the detection method, and is an important indicator for measuring the detection ability of editing sites.

[0281] Precision

[0282] As used herein, the term “Precision” refers to the proportion of true positive events in all RNA editing events detected by a model or method.

[0283]

[0284] wherein: TP (True Positives) is the number of detected true positive events; FP (False Positives) is the number of events that are falsely detected as positive. Precision is an important parameter for evaluating the reliability of the detection method, and the higher the Precision indicates the fewer false positive events and the more reliable the results.

[0285] Accuracy

[0286] As used herein, the term “Accuracy” refers to the overall judgment accuracy of a model or method for all editing events (including positive and negative), defined as the proportion of the number of correctly classified events to the number of all classified events.

[0287]

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

[0289] F-score

[0290] As used herein, the term "F-score" is an index that combines Precision and TPR, defined as the harmonic mean of Precision and TPR.

[0291]

[0292] wherein: Precision, TPR are as previously described. 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.

[0293] APOBEC3B (A3B)

[0294] As used herein, the term "APOBEC3B (A3B)" is a member of the cytidine deaminase family, whose functions include deamination modification of cytosine residues in DNA or RNA molecules. APOBEC3B plays an important role in viral defense and tumor formation, and its abnormal expression can lead to genomic instability and accumulation of a large number of mutations, thereby promoting tumor occurrence and development. Therefore, in the RNA editing research involved in the present application, APOBEC3B is used as a typical model enzyme for C>U editing events at the RNA level.

[0295] DNA / RNA joint mutation detection (DNA / RNA Joint Variant Calling)

[0296] In the present application, "DNA / RNA joint mutation detection" refers to simultaneously using DNA sequencing data and RNA sequencing data in the same sample to identify and distinguish base variations at the RNA level (such as RNA editing) from single nucleotide variations (SNVs) or insertion-deletion variations carried by the genome itself, etc. Through this process, false positives caused by genomic variations can be accurately excluded, and the specificity and accuracy of RNA editing site detection can be improved.

[0297] Single nucleotide variant (SNV)

[0298] As used herein, the term "SNV" refers to a change in a single nucleotide in a genomic sequence, with only a single base difference compared to the reference sequence. SNV can occur in coding regions, non-coding regions, or regulatory regions, and can cause changes in protein coding (such as missense mutations) or functional effects.

[0299] RNA editing site

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

[0301] Differentially Edited RNA Site (DVR)

[0302] As used herein, the term "differentially edited RNA site" refers to an editing site whose RNA editing level (base substitution ratio) shows significant difference between two or more groups, different experimental conditions or different time points. The present application quantifies and compares the editing depth of each sample by generalized linear model or other robust statistical methods, screens the RNA editing sites significantly affected by experimental or environmental conditions, and establishes the "differentially edited RNA site (DVR)" described in the present application.

[0303] Signal optimization pre-scanning

[0304] In the process of DNA / RNA combined analysis, the present application introduces signal optimization pre-scanning to improve the signal-to-noise ratio of the process, and improve the accuracy and specificity of detection.

[0305] In signal optimization pre-scanning, through multi-level combined analysis and pre-labeling strategy, the signal retention and false positive filtering effect in RNA editing site detection are significantly improved.

[0306] Firstly, it performs preliminary labeling of candidate RNA editing sites in the preliminary screening stage of RNA sequencing data through combined DNA / RNA mutation detection. This preprocessing strategy can maximize the identification of potential true RNA editing signals and avoid being down-weighted or lost in the subsequent BQSR step.

[0307] Secondly, by combining DNA and RNA mutation analysis, signal optimization pre-scanning selects RNA editing events with high confidence rather than DNA editing events, reducing false positive sites caused by background noise or artifacts. Especially in low complexity sequences and genomic repeat regions, precise filtering can reduce the interference of false positives. For complex editing types such as C-to-U, signal optimization pre-scanning combined with DNA and RNA data can effectively improve the detection efficiency of low-abundance editing signals. In the preliminary screening stage, these complex editing sites are highlighted, providing a strong foundation for subsequent differential analysis.

[0308] Base Quality Score Recalibration (BQSR)

[0309] As used herein, the term "base quality correction" and "BASR" can be used interchangeably, which is a common step before detecting variants. Base quality correction is to correct systematic bias that may be introduced during sequencing, thereby improving the accuracy of variant detection. When the sequencer reads the DNA sequence, the quality score of some bases may be low due to factors such as chemical reactions and instrument performance. These low-quality bases may be misjudged as mutations, affecting the reliability of subsequent analysis. By correcting these low-quality bases in the BQSR step, the interference of sequencing errors on mutation detection results can be effectively reduced, and the accuracy and reliability of mutation detection can be improved. Currently, BQSR is mainly based on known reference sequences to adjust (usually lower) the quality score of mutant bases, reducing the interference of sequencing noise on mutation detection.

[0310] The main advantages of the present application include:

[0311] (1) The present application develops a method for detecting RNA editing sites, and for the first time introduces a signal optimization pre-scanning step in the method, which reduces the number of bases misjudged as noise by re-correcting the base quality score, thereby significantly improving the retention rate of editing signals, and further significantly improving the number of RNA editing events detected, enhancing the identification of true RNA editing signals, greatly reducing the interference of background noise, and having high sensitivity.

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

[0313] (3) The present application uses a generalized linear mixed statistical analysis method to analyze the difference of RNA editing sites, which can significantly improve the detection rate and reliability of differential RNA editing sites compared with traditional methods, and has high accuracy.

[0314] (4) The present application provides a method for detecting RNA editing sites and / or RNA editing mechanisms in different biological conditions, different experimental treatments or disease states at the whole transcriptome and whole genome level.

[0315] The application will be further described in conjunction with specific examples. It should be understood that these examples are only used to illustrate but not limit the scope of the application. The experimental methods in the following examples, if not otherwise specified, are generally carried out according to the conventional conditions, for example, the conditions described in Sambrook et al., Molecular Cloning: A Laboratory Manual (New York: Cold Spring Harbor Laboratory Press, 1989), or the conditions recommended by the manufacturer. Unless otherwise specified, percentages and parts are weight percentages and weight parts.

[0316] Materials and methods

[0317] 1. Data acquisition and sequence alignment

[0318] 1.1 Sample preparation

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

[0320] Library construction:

[0321] RNA sequencing library: A strand-specific library construction method (such as First-Strand method) is preferably used 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.

[0322] DNA sequencing library: A standard whole-genome DNA library construction method is used to ensure that the sequencing depth reaches more than 30x to provide high-quality genomic data.

[0323] Sequencing platform: A high-throughput sequencing platform (such as BGISeq-500, MGISeq, or Illumina HiSeq) is selected, and the RNA sequencing read length is preferably double-end 100bp, and the sequencing depth is 60 million reads per sample.

[0324] 1.2 Sequence alignment and preprocessing steps

[0325] RNA sequencing data alignment:

[0326] The STAR software (version 2-pass mode) is used to align the RNA-seq data to the reference genome (such as GRCh38) and the reference transcriptome (such as Gencode).

[0327] Perform initial de-duplication, low quality sequence filtering, and repeat sequence marking to generate corrected RNA alignment files (.bam format).

[0328] DNA sequencing data alignment:

[0329] Align WGS data to the same reference genome using BWA-MEM algorithm.

[0330] Perform de-duplication, low quality filtering, and repeat sequence marking to generate corrected DNA alignment files (.bam format).

[0331] 2. Signal optimization pre-scanning step

[0332] Tool selection: The present application performs DNA / RNA joint mutation detection by preferably using the Mutect2 function in the GATK (version 4 and above) tool, in which the RNA sequencing data is treated as a "tumor" sample, while the DNA sequencing data is treated as a "normal" control. This joint mutation detection strategy can effectively distinguish RNA editing sites and genomic DNA mutations, optimizing the detection results of RNA editing events.

[0333] Mutation detection process:

[0334] Perform gatk Mutect2 mutation detection for each sample separately, from which the preliminary RNA mutation sites are generated.

[0335] Use the gatk FilterMutectCalls tool to perform quality control on mutation sites, setting the maximum number of allowed variation events within a single assembly region to 4 to prevent real variations from being mis-filtered.

[0336] Use bcftools to retain single nucleotide variations with the "PASS" label that pass quality control to form the RNA mutation site dataset Z.

[0337] Use gatk IndexFeatureFile to index dataset Z in preparation for the subsequent BQSR step.

[0338] Purpose: The main purpose of the signal optimization pre-scanning is to accurately identify potential RNA editing sites before performing BQSR correction, and to label these sites as known variations to ensure that they are not mis-lower quality scores during the subsequent base quality score correction process. This step significantly reduces false positives caused by sequencing noise or non-editing mutations, improving the reliability and accuracy of RNA editing sites.

[0339] 3. Formal DNA / RNA joint mutation detection

[0340] 3.1 BQSR correction

[0341] Tool selection: In the process of RNA editing detection, BQSR correction is a key step to ensure high accuracy of mutation sites and reduce sequencing noise. The present invention uses BaseRecalibrator and ApplyBQSR modules in GATK tool chain to correct the quality of read data when analyzing RNA-seq data, thereby excluding false mutations caused by errors in the sequencing process.

[0342] Correction process: Use the list of known variant sites generated by the "signal optimization pre-scan" step as the "known sites" part of the reference input for BQSR. Perform BQSR correction.

[0343] 3.2 DNA / RNA joint mutation detection

[0344] Tool selection: The present invention uses GATK Mutect2 tool for DNA / RNA joint mutation detection, combining RNA sequencing and genomic DNA sequencing data, which can effectively distinguish RNA editing from DNA mutations.

[0345] Process:

[0346] Joint mutation detection: Use Mutect2 tool, use RNA sequencing data as "tumor" sample, DNA sequencing data as "normal" control, to identify RNA editing events and DNA mutations.

[0347] Filter mutation sites: Use FilterMutectCalls tool to filter the variant data generated by Mutect2, only keep the mutation sites that meet the standard, and exclude low-quality data.

[0348] Select single nucleotide variant sites: Use SelectVariants tool to screen out single nucleotide variant sites, output new high-quality variant data file.

[0349] Subsequent processing: Further filter the screened single nucleotide variant site data, remove homologous polynucleotide interference information, and ensure the accuracy of mutation sites. Specifically, homopolymer sequence filtering: remove mutations located in homopolymer repeat sequences to prevent false positives caused by sequencing errors; genomic repeat sequence filtering: remove mutations located in highly repetitive regions of the genome to reduce multiple alignment errors. Finally, perform intersection analysis of RNA and DNA mutation sites to generate the final variant result.

[0350] Objective: By DNA / RNA combined mutation detection, the application effectively distinguishes RNA editing events and DNA mutations, significantly improves the detection specificity of RNA editing sites, avoids misjudgment caused by DNA mutations, and ensures more accurate RNA editing detection.

[0351] 4. Editing depth sampling and difference statistical analysis

[0352] 4.1 Editing depth sampling

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

[0354] Implementation steps:

[0355] Split the RNA-seq alignment file by positive and negative strands, and generate n positive strand subsets and n negative strand subsets respectively.

[0356] For each candidate RNA editing site, calculate the depth of the reference base and the edited base in each sample.

[0357] Calculate the editing proportion (editing base depth / total depth) of each site in each sample.

[0358] 4.2 Difference statistical analysis

[0359] Statistical method: Use generalized linear mixed model (Generalized Linear Mixed Model, GLMM) or Dirichlet-multinomial distribution model to perform difference analysis on RNA editing proportions between multiple biological conditions or experimental groups.

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

[0361] Implementation steps:

[0362] Organize the editing depth data of each candidate RNA editing site in each sample into a matrix file.

[0363] Apply statistical models to matrix data to assess whether there is a significant difference in editing proportion of each site under different conditions.

[0364] Perform multiple hypothesis testing (such as Benjamini-Hochberg FDR correction) to screen out statistically significant DVRs (such as FDR≤0.05).

[0365] 4.3. Confirmation and output of final DVR

[0366] Criteria for validation: RNA editing sites with sequencing depth ≥10 and FDR ≤0.05 in differential statistical analysis were selected as the final differential RNA editing sites (DVR).

[0367] Output content: A DVR list containing the following information was generated: chromosomal location (chromosome number, base position); reference base type (REF base); edited base type (ALT base); chain (positive or negative); editing depth and editing proportion in each sample; statistical analysis results (such as T value, FDR, etc.).

[0368] Example 1: Construction of virtual data required for evaluation of differential RNA editing site detection efficiency

[0369] The present description evaluates the effectiveness of the method of the present application in the detection of differential RNA editing sites (DVR) 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 whole genome sequencing (WGS) and transcriptome sequencing (RNA-seq) virtual data.

[0370] In terms of data construction, this embodiment uses WGS and RNA-seq data types for simulation. The WGS data uses a 2x150nt sequencing strategy, with an average sequencing depth of 33x, and 2 biological replicates are set; the RNA-seq data uses a 2x100nt strand-specific sequencing strategy, with a total read length of about 16 million, and 3 to 5 biological replicates are set. In these two types of data, 50,000 SNVs located in DNA and RNA are inserted to simulate the presence of single nucleotide polymorphisms (SNPs) in actual genomes. In addition, 50,000 RNA editing sites that only exist in the RNA layer are additionally inserted in the RNA-seq data to further simulate the environment of real RNA editing events, thereby providing basic data support for subsequent differential RNA editing site detection.

[0371] In order to evaluate the detection effect of differential 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 the "control group" of the two experimental conditions, about 6,000 sites meet the DVR judgment criteria, which are based on the thresholds defined by rMATS GLMM or JACUSA methods. Sites that meet the above conditions are defined as "successfully detected" true positives, which are used to evaluate the sensitivity and accuracy of the method described in the present application in DVR detection.

[0372] To further close to the real application scenario, about 2,000 SNVs meeting the requirements of RNA-RNA detection method are inserted in WGS and RNA-seq data to simulate false positive results. This operation aims to investigate the differences in theoretical precision and accuracy between pure RNA-RNA detection mode (without reference to DNA sequencing) and DNA-RNA combined detection mode. Through this design, the upper limit of the theoretical precision and accuracy of different detection methods is further limited, for example, the upper limit of the precision and accuracy of pure RNA-RNA detection mode is limited to about 0.75 and 0.12, respectively, while the upper limit of the precision and accuracy of DNA-RNA combined detection mode is about 0.98 and 0.56, respectively (see Table 2 for details).

[0373] Based on virtual data, the performance of the method described in the application is comprehensively evaluated through a series of indicators. In simulated data, if the detection method can detect RNA editing sites consistent with the set differential RNA editing sites, it is considered positive; false positives are considered false positives beyond the DVR definition range. Based on the number of true positives and false positives, the precision, accuracy and sensitivity of each detection method are calculated to quantify its detection efficiency and reliability.

[0374] The virtual data set and evaluation system established through the above steps can comprehensively compare the differences between the method described in the application and other mainstream DVR detection tools in terms of detection efficiency, precision, etc. Table 1 and Table 2 detail the sequencing parameters used and the statistical information of the inserted RNA editing sites, which can provide data support for subsequent result analysis and comparison.

[0375] Table 1. Overview of simulated data used for evaluation

[0376]

[0377] Table 2. Involving internal reference in simulated data used for evaluation

[0378]

[0379] Example 2: Comparison of methodologies involved in transverse evaluation

[0380] To further verify the effect of the signal optimization pre-scanning step (hereinafter referred to as "the correction step") introduced in the method of the present application on DVR detection, the detection effect of retaining or omitting the correction step in the CADRES process is compared in this embodiment. At the same time, based on the same set of simulation data, the existing JACUSA2 is also used to analyze in RNA-RNA difference (RRD) mode and RNA-DNA difference (RDD) mode, respectively, and then the intersection of the variations obtained by RRD and RDD (i.e. the RDD-RRD joint detection method of JACUSA2) is taken to maximize the potential of DVR detection. In addition, the existing VaDiR and rMATS-DVR methods are also evaluated to investigate the detection efficiency and applicability of each method from multiple angles. The following table (Table 3) provides the comparison results of each method in this embodiment and an overview of their characteristics.

[0381] Table 3 lists the functional characteristics involved in each methodology compared and evaluated in this embodiment, including: whether to perform base quality score recalibration (BQSR), the way to compare RNA-RNA difference (RRD) and RNA-DNA difference (RDD), the degree of support for multi-condition or multi-sample, the number of biological replicates required, the variant detection engine used, the variant depth statistical analysis method, and additional alignment artifact filtering measures, etc. Through this table, one can intuitively understand the differences in the settings of each analysis process at different detection links, which facilitates subsequent comparative analysis of the detection results.

[0382] Table 3. Comparison of methodologies involved in the simulation data horizontal evaluation of the embodiment

[0383]

[0384] Example 3: Evaluation of the effectiveness of differential site detection methods using virtual data

[0385] In this embodiment, the sensitivity and specificity of different detection methods in DVR detection are evaluated through virtual data simulation experiments, and the effect of the number of repeated experiments on detection performance is also considered.

[0386] The experimental results are as follows Figure 2The same virtual dataset was used in the experimental design, which was applied to the method described in the present application (with or without the signal optimization pre-scan step), the JACUSA2 method (RDD-RRD joint detection mode), and the traditional method only considering the difference between RNA and DNA (such as VaDiR). These methods were compared under the same conditions, and their performance was quantitatively evaluated using indicators such as true positive rate (TPR), precision, accuracy, and F1 score. The true positive rate (TPR) is used to measure the ability of the method to correctly identify DVR events, the precision evaluates the proportion of sites in the detection results that are true DVR events, the recall reflects the proportion of true DVR sites that the method can detect, and the F1 score is a performance indicator that balances precision and recall.

[0387] The experimental results Figure 2 show that all detection methods have similar performance in terms of true positive rate (TPR), but there are significant differences in precision and accuracy. In particular, the method described in the present application and the RDD-RRD joint detection method of JACUSA2 effectively reduce the number of false positive sites through the joint analysis of RDD and RRD, and significantly outperform the traditional method (such as VaDiR) that only considers RDD in terms of precision and detection reliability. Specifically, the precision of the method described in the present application is 10-15% higher than that of the RDD-RRD joint detection method of JACUSA2, indicating that the method described in the present application has a significant advantage in reducing false positives. However, the recall of the method described in the present application is slightly lower than that of the RDD-RRD joint detection method of JACUSA2 by 3-5%, indicating that the method described in the present application may have certain limitations in the detection of some marginal editing events.

[0388] After introducing the signal optimization pre-scan step, the performance of the method described in the present application in terms of true positive rate (TPR) and precision is significantly improved. The signal optimization pre-scan can effectively identify and retain true RNA editing signals, while reducing false positives and improving the overall sensitivity of the detection. This result fully demonstrates the key role of the signal optimization pre-scan in the method described in the present application, making it have a significant advantage in distinguishing between true editing events and background noise. The experiment further analyzed the effect of the number of experimental repeats on the detection performance, and the results showed that when the number of experimental repeats exceeds 4, the improvement effect of the evaluation indicators tends to be flat. This indicates that increasing the number of experimental repeats within a certain range can improve the stability and reliability of the detection, but the increase in the number of repeats has limited optimization effect on the results beyond a certain threshold.

[0389] This example verifies the high efficiency of the method of the present application in detecting differential RNA editing sites. In particular, after combining the optimized pre-scanning step, the method of the present application is superior to other mainstream methods in terms of accuracy and true positive rate (TPR). After optimizing the number of experimental repetitions, the method of the present application demonstrates its high precision and high sensitivity in RNA editing event recognition. By effectively identifying RNA editing sites with statistical significance, the method of the present application not only improves the reliability of DVR detection, but also proves its potential as an important tool in RNA editing research. This experimental result provides strong support for the scientificity and practical application value of the method of the present application.

[0390] Example 4: Construction of cell lines induced to express cytidine deaminase

[0391] This example describes the construction of a lentivirus-mediated cell line induced to express cytidine deaminase APOBEC3B, which is used to study the function of cytidine deaminase-mediated RNA editing and to verify its feasibility as an evaluation tool for differential RNA editing site detection methods. The T-47D cell line and the SK-OV-3 cell line were selected as model systems, where the T-47D cell line has moderate levels of endogenous A3B expression, while the SK-OV-3 cell line has almost no endogenous A3B expression. The purpose of selecting these two cell lines is to compare A3B-mediated DVR events and provide a comparison between endogenous and exogenous A3B expression.

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

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

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

[0395] This example describes the steps of collecting samples and extracting nucleic acids from the constructed T-47D and SK-OV-3 cells to ensure that the data obtained are suitable for subsequent high-throughput sequencing. The cells were cultured under conventional conditions, and the subculture medium was RPMI-1640 medium, 10% fetal bovine serum and 1% pen / strep for T-47D and SK-OV-3 cells, respectively. In the induced group of cells, doxycycline was added at a concentration of 1 pg / mL, and the cells were induced for 72 hours to ensure sufficient expression of APOBEC3B. The uninduced cells served as a control group and no doxycycline was added.

[0396] After the induction was completed, the cells were separated from the culture plates using trypsin digestion and washed with PBS buffer to remove residual substances in 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 according to the manufacturer's instructions to ensure the integrity and purity of the RNA sample. Genomic DNA extraction was performed using the QIAamp DNA Mini Kit (Qiagen) according to the manufacturer's standard operating procedures to ensure that high-quality DNA samples were obtained. After extraction, the RNA and DNA samples were subjected to quality control. The concentration and purity of the RNA and DNA were determined using a NanoDrop instrument to ensure that the 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 extracted RNA had a RIN value greater than 9.0. The RNA samples were subjected to mRNA enrichment and library construction by ribosome removal. Subsequently, the samples were subjected to high-throughput RNA sequencing using the BGI-SEQ500 platform (Shenzhen Huada Gene).

[0397] The RNA-seq results were counted using the Subread software package, and the results were normalized using the TPM method and further corrected using the TBP gene expression content. The results are shown in Table 1. Table 1: RNA-seq resultsFigure 3 As shown, the level of exogenously induced APOBEC3B expression was significantly increased in T-47D and SK-OV-3 cells after 72 hours of doxycycline treatment, indicating that the induced expression cells were successfully constructed.

[0398] Example 6: Detection of RNA differential editing sites using the method described in the present application

[0399] This example describes the detection of RNA differential editing sites in T47D and SK-OV-3 cells using the method described in the “Detailed implementation method” in the present application. Statistical analysis results show (Tables 4 and 5) that the C>U editing frequency changed significantly after 72 hours of doxycycline induction of cells to express APOBEC3B, and the RNA editing sites were successfully detected and screened.

[0400] Table 4. Comparison of the number of differential RNA editing sites (DVR) detected in APOBEC3B-induced expression T-47D cells using different methods

[0401]

[0402] Table 5. Comparison of the number of differential RNA editing sites (DVR) detected in APOBEC3B-induced expression SK-OV-3 cells using different methods

[0403]

[0404] Example 7: Detection of RNA differential editing sites using control methods (JACUSA, rMATS-DVR, VaDiR)

[0405] This example aims to illustrate the specific process of detecting DVRs under the same or similar experimental conditions as the method described in the present application, starting from the sequencing data of T-47D and SK-OV-3 induced expression cell lines, to verify the improvement and advantage of the present application in detection effect and precision. The results of RNA differential editing sites obtained from the sequencing data of the induced expression cell treatment group / control group described in Examples 4-6 using JACUSA, rMATS-DVR, and VaDiR methods were collected respectively. The number of different types of editing events such as C>U, A>G(I) detected by each method is shown in Tables 4 and 5.

[0406] The data in Table 4 shows that the method of the present application is significantly better than the version without the signal optimization pre-scanning step in the detection of RNA editing events. Specifically, the signal optimization pre-scanning step increases the number of C>U and A>G(I) type of RNA editing events, with C>U DVR sites increasing from 739 to 816 and A>G(I) DVR sites increasing from 1451 to 3764. Such an increase indicates that the signal optimization pre-scanning step can more effectively identify low-abundance RNA editing events, especially complex editing types.

[0407] The experimental results show that there is a significant overlap between the SNVs in T-47D cells and the DVRs detected by rMATS-DVR (see Table 4). This result is as expected because the rMATS-DVR method does not employ a mechanism to filter SNVs. Among the rMATS-DVR analysis containing SNVs, 33% of C>U DVRs showed a decrease in alternative allele frequency after APOBEC3B induction, while this pattern was not observed in the DVRs obtained using the method of the present application. Further testing was performed using the SK-OV-3 cell model, in which the exogenous expression level of APOBEC3B is lower but biologically relevant. The results show that the method of the present application is slightly less sensitive than the RDD-RRD combined detection method of JACUSA2 in the detection of DVRs (see Table 4 and Table 5). However, the overall difference is not significant. This result indicates that the method of the present application can significantly improve the false positive problem caused by SNVs compared to methods that only consider RNA mutations.

[0408] Example 8: Comparison of sequence characteristics of C>U RNA differential editing sites detected by different methods

[0409] This example verifies the ability of different methods in the accuracy of differential editing site detection by comparing the sequence characteristics of C>U RNA differential editing sites detected by the method of the present application and the control method, including upstream and downstream base preference and enrichment of enzyme-specific motif. The sequence characteristics of RNA editing sites are particularly important, and the RNA editing mechanism of APOBEC3B has a high dependence on the sequence context of the surrounding bases. APOBEC3B has a known RNA editing target site preference for specific sequence characteristics (motif), i.e., 5'-UUCM (where M represents A or C). By detecting whether C>U events are enriched in 5'-UUCM motif and specific sequence environment upstream and downstream thereof, it can be more accurately verified whether the editing site is mediated by APOBEC3B, thereby excluding false positive signals.

[0410] In the experiment, first, the sequence of a certain base range (such as ±5 or ±10 nt) upstream and downstream of the DVR site detected by different methods is extracted, and unified filtering and quality control is carried out to ensure the consistency and reliability of all sequence data. For these sequence data, nucleotide spectrum analysis and Logo visualization are carried out using WebLogo and other tools to clearly show the base enrichment pattern of the DNR detected by each method around the editing site. The analysis results Figure 4 , Figure 5 ) show that the C>U sites detected by the method of the present application significantly reflect the enrichment characteristics of the APOBEC3B specific target site motif (5'-UUCM) in the editing core region (such as near the C site), while the control methods (including JACUSA2 RDD-RRD combined detection method, JACUSA2 RDD, JACUSA2 RRD and VaDiR method) fail to fully filter DNA SNV or sequencing noise, and the sequence preference of the C>U sites detected is weak and deviates to a certain extent. In addition, the sequence characteristics of the APOBEC3B differential RNA editing site detected by the method of the present application are significantly different from the sequence characteristics of the DNA editing mediated by APOBEC3B (obtained from dataset GSE193225) Figure 4 ). Through the above sequence characteristic analysis, it is proved that the method of the present application exhibits higher accuracy and effectiveness in the specific recognition of RNA editing events than the existing methods (including JACUSA2 RDD-RRD combined detection, JACUSA2 RDD, JACUSA2 RRD, rMATS-DVR and VaDiR method),

[0411] Example 9: Sequence feature comparison of A>I RNA differential editing sites detected by different methods

[0412] The present embodiment compares the performance of the method of the present application and other detection methods (JACUSA2, rMATS-DVR, VaDiR) in terms of specificity and accuracy by detecting and analyzing A>I RNA editing sites, and further verifies the advantage of the method of the present application in enriching A>I RNA editing site-specific motif (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, and its typical feature is to deaminize adenosine (A) to inosine (I). This editing event usually appears as an A>G substitution in high-throughput RNA sequencing. ADAR enzymes have high dependence on target sequences, and their editing preferences are significantly affected by the adjacent base environment, especially the 5'-YAS motif (where Y represents U or C, and S represents G or C) in the sequence before and after the editing site, which is considered an important marker of ADAR activity. Detecting whether this specific motif can be enriched is a key to verifying the biological authenticity and specificity of RNA editing events.

[0413] The experiment uses two independent cell model data sets, including T-47D induced expression cell lines and SK-OV-3 cell lines. The experimental results show that the method of the present application has higher specificity and accuracy than other methods. Figure 6 、 Figure 7)indicates that the method of the present application shows significant detection ability for A>I editing sites in both T-47D and SK-OV-3 cells. The A>I sites detected by the method of the present application are enriched in 5'-YAS motif in the upstream and downstream sequences, which is highly consistent with the 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, they fail to effectively exclude false positive signals (such as SNV or sequencing noise), resulting in a significant weakness in the enrichment of 5'-YAS motif in the A>I sites detected by the method of the present application. In addition, the A>I sites detected by these control methods contain some signals interfered by SNV, further weakening the biological relevance of the editing sites. In further comparison, although the sensitivity of the method of the present application in the SK-OV-3 cell model is slightly lower than that of the JACUSA combined RDD-RRD combined detection, its significant advantage in the motif enrichment ability of A>I sites indicates that it is more suitable for specific detection. In particular, through the "signal optimization pre-scanning" step, the method of the present application can separate the true RNA editing signal from the noise, significantly reduce the interference of false positive signals, and ensure the high consistency of the detection results of A>I sites with the characteristics of ADAR activity. Through the above analysis, this embodiment verifies the significant advantages of the method of the present application in specificity and accuracy in detecting A>I RNA editing events, especially the accuracy and effectiveness in enriching ADAR preferred motif (5'-YAS).

[0414] Example 10: Comparison of the physicochemical characteristics of the detected RNA differential editing sites

[0415] This embodiment verifies the biological relevance and accuracy of the detected sites by comparing the physicochemical characteristics of the RNA secondary structure energy (Minimum Folding Energy, MFE) of the RNA differential editing sites detected by the method of the present application and the control methods (JACUSA2, rMATS-DVR, VaDiR). The same high-throughput sequencing data of the induced APOBEC3B cells (T-47D and SK-OV-3) in the previous embodiments were used, and the method of the present application and the control methods were run under consistent analysis conditions to obtain the DVR sites output by them. The MFE of the sequence within a certain range (±100 nt) upstream and downstream of each DVR site was calculated using RNA secondary structure prediction tools such as RNAFold, and the MFE distribution characteristics of the DVR obtained by each method were compared.

[0416] Analysis results Figure 8 、 Figure 9)indicate that the DVRs detected by the method of the present application and the control method both exhibit significant RNA secondary structure folding characteristics, and the MFE distribution is consistent with the characteristics of known RNA editing sites. This indicates that the detected DVRs are more likely to be in the region where RNA secondary structure can be formed, which is consistent with the biological characteristics of RNA editing. Overall, this embodiment verifies that the method of the present application and the control method can both identify DVR sites with biological characteristics when detecting RNA differential editing sites, which further confirms that the detected DVRs are real events at the RNA level, rather than DNA SNV or background noise interference.

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

[0418] This embodiment uses enhanced cross-linking and immunoprecipitation sequencing (eCLIP-seq) technology to verify the accuracy and specificity of the method of the present application and the control method in detecting differential RNA editing sites. eCLIP-seq is an experimental technique for studying the interaction of RNA-binding proteins (RBPs) with RNA. By ultraviolet cross-linking, RNA is fixed with its binding protein to form a complex, and then immunoprecipitation is performed to enrich specific protein-nucleic acid complexes. RNA is then isolated, reverse transcribed and analyzed by high-throughput sequencing to determine the binding sites of RNA-binding proteins and the distribution of target RNA. In this study, this technology was used to specifically capture APOBEC3B and RNA complexes by immunoprecipitation, and further sequencing of the bound RNA fragments to obtain the binding sites of APOBEC3B on RNA.

[0419] In this study, the APOBEC3B protein binding signal around the DVRs identified by different detection methods (including the method of the present application, JACUSA2 RDD-RRD joint detection, rMATS-DVR and VaDiR method) was analyzed using published eCLIP-seq data (GSE193225). The results show that the C>U DVR sites detected by the method of the present application are significantly superior to other methods in terms of eCLIP-seq signal enrichment. Specifically, the DVRs detected by the method of the present application exhibit strong APOBEC3B binding signal enrichment in the upstream region of the editing site, which is consistent with the biological mechanism of APOBEC3B playing a role by binding to RNA in the C>U RNA editing process Figure 10 )。

[0420] In contrast, JACUSA2 RDD-RRD combined detection and rMATS-DVR and other methods also detect a certain number of C>U DVR, but the signal enrichment degree is low, and the consistency of the binding mode with APOBEC3B is poor. Compared with the method described in the application, there is a statistical difference in the improvement of signal enrichment degree Figure 11 This result shows that the method described in the application can not only identify RNA editing events, but also accurately capture the signals of RNA binding proteins related to editing events. Especially in capturing the binding signals of APOBEC3B and RNA complex, it has a significant advantage compared with other methods.

[0421] Example 13: Difference between the RNA differential editing sites detected by the method described in the application and the optimal method of the prior art

[0422] This example verifies the improvement and technical advantage of the method in detecting DVR by directly comparing with the currently recognized better RNA differential editing detection method (JACUSA2 RDD-RRD combined detection). The experimental data is derived from high-throughput sequencing of induced APOBEC3B cells (T-47D and SK-OV-3), and the same quality control and alignment strategy is used for analysis to ensure the fairness of the comparison. Through the intersection, union and difference set analysis of the DVR sets detected by CADRES and JACUSA2 RDD-RRD combined detection, the specificity and biological significance of the overlapping sites and unique sites of the two methods in detecting C>U and A>G(I) and other editing types are further evaluated.

[0423] The results show that Figure 12), compared with JACUSA2 RDD-RRD combined detection, the method of the present application shows higher specificity and lower false positive rate, especially in C>U type editing events, the method of the present application can more effectively exclude interference sites derived from sequencing noise. When analyzing overlapping sites, it was found that the common DVRs detected by the method of the present application and JACUSA2 were mostly sites with significantly changed editing rates and high coverage, which were highly consistent with the sequence characteristics of known APOBEC3B RNA editing sites, and eCLIP-seq signals were also significantly enriched. For the RNA differential editing sites unique to the method of the present application, these sites showed stronger enrichment on A3B specific sequence preference (such as 5'-UUCM), and the protein signal near the binding site was more significant, further verifying that these sites are real editing events mediated by A3B. In contrast, the RNA differential editing sites unique to the JACUSA2 RDD-RRD combined detection no longer have the 5'-UUCM sequence feature, but have the 5'-UCA sequence feature. This sequence feature is consistent with the sequence feature of the DNA editing site of APOBEC3B (5'-TCA), and part of the sites are related to genomic variation background or low coverage region, and eCLIP-seq signal does not appear significant enrichment phenomenon, suggesting that the result may contain more false positive signals caused by SNV.

[0424] Through comparison with the current optimal RNA editing detection scheme, the present embodiment shows that the method of the present application can show higher specificity and accuracy in detecting real DVRs, especially in the recognition of C>U type RNA editing events. This method not only provides more reliable basic data for the study of RNA editing mechanism, but also shows higher practical value in subsequent biological function analysis and clinical marker screening.

[0425] All the documents mentioned in the present application are cited as references in the present application, as if each document is cited as a reference individually. In addition, it should be understood that those skilled in the art can make various modifications or changes to the present application after reading the above teachings of the present application, and these equivalent forms also fall within the scope defined by the claims attached to the present application.

Claims

1. A method of detecting differential RNA editing sites, characterized by, The method comprises the steps of: A) providing independent sample sets containing N samples to be detected for differential RNA editing sites; wherein N is a positive integer greater than or equal to 2; wherein each of the independent sample sets comprises: (i) RNA alignment data of each of the samples, which is obtained by aligning RNA sequencing data of a single sample with a reference genome, denoted as a first data set; (ii) DNA alignment data of each of the samples, which is obtained by aligning DNA sequencing data of a single sample with a reference genome, and is subjected to necessary sequence information preprocessing and base quality score correction, denoted as a second data set; B) performing joint mutation detection on the first data set and the second data set of each sample to be detected for differential RNA editing sites, denoted as a third data set; combining the third data set of each sample to obtain a known RNA editing site data set, denoted as a fourth data set; C) inputting the fourth data set as a reference to perform base quality score correction on the first data set of all samples; collecting the first data set of all samples after base quality score correction while retaining independent sample information, denoted as a fifth data set; collecting the second data set in all sample sets while retaining independent sample information, denoted as a sixth data set; D) performing joint mutation detection on the fifth data set and the sixth data set to obtain a candidate RNA editing site data set, denoted as a seventh data set; E) obtaining allele depth of each candidate RNA editing site in the seventh data set from the first data set of each sample to be detected for differential RNA editing sites, thereby obtaining allele depth table of each sample containing candidate RNA editing site information and editing site depth data, denoted as an eighth data set; F) combining the eighth data set of each sample to obtain a total allele depth table, denoted as a ninth data set; and performing statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth data set to determine differential RNA editing sites.

2. The method of claim 1, wherein, Step (A) specifically comprises: (A1) providing RNA sequencing data from N samples, and DNA sequencing data from the N samples, wherein N is a positive integer greater than or equal to 2; (A2) aligning the RNA sequencing data from N samples with a reference genome respectively to obtain N RNA alignment data and correct; and (A3) aligning the DNA sequencing data from N samples with a reference genome respectively to obtain N DNA alignment data and correct; wherein steps (A2) and (A3) can be performed interchangeably, sequentially or simultaneously.

3. The method of claim 1, wherein, In step (B), the following steps are further included: (B1) quality control to obtain data passing quality control; (B2) screening effective variations to obtain the known RNA editing site data set; (B3) data indexing to obtain the indexed known RNA editing site data set.

4. The method of claim 3, wherein, In step (B1), the maximum number of variant events allowed in each assembly region (max-events-in-region) is 4.

5. The method of claim 1, wherein, In step (D), further comprising steps: (D1) mutation site filtering; (D2) single nucleotide variant site screening; (D3) data processing.

6. The method of claim 5, wherein, In step (D3), further comprising steps: (D3a) homopolymer nucleotide sequence filtering; (D3b) genomic repeat sequence filtering.

7. The method of claim 1, wherein, Step (E) specifically comprises: (E1) splitting the first data set from N samples by positive and negative strands to obtain N positive strand subsets and N negative strand subsets; (E2) calculating the editing site depth data of each candidate RNA editing site in the seventh data set in the two subsets respectively.

8. An apparatus or system for detecting differential RNA editing sites, comprising: The device or system comprises: (M1) an input module configured to input RNA alignment data and DNA alignment data of N samples to be tested; (M2) a scanning module configured to perform the following operations: performing joint mutation detection on the RNA alignment data and DNA alignment data of each sample to be tested to obtain a known RNA editing site data set; inputting the known RNA editing site data set to perform base quality score correction on the RNA alignment data of each sample to be tested to obtain corrected RNA alignment data of each sample to be tested; (M3) a detection module configured to perform the following operations: collecting the corrected RNA alignment data of each sample to be tested to obtain a fifth data set; collecting the DNA alignment data of each sample to be tested to obtain a sixth data set; performing joint mutation detection on the fifth data set and the sixth data set to obtain a seventh data set; obtaining allele depth of each candidate RNA editing site in the seventh data set from the RNA alignment data of each sample to be tested to obtain an eighth data set; merging the eighth data set of each sample to be tested to obtain a ninth data set; statistically analyzing the candidate RNA editing site information and editing site depth data in the ninth data set to determine differential RNA editing sites; (M4) an output module configured to output information of the differential RNA editing sites.

9. The apparatus or system of claim 8, wherein, The scanning module (M2) comprises the following sub-modules: (M2.1) a preliminary detection sub-module configured to perform the following operations: performing joint mutation detection on the RNA alignment data and DNA alignment data of each sample to be tested to obtain a preliminary RNA editing site data set; (M2.2) a quality control sub-module configured to perform the following operations: performing quality control on the preliminary RNA editing site data set to obtain a quality-controlled RNA editing site data set; (M2.3) an effective variant screening sub-module configured to perform the following operations: screening by tags in the quality-controlled RNA editing site data set to obtain a known RNA editing site data set; (M2.4) a data indexing submodule, configured to index the known RNA editing site dataset, thereby obtaining an indexed known RNA editing site dataset; (M2.5) a quality correction submodule, configured to perform base quality score correction on the RNA alignment data, taking the indexed known RNA editing site dataset as input, thereby obtaining corrected RNA alignment data.

10. The apparatus or system of claim 8, wherein, The detection module (M3) comprises the following submodules: (M3.1) an official detection submodule, configured to perform joint mutation detection on the fifth dataset and the sixth dataset, thereby obtaining an official RNA editing site dataset; (M3.2) a mutation site filtering submodule, configured to perform mutation site-based quality filtering on the official RNA editing site dataset, thereby obtaining a filtered RNA editing site dataset; (M3.3) a screening submodule, configured to extract single nucleotide variant sites in the filtered RNA editing site dataset, thereby obtaining a screened RNA editing site dataset; (M3.4) a data processing submodule, configured to perform filtering on the screened RNA editing site dataset, thereby obtaining a candidate RNA editing site dataset, denoted as a seventh dataset; (M3.5) a depth sampling submodule, configured to obtain 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 an eighth dataset; and merging the eighth dataset of each of the test samples to obtain a ninth dataset (M3.6) a statistical analysis submodule, configured to perform statistical analysis on the candidate RNA editing site information and editing site depth data in the ninth dataset, thereby determining differential RNA editing sites.

Images