A model for identifying asm6a modifications and uses thereof

Abstract

The application discloses a model for identifying ASm6A modification, which is based on MeRIP-seq sequencing data of a to-be-tested sample and can accurately identify ASm6A modification peaks in the to-be-tested sample. When the model identifies the ASm6A modification peaks in the to-be-tested sample, the error rate is always lower than 10% no matter what MeRIP-seq sequencing data of the to-be-tested sample is used, which indicates that the model has excellent stability in identifying the ASm6A modification of the sample and is suitable for identifying MeRIP-seq sequencing data of the sample under different conditions. Compared with the prior art for identifying the ASm6A modification, the AUC value is obviously improved when the model is used to identify the ASm6A modification in the sample, and the identification effect is obviously better.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of biomedical technology, in particular, to a model for identifying ASm6A modification and application thereof. BACKGROUND

[0002] N6-methyladenosine (m6A) is the most common dynamic reversible RNA modification in eukaryotic cells. Through the joint action of a series of methylation enzymes, demethylation enzymes and methylation recognition proteins, m6A modification forms a precise dynamic regulatory cycle in the body, and is widely involved in embryonic development, cell apoptosis, sperm development, immunity and stress response and other life activities. Abnormal m6A modification can cause serious body pathology. For example, existing technologies find that hypoxia can induce breast cancer stem cell phenotype through m6A demethylation. In addition, m6A can also be used as a target for inhibiting cancer under certain conditions. Therefore, identifying the precise site of m6A and studying its internal regulation mechanism are the key basis for elucidating the pathogenesis of related diseases and developing new therapeutic targets.

[0003] At present, due to the resolution limit of experimental methods, it is still very difficult to analyze the precise regulation mechanism of m6A modification at the omics level. However, the method of immunoprecipitation combined with new generation sequencing technology (MeRIP-seq) can identify tens of thousands of RNA methylation modification signal peaks with a length of 100-200 nt at the transcriptome level, which plays an important role in revealing the physiological function of RNA methylation modification.

[0004] At present, with the proposal of the method of ultraviolet light cross-linking m6A antibody combined with immunoprecipitation (m6A individual-nucleotide resolution cross-linking and immunoprecipitation, miCLIP), it is possible to identify single-base m6A modification sites from the whole transcriptome at the experimental level. Existing technologies have high-throughput identified m6A modification sites of mammals at single-base precision, but the results obtained by this technology are still not very stable, and it has not been widely used in the field of m6A research. Therefore, improving the quality of the results of identifying m6A modification regions from MeRIP-seq data will be the focus of current m6A research.

[0005] On the basis of accurately identifying the m6A modification region of the transcriptome, another important research purpose is to analyze the allele-specific modification time. In non-haploid biological clock, the expression of certain modification sites or certain genes tends to occur on a specific chromosome, which is called allele-specificity. Some important life activities, such as gene transcription and epigenetic modification, have significant allele-specificity. These allele-specific regulations affect cell life activities in many ways, including imprinting regulation of embryonic development, inactivation of chromosomes, and regulation of gene expression in specific temporal and spatial contexts. Disturbed allele-specific regulation events often cause serious body disorders, such as abnormal biological development, cardiovascular and cerebrovascular dysfunction, and even cancer. Therefore, the identification and research of allele-specific regulation events can help researchers better understand the mechanisms of various diseases and explore potential therapeutic targets.

[0006] However, due to technical limitations, for a long period of time, the study of allele-specific regulation events was limited to a limited number of genes or chromosome segments. With the development of high-throughput sequencing technology, the study of allele-specific regulation events has been carried out at the omics level. Sequencing technology can identify single nucleotide polymorphism sites (SNPs) present on DNA or RNA. In non-haploid organisms, chromosomes from the mother and father generally carry different SNP sites, so the source of sequencing reads can be inferred by heterozygous SNP sites on the sequence, and then the regulation signals from different alleles can be identified.

[0007] Allele-specific modification is considered to be one of the important reasons for the generation of allele-specific expression. Existing researches mainly focus on the DNA level. For example, when the existing technology studies mouse brain tissue, it is found that more than 1300 alleles have significant differences in DNA methylation levels on two homologous chromosomes. In addition, the existing technology also reveals that more than 10% of human genes are regulated by allele-specific DNA methylation. These allele-specific methylation modifications have a high correlation with SNP sites on the sequence, and the current research on allele-specificity on RNA modification is relatively less.

[0008] Compared with the above-mentioned allele-specific regulation time, the research on the allele specificity of m6A modification on the transcriptome has not been reported. The m6A modification controls the correct function of many metabolism and growth and differentiation related genes in the body. The prior art finds that knockout of METTL3 in mouse embryonic stem cells can prevent differentiation; in Drosophila, knockout of the homologous protein Ime4 of METTL3 also has a sublethal effect, which causes damage to the NOTCH signaling pathway and affects the fertility of Drosophila; in addition, the Ime4 gene also plays an important role in the meiosis of yeast. All of this shows that m6A modification is essential for the early development of sexual organs and the early development of embryos, so studying whether m6A modification has specificity between alleles is an important entry point for deciphering the mechanism of embryonic development.

[0009] The m6A modification signal is mainly analyzed from MeRIP-seq sequencing data by relying on computational biology methods. Researchers first used Fisher's Exact Test method to construct the first peak identification algorithm for MeRIP-seq, which can achieve high identification specificity, but it ignores the upstream and downstream correlation of m6A signals on transcripts, and lacks accurate modeling of read segments on the whole genome in steps, so the accuracy and sensitivity are insufficient. Exome Peak, HEPeak and MetPeak algorithms mainly use Beta binomial distribution to model the methylation level of each chromosome window, and introduce a Hidden Markov Model (HMM) to calculate the upstream and downstream correlation of methylation events. However, these algorithms must first assume that the MeRIP-seq reads conform to the Beta binomial distribution, so obvious errors will be introduced for samples with low sequencing quality, and the parameter optimization of the HMM model requires a certain number of repeated samples to obtain relatively accurate results, but general MeRIP-seq sequencing will only be repeated a small number of times, which also limits the identification performance of the algorithm to some extent. In addition, researchers also use MACS to identify m6A signals in MeRIP-seq sequencing data, but MACS is a peak identification software specifically designed for ChIP-seq analysis, and its application in MeRIP-seq sequencing analysis with different principles will produce obvious errors, and there is also a coordinate offset.

[0010] Therefore, there is an urgent need for a method that can accurately identify m6A peak signals based on MeRIP-seq sequencing, and then accurately identify allele-specific m6A. SUMMARY

[0011] The model for identifying ASm6A modification and application thereof.

[0012] The first object of the present application is to provide a model for identifying ASm6A modification.

[0013] The second object of the present application is to provide application of the above model in preparation of a product for identifying ASm6A.

[0014] The third object of the present application is to provide an ASm6A modification identification system.

[0015] The fourth object of the present application is to provide a computer device.

[0016] The fifth object of the present application is to provide a computer readable storage medium.

[0017] The sixth object of the present application is to provide a computer program product.

[0018] The seventh object of the present application is to provide a model for identifying differential ASm6A of paired samples.

[0019] In order to achieve the above objects, the present application is realized by the following scheme:

[0020] A model for identifying ASm6A modification, comprising a data acquisition module, a data processing module, a modification identification module and a result output module.

[0021] The data acquisition module is used to acquire MeRIP-seq sequencing data of a to-be-tested sample.

[0022] The data processing module acquires m6A peaks of the to-be-tested sample based on the MeRIP-seq sequencing data of the to-be-tested sample, acquires SNP site information covered by each m6A peak, and acquires SNP site information covered by each m6A peak. j (m) According to formula 10, the odds ratio ρ is calculated j (m) And the logarithmic odds ratio y is obtained by taking the logarithm j (m) ;

[0023] Formula 10:

[0024] Where n j (m) is the read count of SNP j (m) in the to-be-tested sample; x (m) ma,j is the read count of SNP j (m)The number of readings of the main haplotype;

[0025] μ ASE The method for obtaining the SNP sites is as follows: based on the RNA-seq sequencing data of the sample to be tested, the SNP sites of the sample to be tested are obtained, and the SNPs in the sample to be tested are targeted. j And calculate SNP according to Formula I. j The dominance ratio ρ of the major haplotype relative to the minor haplotype j And take the logarithm to get the logarithmic dominance ratio y j ;

[0026] Formula 1:

[0027] Where n j SNPs in the sample to be tested j The number of reads; x ma,j SNP j Number of readings for the primary haplotype; x 0,j 0.5*n j ;

[0028] Log dominance ratio y j By combining the processing with a one-way normal random distribution, we obtain the processed y. j For the processed y j Perform sampling and calculate the mean μ of the posterior distribution. ASE ;

[0029] The log dominance ratio of the SNP sites covered by each m6A peak was calculated. j (m) By combining the processing with a one-way normal random distribution, we obtain the processed y. j (m) ; Regarding the processed y j (m) Sampling is performed and the mean μ of the posterior distribution is calculated. (m) , will μ (m) Taking the antilogarithm gives ρ j (m) ', and calculate the frequency of the major allele for each m6A peak;

[0030] The difference significance analysis module is used to perform difference significance analysis on the frequencies of the major alleles of each m6A peak obtained by the data processing module, and obtain the difference significance analysis results of each m6A peak;

[0031] The modification identification module shown identifies m6A peaks with significant differences as ASm6A modified peaks based on the analysis results of the differences in each m6A score, and identifies m6A peaks without significant differences as non-ASm6A modified peaks.

[0032] The result output module is used to output the results obtained by the modification recognition module.

[0033] Among them, ASm6A is an allele-specific RNA N6-methyladenosine (Allele-specofoc m6A, ASm6A) modification, as detailed in the prior art: Genome Res. 2023 Aug; 33(8): 1369-1380. doi: 10.1101 / gr.277704.123. Epub 2023 Sep 15.

[0034] Preferably, in the data processing module, for the processed y j (m) Sampling is performed using the Metropolis-Hastings sampling method, which is an improvement on the Markov chain Monte Carlo algorithm.

[0035] Preferably, in the data processing module, the SNP j The main haplotype is SNP. j The base with the highest base reading, the SNP j The last haplotype was SNP. j The bases of the last high-base reading.

[0036] Preferably, in the data processing module, the frequency of the major allele of each m6A peak is calculated as follows: frequency of major allele = {ρ j (m) / (1+ρ j (m) )}.

[0037] This invention also claims protection for the use of any of the models described above in the preparation of products for identifying ASm6A.

[0038] The present invention also claims protection for an ASm6A modification recognition system, comprising an acquisition unit, a storage unit, and a processing unit, wherein the acquisition unit is used to acquire MeRIP-seq sequencing data of a sample to be tested;

[0039] The storage unit stores program instructions that can be executed by the processing unit;

[0040] The processing unit includes any of the models described above;

[0041] When the program instructions are executed by the processing unit, the MeRIP-seq sequencing data of the sample to be tested obtained by the acquisition unit is input into any of the models described above to obtain the ASm6A modification recognition result.

[0042] The present invention also claims protection for a computer device, including a memory and a processor, wherein the memory stores a computer program executable on the processor;

[0043] When the computer program is executed by the processor, it implements any of the models described above and / or the ASm6A modification recognition system described above.

[0044] The present invention also claims protection for a computer-readable storage medium having a computer program stored thereon, characterized in that the computer program, when executed by a processor, implements any of the models described above and / or the ASm6A modification recognition system described above.

[0045] The present invention also claims protection for a computer program product comprising the aforementioned computer-readable storage medium.

[0046] The present invention also claims protection for a model for identifying differential ASm6A in paired samples, characterized in that it includes a data acquisition module, an identification module, and a result output module;

[0047] The paired samples are samples A and B from different groups but originating from the same individual;

[0048] The data acquisition module is used to acquire MeRIP-seq sequencing data of paired samples and identify them according to any of the models described above, so as to obtain the ASm6A modification peaks of sample A and sample B respectively.

[0049] The identification module, based on the MeRIP-seq sequencing data of sample A and sample B obtained by the data acquisition module, records ASm6A modification peaks with more than 50% overlap in genomic coordinates between sample A and sample B as m6A modification peaks that coexist in both sample A and sample B.

[0050] For ASm6A modification peaks that exist in sample A but not in sample B, they are denoted as specific ASm6A modification peaks in sample A.

[0051] For ASm6A modification peaks that are absent in sample A but present in sample B, they are denoted as nonspecific ASm6A modification peaks of sample A.

[0052] For ASm6A modification peaks that are present in both sample A and sample B but show different major haplotypes, they are denoted as the specific ASm6A modification peak of sample A.

[0053] For ASm6A modification peaks that are present in both sample A and sample B but show the same major haplotype, the difference between ASm6A modification peaks in sample A and sample B is determined by hierarchical Bayesian model. ASm6A modification peaks with significant differences are recorded as specific ASm6A modification peaks of sample A, and ASm6A modification peaks with no significant differences are recorded as non-specific ASm6A modification peaks of sample A.

[0054] The result output module is used to output the results obtained by the identification module.

[0055] Preferably, the difference is significant when q value < 0.05, and the difference is not significant when q value ≥ 0.05.

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

[0057] This invention provides a model for identifying ASm6A modifications. Based on MeRIP-seq sequencing data of the sample, the model accurately identifies ASm6A modification peaks in the sample. When identifying ASm6A modification peaks in the sample, the model consistently maintains an error rate below 10% regardless of the type of MeRIP-seq sequencing data used, demonstrating excellent stability in identifying ASm6A modifications and applicability to MeRIP-seq sequencing data under different conditions. Furthermore, compared to existing methods for identifying ASm6A modifications, the model of this invention significantly improves the AUC value, resulting in significantly better identification performance. Attached Figure Description

[0058] Figure 1 This is a schematic diagram of the principle of a model for ASm6A modification recognition, as shown in Example 2.

[0059] Figure 2 The graph shows the average error rate results when identifying datasets for each classification category according to the method shown in Example 3; A is the average error rate result for the dataset classified by sequencing read length; B is the average error rate result for the dataset classified by library size; C is the average error rate result for the dataset classified by FPKM value of gene expression; D is the average error rate result for the dataset classified by the number of SNPs in the m6A peak; E is the average error rate result for the dataset classified by the number of biological repeats.

[0060] Figure 3 The graph shows the average AUC results for the experimental group, control group 1, and control group 2 in Example 4.

[0061] Figure 4The figures below show the area under the ROC curves for the experimental group, control group 1, and control group 2 in Example 4; A is the AUC result for identifying datasets classified by the number of SNPs in the m6A peak; B is the AUC result for identifying datasets classified by the FPKM value of gene expression; C is the AUC result for identifying datasets classified by the number of biological repeats; D is the AUC result for identifying datasets classified by sequencing read length; and E is the AUC result for identifying datasets classified by library size.

[0062] Figure 5 This is a visualization of the IGV results from Example 4;

[0063] Figure 6 The image shows the identification results of the differential ASm6A in Example 5. Detailed Implementation

[0064] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. These embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. Unless otherwise specified, the experimental methods used in the following embodiments are conventional methods; the materials and reagents used, unless otherwise specified, are commercially available.

[0065] Example 1: A functional model for identifying genes with significant allele-specific expression.

[0066] The functional model for identifying genes with significant allele-specific expression includes a data acquisition module, a data processing module, a differential significance analysis module, a prediction module, and a result output module.

[0067] The data acquisition module is used to acquire RNA-seq sequencing data of the sample to be tested.

[0068] The data processing module performs the following processing based on the RNA-seq sequencing data of the sample to be tested, specifically including:

[0069] Based on the RNA-seq sequencing data of the sample to be tested, the RNA-seq sequencing data of the sample to be tested is mapped to the corresponding genome of the sample to be tested, and the SNP site information of each gene in the sample to be tested and the number of reads of each SNP are obtained. The base with the highest base reading at a SNP site is recorded as the major haplotype of the site, and the base with the second highest base reading is recorded as the minor haplotype of the site.

[0070] For SNP j Calculate SNPs according to Formula 1 j The dominance ratio ρ of the major haplotype relative to the minor haplotype jAnd according to Formula 2, take the logarithm of the advantage ratio to obtain the logarithmic advantage ratio y. j Among them, SNP j These are the SNP sites in the sample to be tested;

[0071] Formula 1:

[0072] Formula 2:

[0073] Where n j SNPs in the sample to be tested j The number of reads; x ma,j SNP j The number of major haplotype readings; x 0,j 0.5*n j ;

[0074] The log-dominance ratio y obtained from Formula 2 j Combining the single-item normal random effects model (REM) with formulas 3 to 6, the processed y is obtained. j ;

[0075] Formula 3:

[0076] Formula 4: θ j ~N(μ,τ) 2 );

[0077] Formula 5: μ ~ Uniform(-∞, +∞);

[0078] Formula 6:

[0079] Where θ j σ is the mean of the log-odds ratios of the sample to be tested. j 2 This represents the variance of the log-odds ratio of the sample to be tested. s 2 =10.

[0080] For the processed y j The Metropolis-Hastings (MH) sampling method, an improvement on the Markov Chain Monte Carlo (MCMC) algorithm, is used for sampling according to Equation 7. Marginal probabilities are calculated according to Equation 8, followed by the posterior distribution of μ according to Equation 9, and finally the mean of the posterior distribution of μ (μ). ASE ).

[0081] Formula 7: p(θ1, …, θ) j ,μ,τ|y1,y2,…,y j );

[0082] Formula 8: p(μ,τ|y1,y2,…,y j );

[0083] Formula 9: p(τ|y1,y2,…,y j ).

[0084] The difference significance analysis module is based on the mean (μ) of the μ posterior distribution obtained by the data processing module. ASE A significance analysis of the differences was conducted, as follows:

[0085] (1) Construct tail distribution data of major allele frequencies (MAF) under the null hypothesis.

[0086] To distinguish significant allele-specific events, the major allele frequency (MAF) distribution needs to be obtained under null hypothesis conditions. Due to the lack of standard real MeRIP-seq data, it is necessary to simulate sequencing data to obtain read counts of major and minor alleles of SNPs without significant ASE or ASm6A events. The read count distribution of individual SNPs captured by sequencing is fitted by using a negative binomial distribution (NBD).

[0087] Assuming coverage of each SNP site (SNP) i The number of reads for a single allele is shown in Equation 10;

[0088] Formula 10:

[0089] Where μ represents the SNP without allelic-specific events. i The theoretical number of segments read is equal to 0.5N. i N i SNP i The number of readscount; k is the discreteness parameter.

[0090] The dispersion parameter k was obtained as follows: whole-genome sequencing (WGS) data were obtained from the 1000 Genomes database (https: / / www.internationalgenome.org / ), and the number of reads N covering each SNP site in each sample was counted. i The number of allele reads x and N were integrated using the maximum likelihood method. i And x, we can obtain the discreteness parameter k.

[0091] Next, FPKM values ​​for all genes were collected from the Cancer Genome Atlas (TCGA) project, and the fitter package in Python was used to fit the data. Genes were categorized into six classes based on their FPKM distribution type (as described in existing works: de Torrente L, Zimmerman S, Suzuki M, Christopeit M, Greally JM, Mar JC. The shape of gene expression distributions matter: how incorporating distribution shape improves the interpretation of cancer transcriptomic data. BMC Bioinformatics. 2020; 21:562). The overall distribution of FPKM for each class was refitted, and samples were drawn from these distributions to obtain the simulated FPKM for each gene. Based on the simulated FPKM and gene length, the total number of reads N for each SNP on that gene was calculated. i '.

[0092] Based on the discreteness parameter k and the simulated N i ', combined with Formula 10, to obtain the SNP of each metaunit (gene) i The negative binomial distribution of allele read numbers was used, and the number of reads of major and minor alleles was simulated by sampling. The dominance ratio ρ was calculated according to Formula 1. i ', and calculate the MAF for each gene; MAF = {ρ i ' / (1+ρ i ')}.

[0093] By introducing the generalized Pareto distribution (GPD), the tail distribution of MAF' for different types of genes was estimated, and the statistical significance threshold was calculated. The 94th percentile was used as the initial threshold t. For MAFs greater than the initial threshold t, the scale parameter and shape parameter of GPD in Formula 11 were combined to estimate their distribution, and the tail distribution of MAF values ​​in the null hypothesis data was obtained.

[0094] Formula 11:

[0095] The parameter ψ is estimated using a hybrid method and gradient descent algorithm proposed in the existing technology (Wang C, Chen GA, new hybrid estimation method for the generalized pareto distribution. Communications in Statistics-Theory and Methods. 2016; 45: 4285-4294.).

[0096] (2) Analysis of significant differences

[0097] The mean (μ) of the μ posterior distribution obtained by the data processing module ASE Taking its antinomial gives ρ ASE And calculate μ ASE The corresponding MAF value (R); MAF = {ρ ASE / (1+ρ ASE )}.

[0098] μ ASE The corresponding MAF value(R) was used in conjunction with Formula 12 to perform a significance analysis of the difference, and the results of the significance analysis of the difference were obtained.

[0099] Formula 12:

[0100] Among them, t is 94%, R simj This indicates that the MAF value for ASm6A / ASE does not exist in the null hypothesis data, N sim Indicating R in the simulation data simj total.

[0101] The value of P is corrected using the BH method to obtain the q value. A q value less than 0.05 is considered to be statistically significant.

[0102] The prediction module obtains the difference significance analysis results based on the difference significance analysis module, and determines whether the gene is an ASE gene. Specifically, it selects genes with significant differences (i.e., q value < 0.05) as the basis for its prediction. ASE The gene corresponding to the corresponding MAF value was identified as the ASE gene; μ values ​​that did not show significant differences (i.e., q value ≥ 0.05) were considered. ASE The gene corresponding to the MAF value is determined to be a non-ASE gene.

[0103] The result output module is used to output the result obtained by the judgment module.

[0104] Example 2: A model for ASm6A modification recognition

[0105] A schematic diagram of the principle of a model for ASm6A modification recognition is shown below.Figure 1 As shown.

[0106] The model for ASm6A modification identification includes a data acquisition module, a data processing module, a difference significance analysis module, a modification identification module, and a result output module.

[0107] The data acquisition module is used to acquire MeRIP-seq sequencing data of the sample to be tested.

[0108] The data processing module performs the following processing based on the MeRIP-seq sequencing data of the sample to be tested:

[0109] Based on the MeRIP-seq sequencing data of the sample to be tested, the m6A peak in the sample to be tested was obtained, and the SNP site information covered by each m6A peak was obtained.

[0110] SNPs covered by the m6A peak j (m) Calculate its dominance ratio ρ according to Formula 10. j (m) And according to Formula 11, take the logarithm of the advantage ratio to obtain the logarithmic advantage ratio y. j (m) ;

[0111] Formula 10:

[0112] Formula 11:

[0113] Where n j (m) SNPs in the sample to be tested j (m) The number of reads; x (m) ma,j SNP j (m) The number of readings of the main haplotype; μ ASE SNPs were treated using a functional model for identifying genes with significant allele-specific expression, as shown in Example 1. i (m) Obtain;

[0114] The log-dominance ratio y obtained from Formula 11 j (m) Combining the single-term normal random distribution (REM) with Equations 12 and 15, we obtain the processed y. j (m) ;

[0115] Formula 12:

[0116] Formula 13:

[0117] Formula 14: μ (m) ~Uniform(-∞,+∞);

[0118] Formula 15:

[0119] Where θ j σ is the mean of the log-dominance ratio. j 2 The variance of the log-dominance ratio; s 2 =10.

[0120] For the processed y j (m) The Metropolis-Hastings (MH) sampling method, which is an improvement on the Markov Chain Monte Carlo (MCMC) algorithm, is used for sampling, and the mean μ of the posterior distribution is calculated. (m) , will μ (m) Taking the antilogarithm gives ρ j (m) '(i.e., ρ) j (m) The expected value of the m6A peak was calculated, and the frequency of the major allele (MAF) was calculated.

[0121] Where MAF={ρ j (m) / (1+ρ j (m) )}.

[0122] The difference between the difference significance analysis module and the difference significance analysis module in Example 1 is that the difference significance analysis is performed on the major allele frequency (MAF) on the m6A peak according to the difference significance analysis in (2) of the difference significance analysis module in Example 1, and the difference significance analysis result (q value) of the m6A peak is obtained.

[0123] The modification identification module, based on the difference significance analysis results of the m6A peak obtained by the difference significance analysis module, identifies ASm6A modified peaks; m6A peaks with significant differences (i.e., q value < 0.05) are identified as ASm6A modified peaks; m6A peaks without significant differences (i.e., q value ≥ 0.05) are identified as non-ASm6A modified peaks.

[0124] The result output module is used to output the results obtained by the modification recognition module.

[0125] Example 3: A method for identifying differences ASm6A in paired samples

[0126] Paired samples are samples from different groups but originating from the same individual, such as tumor samples and normal samples from the same patient. Taking tumor samples and normal samples from the same patient as an example, the specific method is as follows:

[0127] m6A peaks with more than 50% overlap in genomic coordinates between tumor and normal samples were considered to originate from the same m6A peak (i.e., m6A peaks present in both tumor and normal samples). The tumor and normal samples in the paired samples were processed using the model shown in Example 2 to obtain the ASm6A modification peaks of the tumor and normal samples, respectively, and then classified.

[0128] (1) For ASm6A peaks that are present in tumor samples but not in normal samples, they are denoted as specific ASm6A modification peaks in tumor samples.

[0129] (2) For ASm6A peaks that exist in normal samples but not in tumor samples, they are denoted as non-specific ASm6A modification peaks of tumor samples.

[0130] (3) ASm6A peaks that are present in both tumor and normal samples but show different major haplotypes are denoted as specific ASm6A modification peaks of tumor samples.

[0131] (4) For ASm6A peaks present in both tumor and normal samples and showing the same major haplotype, a hierarchical Bayesian model was used to assess the difference. Consistent heterozygous SNP sites within the overlapping region of the ASm6A peak were considered as usable sites for downstream analysis to ensure consistent haplotype construction between the two samples. For each SNP site, the dominance ratio ρ was calculated according to Formula 16. j s ;

[0132] Formula 16:

[0133] in y tumor，j SNPs in tumor IP samples j Read count of the major allele at the locus, n tumor，j SNPs in tumor IP samples j The number of reads counts at a site, μ tumor，b The expected value of the major allele log dominance ratio for allele-specific expression calculated from tumor input samples (i.e., general RNA-seq), where ρ under the null hypothesis... tumor，j =ρ normal，j ;

[0134] Using ρ j sCalculate the corresponding MAF, and perform a significance analysis of the difference according to the method shown in the significance analysis module of Example 1 to obtain the significance analysis result q value; MAF={ρ j s / (1+ρ j s )}.

[0135] A q value < 0.05 is recorded as a specific ASm6A modification peak in the tumor sample; a q value ≥ 0.05 is recorded as a non-specific ASm6A modification peak in the tumor sample.

[0136] Example 4: Performance Evaluation of the Model for ASm6A Modification Recognition

[0137] I. Stability evaluation of the model used for ASm6A modification recognition

[0138] 1. Experimental Methods

[0139] Simulated data were generated using m6A peaks and loci from real MeRIP-seq sequencing data (GEO: GSM11828594), and were classified based on sequencing read lengths (<75bp, ≥75bp and <100bp, ≥100bp and <150bp, ≥150bp and <300bp, and ≥300bp), library size (<10, ≥10 and <20, ≥20 and <30, ≥30 and <40, ≥40 and <50, ≥50 and <60, ≥60 and <70, ≥70 and <80, ≥80 and <90, and ≥90 million reads), FPKM values ​​of gene expression (0–5, 5–10, and 10–15, all of which are log2(FPKM+1)), the number of SNPs in the m6A peaks (1, 2, 3, 4, and >5), and the number of biological repeats (1, 2, 3, 4, and 6).

[0140] Within each classification category, 50% of the m6A peaks were randomly assigned as allele-specific, i.e., true positive ASm6A (MAF > 0.6), and the remaining m6A peaks were labeled as true negative ASm6A (MAF = 0.5), thus obtaining the datasets for each classification category.

[0141] For each category of dataset, the model shown in Example 2 was used to perform repeated identification 50 times, and the average error rate for identifying each category of dataset was calculated.

[0142] 2. Experimental Results

[0143] The average error rate results when using the model shown in Example 2 to identify datasets for each classification category are as follows: Figure 2As shown, A represents the average error rate of the dataset classified by sequencing read length; B represents the average error rate of the dataset classified by library size; C represents the average error rate of the dataset classified by FPKM value of gene expression; D represents the average error rate of the dataset classified by the number of SNPs in the m6A peak; and E represents the average error rate of the dataset classified by the number of biological repeats.

[0144] The results showed that the average error rate did not change significantly when using the model shown in Example 2 to identify datasets with different sequencing read lengths; however, the average error rate decreased with the increase of library size, FPKM value of gene expression, number of SNPs in m6A peak and number of biological repeats; especially for library size and FPKM value of gene expression, the average error rate decreased significantly.

[0145] However, the error rate remained below 10% when identifying samples for each dataset, indicating that the model shown in Example 2 has excellent stability in identifying ASm6A in samples and is suitable for sequencing data under different conditions.

[0146] II. Evaluation of the model used for ASm6A modification recognition

[0147] 1. Experimental Methods

[0148] Experimental group (M6Allele): Using the datasets of each classification category shown in step one, identification was performed using the model shown in Example 2, and the area under the ROC curve (AUC) of the model shown in Example 2 was calculated based on the true positive and true negative results of ASm6A in the dataset, and the average AUC was calculated.

[0149] Control Group 1 (ASPRIN): Using the datasets of each classification category shown in Step 1, and combining the ASPRIN algorithm developed by ASPRIN et al. (existing technology: https: / / doi.org / 10.1016 / j.ajhg.2019.01.018), identification was performed using default parameter configuration. Based on the identification results, if any SNP within an m6A peak is identified by ASPRIN as having ASm6A modification, then the m6A peak is classified as ASm6A modified; otherwise, it is considered non-ASm6A. The area under the ROC curve (AUC) is calculated based on the true positive and true negative results of ASm6A in the dataset, and the average AUC is calculated.

[0150] Control group 2 (Cao S): Using the datasets of each classification category shown in step one, combined with the algorithm developed by Cao S et al. (existing technology: Genome Res. 2023 Aug; 33(8): 1369-1380. doi: 10.1101 / gr.277704.123. Epub 2023Sep 15), identification was performed using the default parameter configuration. Based on the identification results, if any SNP within an m6A peak is identified by Cao S as having ASm6A modification, then the m6A peak is classified as ASm6A modified; otherwise, it is considered non-ASm6A. The area under the ROC curve (AUC) is calculated based on the true positive and true negative results of ASm6A in the dataset, and the average AUC is calculated.

[0151] 2. Experimental Results

[0152] The mean AUC results for the experimental group, control group 1, and control group 2 are shown in the figure below. Figure 3 As shown in the figure; the area under the ROC curves for the experimental group, control group 1, and control group 2 are shown in the figure. Figure 4 As shown, A is the AUC result of the dataset classified by the number of SNPs in the m6A peak; B is the AUC result of the dataset classified by the FPKM value of gene expression; C is the AUC result of the dataset classified by the number of biological repeats; D is the AUC result of the dataset classified by the sequencing read length; and E is the AUC result of the dataset classified by the library size.

[0153] The results showed that the average AUC value of the model shown in Example 2 was consistently significantly higher than that of the algorithms developed by ASPRIN et al. and Cao S et al.; while the AUC values ​​of the algorithms developed by ASPRIN et al. and Cao S et al. changed significantly with the sequenced conditions during the identification process, and the AUC value increased significantly with the increase of the number of SNPs in m6A minutes, indicating that they could not avoid errors caused by sequencing data noise when identifying data under different sequencing conditions.

[0154] III. Application of the model for ASm6A modification recognition

[0155] 1. Experimental Methods

[0156] The MeRIP-seq dataset (GEO: GSE164151, IP sample: GSM4998285, input sample: GSM4998284) in the GEO database was used as the dataset to be analyzed. The model shown in Example 2 was used for identification, and m6A peaks with significant ASm6A and m6A peaks without significant ASm6A were randomly selected from the identification results and visualized using the IGV tool.

[0157] 2. Experimental Results

[0158] The visualization results of IGV are shown in the figure below. Figure 5 As shown, the results indicate that the model shown in Example 2 can accurately identify the ASm6A peak in the sample and will not incorrectly identify non-ASm6A peaks as ASm6A.

[0159] Example 5: Application of the method for identifying differences in ASm6A between paired samples

[0160] I. Experimental Methods

[0161] Two simulated datasets with the same genotyping were generated using m6A peaks and loci from real MeRIP-seq sequencing data (GEO: GSM11828594), denoted as Sample1 and Sample2, respectively. In each simulated dataset, 50% of the m6A peaks were randomly assigned as allele-specific, i.e., true positive ASm6A (MAF > 0.6), and the remaining m6A peaks were labeled as true negative ASm6A (MAF = 0.5).

[0162] The ASm6A events in Sample1 and Sample2 were identified using Example 2, and the differences in ASm6A between the two datasets (Sample1 and Sample2) were identified according to the method shown in Example 3.

[0163] II. Experimental Results

[0164] The identification results of the differential ASm6A are shown in the figure below. Figure 6 As shown, Figure 6 The results on the left side represent the predicted locations of ASm6A, including those not present in Sample1 and Sample2, present in both Sample1 and Sample2, present only in Sample1, and present only in Sample2. Figure 5 The results on the right side show the actual location of ASm6A. The results show that most of the ASm6A in Sample1 and Sample2 can be accurately identified and identified, while only a small number of ASm6A are incorrectly identified and identified.

[0165] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the scope of protection of the present invention. For those skilled in the art, other variations or modifications can be made based on the above description and ideas, and it is neither necessary nor possible to exhaustively describe all implementation methods here. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the claims of the present invention.

Claims

1. A model for identifying ASm6A modifications, characterized in that, It includes a data acquisition module, a data processing module, a difference significance analysis module, a modification module, and a result output module; The data acquisition module is used to acquire MeRIP-seq sequencing data of the sample to be tested; The data processing module obtains the m6A peak of the sample based on the MeRIP-seq sequencing data of the sample, and obtains the SNP information covered by each m6A peak, and then processes the SNPs covered by each m6A peak. j （m） Calculate the dominance ratio ρ according to formula 10. j （m） And take the logarithm to get the logarithmic advantage ratio y j （m） ; Formula 10: ; Where n j （m） SNPs in the sample to be tested j （m） The number of reads; x （m） ma,j SNP j （m） The number of readings of the main haplotype; μ ASE The method for obtaining the SNP sites is as follows: based on the RNA-seq sequencing data of the sample to be tested, the SNP sites of the sample to be tested are obtained, and the SNPs in the sample to be tested are targeted. j And calculate SNP according to Formula 1. j The main haplotype relative to SNP j Last time, the haplotype dominance ratio ρ j And take the logarithm to get the logarithmic advantage ratio y j ; Formula 1: ; Where n j SNPs in the sample to be tested j The number of reads; x ma,j SNP j Number of readings for the primary haplotype; x 0,j 0.5*n j ; Log dominance ratio y j By combining the processing with a one-way normal random distribution, we obtain the processed y. j For the processed y j Perform sampling and calculate the mean μ of the posterior distribution. ASE ; The log dominance ratio of the SNP sites covered by each m6A peak was calculated. j （m） By combining the processing with a one-way normal random distribution, we obtain the processed y. j （m） ; Regarding the processed y j （m） Sampling is performed and the mean μ of the posterior distribution is calculated. （m） , will μ （m） Taking the antilogarithm gives ρ j （m） And calculate the frequency of the major alleles of each m6A peak; The difference significance analysis module is used to perform difference significance analysis on the frequencies of the major alleles of each m6A peak obtained by the data processing module, and obtain the difference significance analysis results of each m6A peak; The modification identification module shown identifies m6A peaks with significant differences as ASm6A modified peaks based on the analysis results of the differences in each m6A score, and identifies m6A peaks without significant differences as non-ASm6A modified peaks. The result output module is used to output the results obtained by the modification recognition module.

2. The model according to claim 1, characterized in that, In the data processing module, for the processed y j （m） Sampling is performed using the Metropolis-Hastings sampling method, which is an improvement on the Markov chain Monte Carlo algorithm.

3. The model according to claim 1, characterized in that, In the data processing module, the SNP j The main haplotype is SNP. j The base with the highest base reading, the SNP j The last haplotype was SNP. j The bases of the last high-base reading.

4. An ASm6A modification recognition system, comprising an acquisition unit, a storage unit, and a processing unit, characterized in that, The acquisition unit is used to acquire MeRIP-seq sequencing data of the sample to be tested; The storage unit stores program instructions that can be executed by the processing unit; The processing unit comprises the model according to any one of claims 1 to 3; When the program instructions are executed by the processing unit, the MeRIP-seq sequencing data of the sample to be tested obtained by the acquisition unit is input into the model described in any one of claims 1 to 3 to obtain the ASm6A modification recognition result.

5. A computer device, comprising a memory and a processor, characterized in that, The memory stores computer programs that can be executed on the processor; When the computer program is executed by the processor, it implements the model described in any one of claims 1 to 3 and / or the ASm6A modification recognition system described in claim 4.

6. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the model described in any one of claims 1 to 3 and / or the ASm6A modification recognition system described in claim 4.

7. A computer program product, characterized in that, The computer program product includes the computer-readable storage medium of claim 6.

8. A model for identifying differences ASm6A in paired samples, characterized in that, It includes a data acquisition module, an identification module, and a result output module; The data acquisition module is used to acquire MeRIP-seq sequencing data of paired samples and identify them according to the model described in any one of claims 1 to 3, respectively obtaining the ASm6A modification peaks of sample A and sample B; the paired samples are sample A and sample B from different groups but originating from the same individual; The identification module, based on the MeRIP-seq sequencing data of sample A and sample B obtained by the data acquisition module, records the m6A modification peaks where the genomic coordinates of sample A and sample B overlap by more than 50% as m6A modification peaks that exist simultaneously in sample A and sample B. For ASm6A modification peaks that exist in sample A but not in sample B, they are denoted as specific ASm6A modification peaks in sample A. For ASm6A modification peaks that are absent in sample A but present in sample B, they are denoted as nonspecific ASm6A modification peaks of sample A. For ASm6A modification peaks that are present in both sample A and sample B but show different major haplotypes, they are denoted as the specific ASm6A modification peak of sample A. For ASm6A modification peaks that are present in both sample A and sample B but show the same major haplotype, the difference between ASm6A modification peaks in sample A and sample B is determined by hierarchical Bayesian model. ASm6A modification peaks with significant differences are recorded as specific ASm6A modification peaks of sample A, and ASm6A modification peaks with no significant differences are recorded as non-specific ASm6A modification peaks of sample A. The result output module is used to output the results obtained by the identification module.

9. The model according to claim 8, characterized in that, The difference is considered significant when q < 0.05, and the difference is considered insignificant when q ≥ 0.05.

Images