Method for quantitatively calculating RNA editing difference

Abstract

The invention relates to the technical field of molecular biology, and discloses a method for quantitatively calculating RNA (Ribonucleic Acid) editing difference, which comprises the following steps: after initial data is generated by high-throughput sequencing, carrying out mass filtering to remove low-mass readings; comparing the filtered data with a reference genome to generate a comparison result; carrying out base quality correction and variation site identification by utilizing GATK so as to identify an RNA editing site; then, carrying out quantitative analysis on the RNA editing event from the aspects of sites, genes and samples to obtain difference analysis data; carrying out difference analysis on the samples under different conditions, and calculating differences of RNA editing sites, genes and the samples to form an RNA editing difference matrix; finally, the difference matrix is analyzed through statistical examination, and RNA editing sites and genes with significant differences are screened out. The invention provides a systematic analysis framework, and RNA editing changes among different samples can be accurately evaluated and compared.

Description

Technical Field

[0001] The present invention relates to the technical field of molecular biology, and in particular to a method for quantitatively calculating RNA editing differences. Background Art

[0002] RNA editing is a biochemical process that is widely present in eukaryotes. It involves changes in the base sequence on the mRNA molecule. These changes are mainly caused by errors in the transcription process or post-transcriptional modifications. RNA editing is mainly mediated by deaminases, which can convert specific bases, such as cytosine (C) to uracil (U), or adenosine (A) to inosine (I). This base conversion can not only affect the splicing and stability of mRNA, but also affect the amino acid composition of proteins by changing the coding sequence, thereby regulating the structure and function of proteins. RNA editing can not only affect the coding function of proteins, but also involves editing of non-coding regions, changes in RNA structure, RNA shearing caused by editing of intron regions, and editing of 3' UTR and miRNA, targeting and regulating gene expression.

[0003] RNA editing plays a vital role in many biological processes, involving cell differentiation, neurotransmission, and immune response. In addition, abnormal RNA editing is closely related to the occurrence and development of various diseases, including cardiovascular disease, nervous system diseases, and cancer. For example, abnormal RNA editing of Apo-B protein is associated with lipoprotein metabolism disorders, while RNA editing of GluR2 receptors affects neural signal transmission. Therefore, RNA editing is of great significance for basic biological research, and it also has significant value in analyzing disease mechanisms and clinical treatment.

[0004] With the continuous advancement of high-throughput sequencing technology, researchers are now able to systematically identify RNA editing events across the entire genome. The development of these technologies has enabled RNA editing research to shift from the discovery of single genes to a comprehensive exploration of editing patterns across the transcriptome. By comparing RNA editing patterns under different conditions, researchers have been able to reveal the role of RNA editing in physiological and pathological processes. Despite this, compared to research areas such as DNA mutation, RNA expression, and RNA splicing, research in the field of RNA editing is still relatively small. Other research directions, such as the prediction of disease molecular markers and the study of disease progression, can also be applied to the field of RNA editing to reveal its potential value in precision medicine research.

[0005] Although some progress has been made, there is still a lack of an efficient and accurate method to quantify RNA editing differences. Existing methods are often limited by sequencing depth, the accuracy of identifying editing sites, and comparative analysis between different samples. Therefore, developing a new method that can accurately quantify RNA editing differences is crucial for a deeper understanding of the role of RNA editing in biology and disease. Summary of the invention

[0006] In view of this, the present invention proposes a method for quantitatively calculating RNA editing differences, aiming to provide a high-precision, high-sensitivity and high-repeatability RNA editing quantification technology to meet the needs of clinical sample analysis and basic research.

[0007] The present invention provides a method for quantitatively calculating RNA editing differences, comprising:

[0008] Perform high-throughput sequencing on samples to generate initial data;

[0009] Performing quality filtering on the initial data to remove low-quality reads;

[0010] Compare the filtered data with the reference genome to generate alignment results;

[0011] GATK was used to align the results, perform base quality correction, and identify variant sites to obtain RNA editing sites;

[0012] Quantitatively analyze RNA editing events at the site level, gene level, and sample level to generate corresponding differential analysis data;

[0013] Perform differential analysis on samples under different conditions, calculate the differences among RNA editing sites, genes, and samples, and generate an RNA editing differential matrix;

[0014] The difference matrix was analyzed by statistical tests to screen out RNA editing sites and genes with significant differences.

[0015] Furthermore, the RNA editing difference analysis includes site-level differences, gene-level differences, and sample-level differences;

[0016] Among them, site-level difference data include differences in site editing read numbers, site editing frequency, and RNA editing site load;

[0017] Gene-level difference data include differences in gene editing read numbers, gene editing frequency, and RNA editing gene load;

[0018] Sample-level difference data include differences in total editor load.

[0019] Further, the RNA editing site load is a value obtained by dividing the site editing frequency by the sequencing length and sequencing depth of the coding region;

[0020] The gene editing load is the value obtained by dividing the gene editing frequency by the sequencing length and sequencing depth of the gene coding region;

[0021] The total editing load is the total number of editing sites in a sample divided by the total length of the protein coding region in the sample.

[0022] Furthermore, the variance analysis was performed by using nonparametric statistical methods to ensure that the data did not rely on normal distribution;

[0023] At the site level, differential analysis was performed using the editing matrix of site editing read numbers and site editing frequency;

[0024] At the gene level, differential analysis was performed using the editing matrix of gene editing read numbers and gene editing frequency;

[0025] At the sample level, differential analysis was performed by total editing load of the samples.

[0026] Furthermore, the analysis results of RNA editing differences are displayed through visualization tools, including heat maps of editing site differences, difference volcano maps, and PCA cluster analysis maps.

[0027] Furthermore, the step of identifying the RNA editing site comprises:

[0028] GATK was used to calibrate the variant quality scores of the detected variant sites to filter out false positive sites;

[0029] Remove sites shared with the dbSNP database and exclude DNA mutation sites;

[0030] GATK was used to correct SNPs within 30 bp near indels to exclude false positive sites caused by indel alignment errors;

[0031] According to the criteria of alignment depth greater than 3 and mutation number greater than 3, mutation sites with less than 5 reads were filtered out; sites with mutation frequency less than 0.1 or greater than 0.9 were filtered out.

[0032] Furthermore, the RNA editing difference analysis also includes:

[0033] The differential RNA editing sites were functionally annotated by GO functional enrichment analysis, KEGG pathway enrichment analysis and GSEA enrichment analysis.

[0034] Furthermore, by performing association analysis between RNA editing differential genes and RNA expression differential genes, a nine-quadrant diagram was used to display the relationship between RNA editing differential genes and expression differential genes.

[0035] Compared with the prior art, the present invention has the following beneficial effects:

[0036] The present invention provides a more comprehensive and accurate RNA editing event analysis method, which can reveal the dynamic changes of RNA editing under different conditions or disease states, and provide a new perspective for understanding the role of RNA editing in biological processes and diseases. In addition, the present invention ensures the robustness of the analysis results through non-parametric statistical methods, and can accurately identify editing sites and genes with significant differences even when the data does not meet the normal distribution. Through functional annotation and association analysis, the present invention can also reveal the complex relationship between RNA editing and gene expression, and provide a powerful tool for in-depth study of the role of RNA editing in gene regulatory networks. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Various other advantages and benefits will become apparent to those of ordinary skill in the art by reading the detailed description of the preferred embodiments below. The accompanying drawings are only for the purpose of illustrating the preferred embodiments and are not to be considered as limiting the present invention. Moreover, the same reference symbols are used throughout the accompanying drawings to represent the same components. In the accompanying drawings:

[0038] Figure 1 This is the RNA editing clustering PCA diagram in Example 1;

[0039] Figure 2 This is the distribution diagram of RNA editing species in Example 1;

[0040] Figure 3 This is the RNA editing clustering PCA diagram in Example 1;

[0041] Figure 4 This is the RNA editing region distribution map in Example 1;

[0042] Figure 5 The volcano plot of RNA editing sites / gene differences in Example 1;

[0043] Figure 6 The heat map of RNA editing sites / gene differences in Example 1;

[0044] Figure 7 The differential GO function enrichment of RNA editing in Example 1;

[0045] Figure 8 The RNA editing difference KEGG functional enrichment in Example 1;

[0046] Fig. 9 The RNA editing difference GSEA functional enrichment in Example 1;

[0047] Fig.10 The nine-quadrant diagram of RNA editing differentially expressed genes and differentially expressed genes in Example 1;

[0048] Fig.11 The total RNA editing load (REB) of the sample in Example 1;

[0049] Fig.12 This is the RNA editing type distribution diagram in Example 2;

[0050] Fig.13 This is the distribution diagram of RNA editing species in Example 2;

[0051] Fig.14 This is the RNA editing region distribution map in Example 2;

[0052] Fig.15 The volcano plot of RNA editing sites / gene differences in Example 2;

[0053] Fig.16 The heat map of RNA editing sites / gene differences in Example 2;

[0054] Fig.17 The differential GO function enrichment of RNA editing in Example 2;

[0055] Fig.18 The RNA editing differential KEGG functional enrichment in Example 2;

[0056] Fig.19 GSEA functional enrichment of RNA editing differences in Example 2;

[0057] Fig. 20 The nine-quadrant diagram of RNA editing differentially expressed genes and differentially expressed genes in Example 2;

[0058] Fig.21 This is the total RNA editing load (REB) of the samples in Example 2. DETAILED DESCRIPTION

[0059] Exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments described herein. On the contrary, these embodiments are provided in order to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art. It should be noted that, in the absence of conflict, the embodiments of the present invention and the features described in the embodiments can be combined with each other. The present invention will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0060] The present invention provides an innovative RNA editing omics detection technology, which uses high-throughput sequencing data processing and comprehensive quantification strategies to achieve accurate quantification and analysis of RNA editing sites, genes and sample-level differences. The method covers efficient identification of RNA editing sites, editing quantification at the site and gene levels, calculation of total sample editing load, and functional annotation and visualization analysis of differential sites and genes. In view of the deficiencies of the prior art in quantifying and comparing RNA editing events, the present invention introduces unique quantitative indicators, such as the number of reads of site editing, editing frequency and RNA editing site load, and combines advanced statistical methods and bioinformatics tools to improve the accuracy and repeatability of RNA editing differential analysis. Through these technical means, the present invention not only promotes the application of RNA editing in the fields of biology and disease research, but also provides a key scientific basis and technical support for precision medicine and personalized treatment. The technical solution of the present invention is described below by specific embodiments.

[0061] Example 1

[0062] RNA editing and neuroinflammation

[0063] 1.1 Sample and data processing

[0064] This example case is intended to explore the potential mechanism of action of ropivacaine combined with dexamethasone in the sciatic nerve. In the experiment, 9 rats were randomly divided into three groups: Control group (control group), Rop group (ropivacaine group), and Rop+Dex group (ropivacaine and dexamethasone combined group), each group containing three rats. The Control group received a saline injection, the Rop group received a diluted ropivacaine (concentration of 0.5%), and the Rop+Dex group received a mixed solution injection of ropivacaine and dexamethasone. Two days after the injection, the rats were sampled and RNA-seq sequencing was performed to analyze their RNA editing site data and its annotation information. The statistics of the distribution of RNA editing types reflect the base changes at the editing site, among which C-to-U and A-to-I editing types are the most common. The statistics of the distribution of RNA editing types focus on the mutation type of the editing site, and the results show that most RNA editing sites are synonymous mutations, that is, although the RNA base has changed, it has not caused changes in the protein sequence. However, there are also some non-synonymous mutations, start codon deletions, stop codon deletions or gains, etc. These mutations directly affect the protein sequence and have significant biological significance. The statistics of RNA editing region distribution focus on the specific regions where the editing sites are located, especially the site changes in the exon region, because they directly affect the protein structure and are therefore the most intuitive to observe. The RNA editing cluster PCA diagram was analyzed using the prcomp R package and visualized using the ggplot2 R package. The specific RNA editing cluster PCA diagram is shown below: Figure 1 shown.

[0065] 1.2 RNA editing site quantification and differential analysis results

[0066] After obtaining RNA editing sites and their annotation information, we calculated the editing matrix and performed subsequent differential analysis according to different quantification strategies. The analysis results showed that although the differential analysis methods used varied depending on the quantification strategy, all differentially expressed genes met the conditions of a corrected p-value (padj) less than 0.05 and an absolute fold change (|FC|) greater than or equal to 2. The differential volcano plot was drawn using the ggplot2R package, while the differential heatmap was generated using the heatmap.2R package. The GOseq R package was used to identify the Gene Ontology (GO) terms that were significantly enriched in the RNA editing differentially expressed genes. KEGG (Kyoto Encyclopedia of Genes and Genomes) PATHWAY (http: / / www.genome.jp / kegg / ) is a major public pathway database for studying genetic biological functions. We used KOBAS (v2.0) software to perform KEGG pathway enrichment analysis on RNA editing differentially expressed genes, where GO terms and KEGG pathways that met padj less than 0.05 were considered to be significantly enriched. Gene Set Enrichment Analysis (GSEA) is a method used to study whether the expression of a predefined set of genes changes in a coordinated manner under a specific biological state. We used the fgseaR package to perform GSEA enrichment analysis on RNA editing differential genes, where the p-value was less than 0.05 and the normalized enrichment score (|Nes|) was greater than 1. Figure 2-4 They are the RNA editing type distribution map, RNA editing cluster PCA map and RNA editing region distribution map. Figure 5-9 They are RNA editing site / gene difference volcano map, RNA editing site / gene difference heat map, RNA editing difference GO functional enrichment, RNA editing difference KEGG functional enrichment and RNA editing difference GSEA functional enrichment.

[0067] The nine-quadrant plot is an effective tool for revealing the correlation between differentially edited RNA genes and differentially expressed genes. It helps us identify and interpret patterns and trends in the data set. When processing the data, we used the dplyr and ggplot2 packages in the R language to merge and organize the information. Then, using the geom_point function in the ggplot2 package, we drew a scatter plot and labeled different data points by distinguishing between colors and shapes. The nine-quadrant plot of differentially edited RNA genes and differentially expressed genes is shown in Figure 1. Fig.10 shown.

[0068] The total editing burden (REB) of each sample was calculated by the following formula: the total number of editing sites divided by the total length of the protein coding region (in million base pairs). Next, we used a t-test to assess whether there were significant differences in REB between different groups. However, since there were only three samples in each group in the neuroinflammation study and the individual differences were small, no significant differences in the total RNA editing burden between different groups were detected. The total RNA editing burden (REB) of the samples is shown in the figure. Fig.11 shown.

[0069] Example 2 RNA Editing and Systemic Lupus Erythematosus

[0070] 2.1 Sample and data processing

[0071] This case study is dedicated to exploring the mechanism of RNA editing in systemic lupus erythematosus (SLE). We analyzed the CD19 + B cell samples were sequenced and analyzed, and the samples were divided into three groups: AC group (healthy control group, 92 cases), AS group (SLE group without renal damage, 33 cases) and AL group (SLE group with renal damage, 47 cases). Through analysis, we obtained the RNA editing site data and annotation information of these groups. The RNA editing type distribution statistics revealed the base changes at the editing site, mainly C-to-U and A-to-I editing. In the RNA editing type distribution statistics, we found that most editing sites were synonymous mutations, which means that the change of RNA base did not lead to changes in protein sequence. However, there are also some non-synonymous mutations, start codon deletions, stop codon deletions or gains, etc. These mutations directly affect the protein sequence and have significant biological significance. The RNA editing region distribution statistics focus on the gene region where the editing site is located, especially the changes in the exon region. These changes are directly related to the protein structure and are therefore the most intuitive in observation. The distribution of RNA editing shows subtle differences in samples from different groups. Figure 12-14 They are the RNA editing type distribution map, RNA editing species distribution map and RNA editing region distribution map.

[0072] 2.2 RNA editing site quantification and differential analysis results

[0073] After obtaining RNA editing sites and their annotation information, we calculated the editing matrix and performed subsequent differential analysis according to different quantification strategies. The analysis results showed that although the differential analysis methods used varied depending on the quantification strategy, all differentially expressed genes met the conditions of a corrected p-value (padj) less than 0.05 and an absolute fold change (|FC|) greater than or equal to 2. The differential volcano plot was drawn using the ggplot2R package, while the differential heatmap was drawn using the heatmap.2R package. The GOseq R package was used to identify the Gene Ontology (GO) terms that were significantly enriched in the RNA editing differentially expressed genes. KEGG (Kyoto Encyclopedia of Genes and Genomes) PATHWAY (http: / / www.genome.jp / kegg / ) is a major public pathway database for studying genetic biological functions. We used KOBAS (v2.0) software to perform KEGG pathway enrichment analysis on RNA editing differentially expressed genes, where GO terms and KEGG pathways that met padj less than 0.05 were considered to be significantly enriched. Gene Set Enrichment Analysis (GSEA) is a method used to study whether the expression of a predefined gene set changes in a coordinated manner under a specific biological state. We used the fgseaR package to perform GSEA enrichment analysis on RNA editing differential genes, where gene sets with a p-value less than 0.05 and a standardized enrichment score (|Nes|) greater than 1 were considered significantly enriched. Figure 15-19 They are RNA editing site / gene difference volcano map, RNA editing site / gene difference heat map, RNA editing difference GO functional enrichment, RNA editing difference KEGG functional enrichment and RNA editing difference GSEA functional enrichment.

[0074] The nine-quadrant plot is an effective tool for revealing the correlation between differentially edited genes and differentially expressed genes. It helps us identify and interpret patterns and trends in the data set. When processing the data, we used the dplyr and ggplot2 packages in the R language to merge the data and used the geom_point function in the ggplot2 package to draw a scatter plot. By distinguishing between colors and shapes, we can clearly show the characteristics of different data points. Fig. 20 This is a nine-quadrant diagram of RNA editing differentially expressed genes and differentially expressed genes.

[0075] The total editing burden (REB) of each sample was calculated by the following formula: the total number of editing sites divided by the total length of the protein coding region (in megabase pairs, mb). Then, a t-test was used to evaluate whether there was a significant difference in REB between different groups. The analysis results showed that there was a significant difference in REB between the AC group and the AS group. Fig.21 It is the total RNA editing burden (REB) of the sample.

[0076] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, ordinary technicians in the relevant field should understand that the specific implementation methods of the present invention can still be modified or replaced by equivalents. Any modification or equivalent replacement that does not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A method for quantitatively calculating RNA editing differences, characterized in that: include: Perform high-throughput sequencing on samples to generate initial data; Performing quality filtering on the initial data to remove low-quality reads; Compare the filtered data with the reference genome to generate alignment results; GATK was used to align the results, perform base quality correction, and identify variant sites to obtain RNA editing sites; Quantitatively analyze RNA editing events at the site level, gene level, and sample level to generate corresponding differential analysis data; Perform differential analysis on samples under different conditions, calculate the differences among RNA editing sites, genes, and samples, and generate an RNA editing differential matrix; The difference matrix was analyzed by statistical tests to screen out RNA editing sites and genes with significant differences.

2. The method for quantitatively calculating RNA editing differences according to claim 1, characterized in that The RNA editing difference analysis includes site-level differences, gene-level differences, and sample-level differences; Among them, site-level difference data include differences in site editing read numbers, site editing frequency, and RNA editing site load; Gene-level difference data include differences in gene editing read numbers, gene editing frequency, and RNA editing gene load; Sample-level difference data include differences in total editor load.

3. The method for quantitatively calculating RNA editing differences according to claim 2, characterized in that: The RNA editing site load is the value obtained by dividing the site editing frequency by the sequencing length and sequencing depth of the coding region; The gene editing load is the value obtained by dividing the gene editing frequency by the sequencing length and sequencing depth of the gene coding region; The total editing load is the total number of editing sites in a sample divided by the total length of the protein coding region in the sample.

4. The method for quantitatively calculating RNA editing differences according to claim 2, characterized in that: The variance analysis was performed by using nonparametric statistical methods to ensure that the data did not rely on normal distribution; At the site level, differential analysis was performed using the editing matrix of site editing read numbers and site editing frequency; At the gene level, differential analysis was performed using the editing matrix of gene editing read numbers and gene editing frequency; At the sample level, differential analysis was performed by total editing load of the samples.

5. The method for quantitatively calculating RNA editing differences according to claim 1, characterized in that The analysis results of RNA editing differences are displayed through visualization tools, including heat maps of editing site differences, difference volcano maps, and PCA cluster analysis maps.

6. The method for quantitatively calculating RNA editing differences according to claim 1, characterized in that: The RNA editing site identification step comprises: GATK was used to calibrate the variant quality scores of the detected variant sites to filter out false positive sites; Remove sites shared with the dbSNP database and exclude DNA mutation sites; GATK was used to correct SNPs within 30 bp near indels to exclude false positive sites caused by indel alignment errors; According to the criteria of alignment depth greater than 3 and mutation number greater than 3, mutation sites with less than 5 reads were filtered out; sites with mutation frequency less than 0.1 or greater than 0.9 were filtered out.

7. The method for quantitatively calculating RNA editing differences according to claim 1, characterized in that: The RNA editing differential analysis also includes: The differential RNA editing sites were functionally annotated by GO functional enrichment analysis, KEGG pathway enrichment analysis and GSEA enrichment analysis.

8. The method for quantitatively calculating RNA editing differences according to claim 1, characterized in that: By performing association analysis between RNA editing differential genes and RNA expression differential genes, a nine-quadrant diagram was used to display the relationship between RNA editing differential genes and expression differential genes.

Images