Model for identifying ASm6A modification and application thereof

Abstract

The invention discloses a model for identifying ASm6A modification, and the model is based on MeRIP-seq sequencing data of a to-be-detected sample, and can accurately identify an ASm6A modification peak in the to-be-detected sample. When the model is used for identifying an ASm6A modification peak in a to-be-detected sample, the identification error rate is always lower than 10% no matter what kind of MeRIP-seq sequencing data of the to-be-detected sample is used, and it is indicated that when the model is used for identifying ASm6A modification of the sample, the model has excellent stability and is suitable for identifying the MeRIP-seq sequencing data of the sample under different conditions. Compared with an existing method for recognizing ASm6A modification, the AUC value is obviously increased when the model is used for recognizing ASm6A modification in a sample, and the recognition effect is obviously better.

Description

Technical Field

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

[0002] N6-methyladenosine (m6A) is the most common dynamic and reversible RNA modification in eukaryotic cells. Through the joint action of a series of methylases, demethylases and methylation recognition proteins, m6A modification forms a precise dynamic regulatory cycle in the body, and is widely involved in life activities such as embryonic development, cell apoptosis, sperm development, immunity and stress response. Abnormal m6A modification can cause serious body lesions. For example, the existing technology has found that hypoxic conditions can induce breast cancer stem cell phenotypes 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 intrinsic regulatory mechanism is the key basis for clarifying the pathogenesis of related diseases and developing new therapeutic targets.

[0003] At present, due to the resolution limitation of experimental methods, it is still very difficult to analyze the precise regulatory mechanism of m6A modification at the omics level. Immunoprecipitation combined with next-generation sequencing technology (MeRIP-seq) can simultaneously identify tens of thousands of RNA methylation modification signal peaks with a length of 100 to 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 introduction of the m6A individual-nucleotide resolution cross-linking and immunoprecipitation (miCLIP) method, it is possible to successfully identify m6A modification sites at a single base level in the whole transcriptome at the experimental level. Existing technologies have performed high-throughput identification of m6A modification sites in mammals at a single base accuracy. However, since the technology is still immature, the identification results obtained are still not very stable. It has not been widely used in the field of m6A research. MeRIP-seq is still the most important high-throughput research method. Therefore, improving the quality of the results of identifying m6A modification regions from MeRIP-seq data will be the key basis of current m6A research.

[0005] On the basis of accurately identifying the m6A modified regions of the transcriptome, another important research goal is to analyze the allele-specific modification time. In non-haploid biological clocks, the expression of certain modification sites or certain genes tends to occur on specific chromosomes, which is called allele specificity. Some regulatory events of important life activities, such as gene transcription and epigenetic modification, have significant allele specificity. These allele-specific regulations affect the life activities of cells in many aspects, including imprinting regulation of embryonic development, inactivation of chromosomes, and regulation of gene expression in a specific time and space background. Disordered allele-specific regulatory events can often cause serious body lesions, such as abnormal biological development, cardiovascular and cerebrovascular function lesions, and even cancer. Therefore, the identification and study of allele-specific regulatory 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 time, the study of allele-specific regulatory events was confined to limited genes or chromosome segments. With the development of high-throughput sequencing technology, the study of allele-specific regulatory events has been carried out at the omics level. Sequencing technology can identify single nucleotide polymorphisms (SNPs) on DNA or RNA. Since in non-haploid organisms, chromosomes from the mother and father generally carry different SNP sites, the source of sequencing reads can be inferred through heterozygous SNP sites in the sequence, thereby identifying regulatory signals from different alleles.

[0007] Allele-specific modification is considered to be one of the important reasons for allele-specific expression. Existing research mainly focuses on the DNA level. For example, when studying mouse brain tissue, existing technology found that the DNA methylation levels of more than 1,300 alleles on two homologous chromosomes were significantly different; in addition, existing technology also revealed that more than 10% of human genes are regulated by allele-specific DNA methylation. These allele-specific methylation modifications are highly correlated with SNP sites on the sequence, and there are relatively few studies on allele-specific RNA modifications.

[0008] Compared with the above-mentioned allele-specific regulation time, there are still no reports on the allele-specific research of m6A modification on the transcriptome. The m6A modification controls the correct functions of many genes related to metabolism and growth differentiation in the body. The existing technology has found that knocking out METTL3 in mouse embryonic stem cells can prevent their differentiation; in fruit flies, the knockout of METTL3's homologous protein Ime4 also has a sublethal effect, which will cause damage to the NOTCH signaling pathway and affect the fertility of fruit flies; in addition, the Ime4 gene also plays an important role in yeast meiosis. This shows that m6A modification is crucial to the early development of biological sex organs and early embryonic development. Therefore, studying whether m6A modification is specific between alleles is an important entry point for deciphering the mechanism of embryonic development.

[0009] The analysis of m6A modification signals from MeRIP-seq sequencing data mainly relies on computational biology methods. Researchers first used Fisher's Exact Test method to construct the first peak identification algorithm for MeRIP-seq. This method 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. Its 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 Hidden Markov Model (Hidden Markov Model). However, these algorithms must first assume that the reads of MeRIP-seq conform to the Beta binomial distribution, so they will introduce obvious errors for samples with low sequencing quality, and the parameter optimization of the HMM model needs to be based on a certain number of repeated samples to obtain more accurate results, but general MeRIP-seq sequencing will only perform a small number of repeated experiments, which also limits the identification performance of the algorithm to a certain extent; in addition, researchers will also use MACS to identify the m6A signal of MeRIP-seq sequencing data, but MACS is a peak identification software designed specifically for ChIP-seq analysis. Applying it to MeRIP-seq sequencing analysis with different principles will produce obvious errors, and there will also be coordinate offsets.

[0010] Therefore, there is an urgent need for a method that can accurately identify mA peak signals based on MeRIP-seq sequencing, thereby achieving accurate identification of allele-specific mA. Summary of the invention

[0011] The purpose of the present invention is to overcome the above-mentioned deficiencies of the prior art and provide a model for identifying ASm6A modification and its application.

[0012] A first object of the present invention is to provide a model for identifying ASm6A modifications.

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

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

[0015] A fourth object of the present invention is to provide a computer device.

[0016] A fifth object of the present invention is to provide a computer-readable storage medium.

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

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

[0019] In order to achieve the above object, the present invention is implemented by the following scheme:

[0020] A model for identifying ASm6A modifications, 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 obtain MeRIP-seq sequencing data of the sample to be tested;

[0022] The data processing module obtains the m6A peak of the sample to be tested based on the MeRIP-seq sequencing data of the sample to be tested, obtains the SNP site information covered by each m6A peak, and j (m) , calculate the odds ratio ρ according to formula 10 j (m) , and take the logarithm to get the log odds ratio y j (m) ;

[0023] Formula 10:

[0024] where n j (m) SNP in the sample to be tested j (m) The number of read counts; x (m) ma,j SNP j (m)The number of reads of the major haplotype on ;

[0025] μ ASE The acquisition method is as follows: based on the RNA-seq sequencing data of the sample to be tested, the SNP site of the sample to be tested is obtained, and the SNP in the sample to be tested is j , and calculate SNP according to formula I j The odds ratio of the major haplotype to the minor haplotype is ρ j , and take the logarithm to get the log odds ratio y j ;

[0026] Formula 1:

[0027] where n j SNP in the sample to be tested j The number of read counts; x ma,j SNP j The number of reads of the main haplotype; x 0,j =0.5*n j ;

[0028] The log odds ratio y j Combined with the one-way normal random distribution, the processed y is obtained j , for the processed y j Sampling and calculating the mean μ of the posterior distribution ASE ;

[0029] The logarithmic odds ratio y of the SNP sites covered by each m6A peak j (m) Combined with the one-way normal random distribution, the processed y is obtained j (m) ; For the processed y j (m) , perform sampling and calculate the mean μ of the posterior distribution (m) , μ (m) Taking the antilog gives ρ j (m) ', and calculate the frequency of the major allele of each m6A peak;

[0030] The difference significance analysis module is used to perform a difference significance analysis on the frequency of the major allele of each m6A peak obtained by the data processing module to obtain a difference significance analysis result of each m6A peak;

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

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

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

[0034] Preferably, in the data processing module, for the processed y j (m) , sampling is performed using the Metropolis-Hastings sampling method improved based 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 read, the SNP j Last minor haplotype is SNP j The base of the last high-base read.

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

[0037] The present invention also claims the use of any of the above-mentioned models in the preparation of products for identifying ASm6A.

[0038] The present invention also claims protection for an ASm6A modification identification 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 above-mentioned models to obtain the ASm6A modification recognition result.

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

[0043] When the computer program is executed by a processor, any of the above-mentioned models and / or the above-mentioned ASm6A modification identification system are implemented.

[0044] The present invention also seeks to protect a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements any of the above-mentioned models and / or the above-mentioned ASm6A modification identification system.

[0045] The present invention also seeks to protect a computer program product, wherein the computer program product comprises the above-mentioned computer-readable storage medium.

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

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

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

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

[0050] For the ASm6A modification peaks present in sample A but not in sample B, they were recorded as the specific ASm6A modification peaks in sample A;

[0051] For the ASm6A modification peak that is absent in sample A but present in sample B, it is recorded as the nonspecific ASm6A modification peak of sample A;

[0052] For ASm6A modification peaks that exist in both sample A and sample B but show different major haplotypes, they are recorded as specific ASm6A modification peaks of sample A;

[0053] For ASm6A modification peaks that exist in both sample A and sample B but show the same major haplotype, the hierarchical Bayesian model is used to judge the difference between the ASm6A modification peaks in sample A and sample B. The ASm6A modification peaks with significant differences are recorded as the specific ASm6A modification peaks of sample A, and the ASm6A modification peaks with no significant differences are recorded as the nonspecific ASm6A modification peaks of sample A.

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

[0055] Preferably, the significant difference is when the q value is less than 0.05, and the insignificant difference is when the q value is greater than or equal to 0.05.

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

[0057] The present invention provides a model for identifying ASm6A modification, and the model is based on the MeRIP-seq sequencing data of the sample to be tested, and can accurately identify the ASm6A modification peak in the sample to be tested. When the model of the present invention identifies the ASm6A modification peak in the sample to be tested, no matter what MeRIP-seq sequencing data of the sample to be tested is used, the error rate of recognition is always less than 10%, indicating that the model of the present invention has excellent stability when identifying the ASm6A modification of the sample, and is suitable for identifying the MeRIP-seq sequencing data of the sample under different conditions. And compared with the existing method for identifying ASm6A modification, the AUC value when using the model of the present invention to identify the ASm6A modification in the sample is significantly improved, and the recognition effect is significantly better. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 This is a schematic diagram of the principle of a model for ASm6A modification identification shown in Example 2;

[0059] Figure 2 The figure is a graph of the average error rate results when the data sets of each classification category are identified according to the method shown in Example 3; A is the average error rate result of the data set classified according to the sequencing read length; B is the average error rate result of the data set classified according to the library size; C is the average error rate result of the data set classified according to the FPKM value of gene expression; D is the average error rate result of the data set classified according to the number of SNPs in the m6A peak; E is the average error rate result of the data set classified according to the number of biological replicates;

[0060] Figure 3 The average AUC result graph of the experimental group, control group 1 and control group 2 in Example 4;

[0061] Figure 4The results are the area under the ROC curve of the experimental group, control group 1 and control group 2 in Example 4; A is the AUC result diagram when the data set classified by the number of SNPs in the m6A peak is identified; B is the AUC result diagram when the data set classified by the FPKM value of gene expression is identified; C is the AUC result diagram when the data set classified by the number of biological replicates is identified; D is the AUC result diagram when the data set classified by the sequencing read length is identified; E is the AUC result diagram when the data set classified by the library size is identified;

[0062] Figure 5 This is the IGV visualization result diagram in Example 4;

[0063] Figure 6 This is a diagram showing the identification results of the differential ASm6A in Example 5. DETAILED DESCRIPTION

[0064] The present invention is further described in detail below in conjunction with the accompanying drawings and specific examples of the specification. The examples are only used to explain the present invention and are not used to limit the scope of the present invention. The test methods used in the following examples are conventional methods unless otherwise specified; the materials and reagents used are reagents and materials that can be obtained from commercial channels unless otherwise specified.

[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 difference significance analysis module, a prediction module and a result output module.

[0067] The data acquisition module is used to obtain 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] Targeting SNP j , calculate SNP according to formula 1 j The odds ratio of the major haplotype to the minor haplotype is ρ j, and take the logarithm of the odds ratio according to formula 2 to obtain the logarithmic odds ratio y j ; Among them, SNP j is the SNP site in the sample to be tested;

[0071] Formula 1:

[0072] Formula 2:

[0073] where n j SNP in the sample to be tested j The number of read counts; x ma,j SNP j The number of reads of the main haplotype on x 0,j =0.5*n j ;

[0074] The log odds ratio y obtained from formula 2 j Combined with the single normal random effects model (REM), the processed data are processed according to formula 3 to formula 6 to obtain the processed y 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 ratio of the sample to be tested, σ j 2 is 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 based on the Markov chain Monte Carlo (MCMC) algorithm is used to perform sampling according to Formula 7, and the marginal probability is calculated according to Formula 8, and then the posterior distribution of μ is calculated according to Formula 9, and the mean of the posterior distribution of μ is calculated (μ 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 ) were analyzed for significant differences, as follows:

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

[0086] In order to distinguish significant allele-specific events, the major allele frequency (MAF) distribution needs to be obtained under the null hypothesis. Due to the lack of real MeRIP-seq data that meet the standards, sequencing data needs to be simulated to obtain the read counts of the major and minor alleles of SNPs without significant ASE or ASm6A events; the read count distribution of single SNPs captured by sequencing is fitted by using negative binomial distribution (NBD).

[0087] Assume that each SNP site is covered (SNP i ) is shown in Formula 10;

[0088] Formula 10:

[0089] Where μ is the SNP when there is no allele-specific event i The theoretical number of read segments is equal to 0.5N i , N i SNP i The number of reads count; k is the discreteness parameter.

[0090] The dispersion parameter k is obtained by obtaining whole genome sequencing (WGS) data from the 1000 Genomes database (https: / / www.internationalgenome.org / ) and counting the number of coverage reads N for each SNP site in each sample. i and the number of reads of the alleles, and integrate N by the maximum likelihood method i and x, and obtain the dispersion parameter k.

[0091] Then, the FPKM values ​​of all genes were collected from the Cancer Genome Atlas Project (TCGA) and fitted using the fitter package in Python. The genes were divided into six categories of genes according to the distribution type of FPKM (as described in the prior art 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.), and the overall distribution of FPKM for each category of genes was refitted, and the simulated FPKM of each gene was obtained by sampling based on the distribution. The total number of reads N for each SNP on the gene was calculated based on the simulated FPKM and gene length of each gene. i '.

[0092] Based on the discreteness parameter k and the N obtained by simulation i ', combined with formula 10 to obtain the SNP of each meta-unit (gene) i The negative binomial distribution of the number of allele reads is used, and the number of reads of the major and minor alleles is simulated by sampling, and the odds ratio ρ is calculated according to Formula 1 i ', and calculate the MAF on each gene; MAF = {ρ i ' / (1+ρ i ')}.

[0093] The generalized Pareto distribution (GPD) was introduced to estimate the tail distribution of MAF' of different types of genes, and the statistical significance threshold was calculated. The 94% quantile was used as the initial threshold t. For MAF greater than the initial threshold t, the scale parameter and shape parameter of GPD in formula 11 were combined to estimate its distribution, and the tail distribution of MAF values ​​in the null hypothesis data was obtained.

[0094] Formula 11:

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

[0096] (2) Analysis of significance of differences

[0097] For the mean of the μ posterior distribution obtained by the data processing module (μ ASE ), taking the antilogarithm of it, we get ρ ASE , and calculate μ ASE Corresponding MAF value (R); MAF = {ρ ASE / (1+ρ ASE )}.

[0098] μ ASE The corresponding MAF value (R) is combined with formula 12 to perform a significant difference analysis to obtain the significant difference analysis result;

[0099] Formula 12:

[0100] Among them, t is 94%, R simj The MAF value representing the null hypothesis that ASm6A / ASE does not exist in the data, N sim Indicates R in the simulated data simj total.

[0101] The P value was corrected by the BH method to obtain the q value, and a q value less than 0.05 was considered to be significantly different.

[0102] The prediction module obtains the difference significance analysis result based on the difference significance analysis module, and determines whether the gene is an ASE gene, specifically: μ ASE The genes corresponding to the corresponding MAF values ​​were judged as ASE genes; the μ ASE The genes corresponding to the corresponding MAF values ​​were judged as non-ASE genes.

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

[0104] Example 2 A model for ASm6A modification identification

[0105] A schematic diagram of a model for ASm6A modification recognition is shown below. Figure 1 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 obtain the 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, specifically including:

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

[0110] Targeting SNPs covered on m6A peaks j (m) , calculate the odds ratio ρ according to formula 10 j (m) , and take the logarithm of the odds ratio according to formula 11 to obtain the logarithmic odds ratio y j (m) ;

[0111] Formula 10:

[0112] Formula 11:

[0113] where n j (m) SNP in the sample to be tested j (m) The number of read counts; x (m) ma,j SNP j (m) The number of reads of the major haplotype on ASE SNPs were processed by the functional model for identifying genes with significant allele-specific expression as described in Example 1 i (m) Get;

[0114] The log odds ratio y obtained from formula 11 j (m) Combined with the single normal random distribution (REM), it is processed according to formula 12 and formula 15 to 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 odds ratio, σ j 2 is the variance of the log odds ratio; s 2 =10.

[0120] For the processed y j (m) , the Metropolis-Hastings (MH) sampling method based on the Markov chain Monte Carlo (MCMC) algorithm is used for sampling, and the mean μ of the posterior distribution is calculated (m) , μ (m) Taking the antilog gives ρ j (m) '(i.e. j (m) The expected value of m6A peak was calculated and the major allele frequency (MAF) of m6A peak 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 described in Example 1 is that: the difference significance analysis is performed on the frequency (MAF) of the major allele on the m6A peak according to (2) difference significance analysis in the difference significance analysis module described in Example 1 to obtain the difference significance analysis result (q value) of the m6A peak.

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

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

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

[0126] Paired samples are samples from different groups but derived 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] The m6A peaks whose genomic coordinates overlap more than 50% between the tumor sample and the normal sample are regarded as coming from the same m6A peak (i.e., the m6A peaks that exist in both the tumor sample and the normal sample). The tumor sample and the normal sample in the paired samples are processed using the model shown in Example 2, respectively, to obtain the ASm6A modified peaks of the tumor sample and the normal sample, and classify them:

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

[0129] (2) The ASm6A peaks that exist in normal samples but not in tumor samples are recorded as nonspecific ASm6A modification peaks of tumor samples.

[0130] (3) For ASm6A peaks that exist in both tumor samples and normal samples but show different major haplotypes, they are recorded as specific ASm6A modification peaks of tumor samples.

[0131] (4) For ASm6A peaks that exist in both tumor samples and normal samples and show the same major haplotype, the hierarchical Bayesian model is used to evaluate the difference. The consistent heterozygous SNP sites in the overlapping region of the ASm6A peak are regarded as available sites for downstream analysis to ensure consistent haplotype construction between the two samples. For each SNP site, the odds ratio ρ is 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 of the locus, n tumor，j SNPs in tumor IP samples j The number of reads count at the site, μ tumor，b is the expected value of the log-odds ratio of the major allele of allele-specific expression calculated for the tumor input sample (i.e., general RNA-seq), where ρ under the null hypothesis tumor，j =ρ normal，j ;

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

[0135] A q value of <0.05 was recorded as a specific ASm6A modification peak of the tumor sample; a q value of ≥0.05 was recorded as a nonspecific ASm6A modification peak of the tumor sample.

[0136] Example 4 Performance evaluation of the model for ASm6A modification identification

[0137] 1. Stability assessment of the model used for ASm6A modification recognition

[0138] 1. Experimental methods

[0139] Simulated data were generated using m6A peaks and sites from real MeRIP-seq sequencing data (GEO: GSM11828594) and classified based on sequencing read length (<75 bp, ≥75 bp and <100 bp, ≥100 bp and <150 bp, ≥150 bp and <300 bp, and ≥300 bp), 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 value of gene expression (0-5, 5-10, and 10-15, all log2(FPKM+1)), number of SNPs in m6A peaks (1, 2, 3, 4, and >5), and number of biological replicates (1, 2, 3, 4, and 6).

[0140] In each classification category, 50% of the mA peaks were randomly assigned as allele-specific, i.e., true positives of ASm6A (MAF>0.6), and the remaining mA peaks were marked as true negatives of ASm6A (MAF=0.5), obtaining data sets for each classification category.

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

[0142] 2. Experimental results

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

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

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

[0146] Evaluation of models for ASm6A modification recognition

[0147] 1. Experimental methods

[0148] Experimental group (M6Allele): Using the data sets of each classification category shown in step 1, 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 data set, and the average AUC was calculated.

[0149] Control group 1 (ASPRIN): Using the datasets of each classification category shown in step 1, combined with the algorithm ASPRIN developed by ASPRIN et al. (prior art: https: / / doi.org / 10.1016 / j.ajhg.2019.01.018), the default parameter configuration was used for identification. For its identification results, if any SNP within an mA peak was identified by ASPRIN as having ASm6A modification, the mA peak was classified as ASm6A modification, otherwise it was regarded as non-ASm6A, and the area under the ROC curve (AUC) was calculated based on the true positive and true negative results of ASm6A in the dataset, and the average AUC was calculated.

[0150] Control group 2 (Cao S): Using the datasets of each classification category shown in step 1, combined with the algorithm developed by Cao S et al. (Prior art: Genome Res. 2023Aug; 33(8): 1369-1380.doi:10.1101 / gr.277704.123.Epub 2023Sep 15), identification was performed using the default parameter configuration. For its identification results, if any SNP within an mA peak is identified by Cao S as having ASm6A modification, the mA peak is classified as ASm6A modification, otherwise it is considered non-ASm6A, and 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 average AUC results of the experimental group, control group 1 and control group 2 are shown in the figure Figure 3 As shown; the area under the ROC curve results of the experimental group, control group 1 and control group 2 are shown in Figure 4 As shown, A is the AUC result graph when identifying a dataset classified according to the number of SNPs in the m6A peak; B is the AUC result graph when identifying a dataset classified according to the FPKM value of gene expression; C is the AUC result graph when identifying a dataset classified according to the number of biological replicates; D is the AUC result graph when identifying a dataset classified according to sequencing read length; E is the AUC result graph when identifying a dataset classified according to library size.

[0153] The results showed that the average AUC value when using the model shown in Example 2 for identification was always significantly higher than that of the algorithm developed by ASPRIN et al. and the algorithm developed by Cao S et al.; during the identification process, the AUC values ​​of the algorithms developed by ASPRIN et al. and Cao S et al. changed significantly with the changes in sequencing conditions, and the AUC values ​​increased significantly with the increase in the number of SNPs in the m6A minute, indicating that it is difficult to avoid errors caused by sequencing data noise when identifying data under different sequencing conditions.

[0154] 3. Application of the model for ASm6A modification identification

[0155] 1. Experimental methods

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

[0157] 2. Experimental results

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

[0159] Example 5 Application of the method for identifying differential ASm6A in paired samples

[0160] 1. Experimental Methods

[0161] The m6A peaks and sites of real MeRIP-seq sequencing data (GEO: GSM11828594) were used to generate two simulated data with the same genotyping, which were recorded as Sample 1 and Sample 2. In each simulated data, 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 marked as true negative ASm6A (MAF=0.5).

[0162] Example 2 was used to identify the ASm6A events in Sample 1 and Sample 2, respectively, and the differential ASm6A between the two data sets (in Sample 1 and Sample 2) was identified according to the method shown in Example 3.

[0163] 2. Experimental Results

[0164] The identification results of differential ASm6A are shown in the figure Figure 6 As shown, Figure 6 The results on the left in the middle are the predicted positions of ASm6A, including those that do not exist in both Sample 1 and Sample 2, those that exist in both Sample 1 and Sample 2, those that exist only in Sample 1, and those that exist only in Sample 2. Figure 5 The result on the middle right is the actual position of ASm6A. The results show that most of the ASm6A in Sample 1 and Sample 2 can be accurately identified and identified, and only a small part of ASm6A is 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 rather than to limit the protection scope of the present invention. For ordinary technicians in this field, other different forms of changes or modifications can be made based on the above descriptions and ideas. It is not necessary and impossible to list all the implementation methods here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope 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 obtain MeRIP-seq sequencing data of the sample to be tested; The data processing module obtains the m6A peak of the sample to be tested based on the MeRIP-seq sequencing data of the sample to be tested, obtains the SNP site information covered by each m6A peak, and j (m) , calculate the odds ratio ρ according to formula 10 j (m) , and take the logarithm to get the log odds ratio y j (m) ; where n j (m) SNP in the sample to be tested j (m) The number of read counts; x (m) ma,j SNP j (m) The number of reads of the major haplotype on ; μ ASE The acquisition method is as follows: based on the RNA-seq sequencing data of the sample to be tested, the SNP site of the sample to be tested is obtained, and the SNP in the sample to be tested is j , and calculate SNP according to formula I j The odds ratio of the major haplotype to the minor haplotype is ρ j , and take the logarithm to get the log odds ratio y j ; Formula 1: where n j SNP in the sample to be tested j The number of read counts; x ma,j SNP j The number of reads of the main haplotype; x 0,j =0.5*n j ; The log odds ratio y j Combined with the one-way normal random distribution, the processed y is obtained j , for the processed y j Sampling and calculating the mean μ of the posterior distribution ASE ; The logarithmic odds ratio y of the SNP sites covered by each m6A peak j (m) Combined with the one-way normal random distribution, the processed y is obtained j (m) ; For the processed y j (m) , perform sampling and calculate the mean μ of the posterior distribution (m) , μ (m) Taking the antilog gives ρ j (m) ', and calculate the frequency of the major allele of each m6A peak; The difference significance analysis module is used to perform a difference significance analysis on the frequency of the major allele of each m6A peak obtained by the data processing module to obtain a difference significance analysis result of each m6A peak; The modification recognition module shown identifies the m6A peaks with significant differences as ASm6A modification peaks based on the significant difference analysis results of each m6A score, and identifies the m6A peaks without significant differences as non-ASm6A modification peaks; The result output module is used to output the result obtained by the modification identification 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 improved based 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 read, the SNP j Last minor haplotype is SNP j The base of the last high-base read.

4. Use of the model described in any one of claims 1 to 3 in the preparation of a product for identifying ASm6A.

5. An ASm6A modification identification 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 described in 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.

6. A computer device comprising a memory and a processor, characterized in that: The memory stores a computer program executable on the processor; When the computer program is executed by a processor, the model described in any one of claims 1 to 3 and / or the ASm6A modification identification system described in claim 5 are implemented.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the model described in any one of claims 1 to 3 and / or the ASm6A modification identification system described in claim 5 are implemented.

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

9. 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; The paired samples are sample A and sample B from different groups and derived from the same individual; 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 to obtain ASm6A modification peaks of sample A and sample B respectively; The identification module records the m6A modification peaks whose genomic coordinates overlap by more than 50% between sample A and sample B as m6A modification peaks that exist simultaneously in sample A and sample B based on the MeRIP-seq sequencing data of sample A and sample B obtained by the data acquisition module; For the ASm6A modification peaks present in sample A but not in sample B, they were recorded as the specific ASm6A modification peaks in sample A; For the ASm6A modification peak that is absent in sample A but present in sample B, it is recorded as the nonspecific ASm6A modification peak of sample A; For ASm6A modification peaks that exist in both sample A and sample B but show different major haplotypes, they are recorded as specific ASm6A modification peaks of sample A; For ASm6A modification peaks that exist in both sample A and sample B but show the same major haplotype, the hierarchical Bayesian model is used to judge the difference between the ASm6A modification peaks in sample A and sample B. The ASm6A modification peaks with significant differences are recorded as the specific ASm6A modification peaks of sample A, and the ASm6A modification peaks with no significant differences are recorded as the nonspecific ASm6A modification peaks of sample A. The result output module is used to output the result obtained by the identification module.

10. The model according to claim 9, characterized in that The significant difference is when the q value is less than 0.05, and the insignificant difference is when the q value is greater than or equal to 0.05.

Images