Ac4C modification site identification method

Abstract

The invention relates to the technical field of biology, and discloses an ac4C modification site identification method which comprises the following steps: obtaining a processing group and a control group based on ac4C-seq high-throughput sequencing data, and identifying ac4C modification sites under single base resolution; based on the ac4C modification site data, performing difference analysis between the processing group and the control group to obtain different ac4C modification sites; processing and analyzing the ac4C modification site data, and drawing distribution characteristics of modification sites; and annotating and analyzing the data of the different ac4C modification sites to obtain biological process information related to ac4C modification. According to the method, the ac4C modification sites are recognized under the single base resolution, the ac4C modification sites with inter-group differences are obtained, the accuracy and efficiency of ac4C modification site identification are improved, the generalization ability of the ac4C modification sites is enhanced, and meanwhile distribution and biological significance of the modification sites can be analyzed.

Description

Technical Field

[0001] The present invention relates to the field of biotechnology, and in particular to a method for identifying ac4C modification sites. Background Art

[0002] More than 170 different types of post-transcriptional modifications have been found on RNA, among which the main modifications on mRNA, such as N6-adenylate methylation (N6-Methyladenosine, m6A), N1-adenylate methylation (N1-Methyladenosine, m1A) and 5-methylcytosine (5-Methylcytosine, m5C), have been widely studied. In addition to these mainstream modifications, recent studies have found that N4-acetylcytidine (N4-acetylcytidine, ac4C) modification also exists on mammalian mRNA, which plays a certain role in enhancing the stability of mRNA and improving translation efficiency. However, there is little research on ac4C acetylation modification, and the research methods are not perfect.

[0003] At present, there are two main methods for detecting ac4C modification sites, namely acRIP-seq based on antibody enrichment and ac4C-seq based on chemical reduction. Among them, acRIP-seq uses ac4C-specific antibodies for co-immunoprecipitation to enrich and recover RNA fragments containing ac4C modifications. This method is relatively mature, but the co-immunoprecipitation method usually results in a resolution of 100-150nt, which cannot accurately locate the ac4C acetylation modification site. ac4C-seq reacts sodium cyanoborohydride (NaCNBH3) with ac4C sites under acidic conditions to generate reduced N4-acetyltetrahydrocytidine. This change in base structure causes base mismatches in the sequence during reverse transcription, and the conversion from C to T is detected during high-throughput sequencing of the cDNA sequence, thereby detecting the ac4C modification site at single-base resolution.

[0004] In order to identify ac4C modification sites, researchers usually develop some specialized tools. The latest tools for identifying ac4C sites are PACES and XG-ac4C. PACES is a method based on RNA sequence features and random forest machine learning, which combines position-specific dinucleotide sequence features and K nucleotide frequencies to extract sequence features. The XG-ac4C method is a computational model based on the extreme gradient boosting method to identify ac4C modification sites, which uses nucleotide chemical properties, nucleotide density, Kmer, single encoding, electron-ion interaction pseudopotential, and trinucleotide electron-ion interaction pseudopotential to represent the mRNA sequences in the dataset.

[0005] However, both methods have certain limitations. First, the models constructed by these two methods are not general enough. For example, they only have good prediction results for human sequencing data and are not applicable to other species. Second, the resolution of the prediction results is not enough. Although ac4C-seq can achieve single-base resolution at the experimental method level, the resolution of the prediction results of these two methods can only reach a dozen bases.

[0006] There is currently no reliable method to accurately identify ac4C modification sites at the single-base level using ac4C-seq high-throughput sequencing data. In addition, there is a lack of statistical methods to identify differential ac4C modification sites between different groups, as well as a lack of tools to link ac4C modification with possible biological processes to further explore its potential biological significance. Therefore, the development of new analytical tools and statistical methods to identify and understand the biological role of ac4C modification is an urgent problem to be solved. Summary of the invention

[0007] The present invention provides an ac4C modification site identification method to achieve the recognition of ac4C modification sites at single-base resolution, obtain inter-group differential ac4C modification sites, propose an analysis method to identify the distribution characteristics of ac4C modification sites and determine the potential connection between ac4C modification and biological processes, thereby improving the accuracy and efficiency of ac4C modification site identification and enhancing its generalization ability, and at the same time being able to analyze the distribution and biological significance of the modification sites.

[0008] The present invention provides a method for identifying ac4C modification sites, comprising:

[0009] Based on the data from ac4C-seq high-throughput sequencing, the treated and controlled groups were obtained, and the ac4C modification sites were identified at single-base resolution.

[0010] Based on the ac4C modification site data, performing difference analysis between the treatment group and the control group to obtain differential ac4C modification sites;

[0011] The ac4C modification site data are processed and analyzed to draw the distribution characteristics of the modification sites; the differential ac4C modification site data are annotated and analyzed to obtain biological process information related to ac4C modification.

[0012] Furthermore, the step of obtaining the treatment group and the control group based on the data from the ac4C-seq high-throughput sequencing and identifying the ac4C modification site at a single-base resolution comprises:

[0013] The data from the ac4C-seq high-throughput sequencing were quality controlled, the set sequences were removed and the PCR bias was corrected. The quality-controlled data were compared with the reference genome of the corresponding species using STAR software to obtain the bam file of the comparison results of each sample after and without chemical reduction treatment; the comparison results after chemical reduction treatment were used as the treatment group, and the comparison results without chemical treatment were used as the control group;

[0014] The pileup mode of JACUSA software was used to count the frequency of occurrence of four bases, A, C, G, and T, at each site in the bam file;

[0015] According to the frequency of bases at each site, the base error rate of the site is calculated;

[0016] For each site, statistical tests were performed based on the number of correct bases and incorrect bases in the Treat sample and the Control sample;

[0017] The tested data were filtered to obtain statistically significant potential ac4C modification sites, and the filtered data were written into the output file to complete the identification of ac4C modification sites.

[0018] Furthermore, in the step of calculating the base error rate of each site according to the frequency of the base at the site, the calculation formula is:

[0019] Correct base number: the frequency of cytosine occurrence in the site;

[0020] Number of incorrect bases: the sum of the frequencies of the four bases at the site minus the number of correct bases;

[0021] Base error rate = number of incorrect bases / sum of the frequencies of occurrence of the four bases.

[0022] Furthermore, in the step of performing statistical testing on each site based on the number of correct bases and the number of incorrect bases in the Treat sample and the Control sample, the statistical testing method is:

[0023] If each sample had multiple replicates, the generalized linear mixed effects model GLMM was used for statistical testing;

[0024] If there was only one replicate for each sample, the chi-square test was used for statistical testing.

[0025] Furthermore, the step of performing difference analysis between the treatment group and the control group based on the ac4C modification site data to obtain differential ac4C modification sites comprises:

[0026] Merging the data of the treatment group and the control group, and filling all missing values ​​in the merged data;

[0027] Normalize the padded data according to the scaling factor;

[0028] For the standardized data, different methods are used to perform difference tests on the treatment group and the control group according to whether there are repeated samples;

[0029] The tested data were filtered to obtain potential differential ac4C modification sites with statistical significance, and the filtered data were written into the output file to complete the identification of differential ac4C modification sites.

[0030] Furthermore, in the step of performing a difference test on the treatment group and the control group using different methods according to whether there are repeated samples, the test method is:

[0031] If each sample had multiple replicates, the generalized linear mixed effects model GLMM was used for difference testing;

[0032] If there was only one replicate for each sample, Fisher's exact test was used for difference testing.

[0033] Furthermore, the generalized linear mixed effects model GLMM simultaneously handles fixed effects and random effects;

[0034] The fixed effect part refers to cytosine information and grouping information, the former is the dependent variable and the latter is the independent variable; wherein the grouping includes Treat and Control;

[0035] The random effect part refers to the sample information, and each sample has a random intercept.

[0036] Furthermore, the steps of processing and analyzing the ac4C modification site data to draw the distribution characteristics of the modification sites; and annotating and analyzing the differential ac4C modification site data to obtain biological process information related to ac4C modification include:

[0037] Analyze the distribution characteristics of ac4C modification sites on chromosomes and gene elements;

[0038] Gene annotation was performed based on the differential information of ac4C modification sites to obtain genes related to ac4C modification;

[0039] The genes are used to perform enrichment analysis to identify biological processes related to ac4C modification; wherein the enrichment analysis includes GO enrichment analysis and KEGG enrichment analysis.

[0040] Furthermore, the step of analyzing the distribution characteristics of ac4C modification sites on chromosomes and gene elements includes:

[0041] According to the length information of the chromosome, it is divided into continuous intervals of the same length. Finally, the number of sites modified by ac4C in each interval is counted, and the statistical results are drawn into a distribution characteristic diagram to analyze the distribution characteristics of ac4C modified sites on the chromosome.

[0042] The genome annotation gff file of the species was parsed to obtain the basic information of each element, and then the distribution of ac4C modification sites on the elements was counted and mapped to analyze the distribution characteristics of ac4C modification sites on gene elements; wherein the elements include UTR5, CDS, and UTR3.

[0043] The beneficial effects of the present invention are:

[0044] The present invention provides a simple method for identifying ac4C modification sites for chemical ac4C-seq data, improves the resolution of ac4C modification site identification to single base, and provides a method for identifying ac4C modification sites with differences between groups. For experiments with repeated samples, a generalized linear mixed effects model (GLMM) is used to test the data, taking into account both the fixed effects of bases and groupings and the random effects between different samples. In addition, the distribution characteristics of the identified ac4C modification sites are analyzed and the role of ac4C modification in cell function and metabolic pathways is discovered through enrichment analysis. Finally, the method for identifying ac4C modification sites is written in Julia language, with fast calculation speed and high efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 Schematic diagram of the process of identifying ac4C modification sites of the present invention.

[0046] Figure 2 Schematic diagram of the specific process of step T1 in a specific example of the present invention.

[0047] Figure 3 Schematic diagram of the specific process of step T2 in a specific example of the present invention.

[0048] Figure 4 Schematic diagram of the specific process of step T3 in a specific example of the present invention.

[0049] Figure 5 Schematic diagram of the specific process of step T4 in a specific example of the present invention.

[0050] Figure 6 Schematic diagram of the specific process of step T5 in a specific example of the present invention.

[0051] Figure 7 Schematic diagram of the specific process of step T6 in a specific example of the present invention.

[0052] Figure 8Schematic diagram of the specific process of step T7 in the specific example of the present invention.

[0053] Fig. 9 It is a schematic diagram of the exon region location information file in a specific example of the present invention.

[0054] Fig.10 A schematic diagram is added for calculating variation information in a specific example of the present invention.

[0055] The realization of the purpose, functional features and advantages of the present invention will be further explained in conjunction with embodiments and with reference to the accompanying drawings. DETAILED DESCRIPTION

[0056] It should be understood that the specific embodiments described herein are only used to explain the present invention, and are not used to limit the present invention.

[0057] In the research field of ac4C modification sites, there are currently several major challenges. First, existing tools require a large amount of data for training when identifying ac4C modification sites, which not only increases the learning cost, but also prolongs the analysis time, and the generalization ability of these tools is generally poor. Secondly, the analysis results of existing tools can usually only reach a resolution of more than a dozen bases, and there is a lack of methods that can accurately identify ac4C modification sites at single-base resolution. In addition, considering the random effects between samples, there is currently a lack of software that can effectively count and identify ac4C modification sites with differences between groups. Finally, for the identified ac4C modification sites, there is currently a lack of effective tools to analyze their distribution characteristics and explore the potential connection between these modification sites and biological processes. Therefore, the development of new tools and methods to improve the accuracy and efficiency of ac4C modification site identification, enhance its generalization ability, and analyze the distribution and biological significance of modification sites at the same time is a key problem that the present invention aims to solve.

[0058] To fill the current research gap in the field, this paper proposes a method for identifying ac4C modification sites based on chemical ac4C-seq data to achieve the identification of ac4C modification sites at single-base resolution and obtain differential ac4C modification sites between groups. At the same time, an analysis method is proposed to identify the distribution characteristics of ac4C modification sites and explore the potential connection between ac4C modification and some biological processes.

[0059] like Figure 1 As shown, the ac4C modification site identification method provided by the present invention comprises:

[0060] S1. Based on the data from ac4C-seq high-throughput sequencing, the treatment group and the control group were obtained, and the ac4C modification sites were identified at single-base resolution. The specific steps are as follows:

[0061] S11. Perform quality control on the data from ac4C-seq high-throughput sequencing, remove the set sequences (connectors, low-quality sequences) and correct the PCR bias. Use STAR software to align the data after quality control with the reference genome of the corresponding species, and obtain the bam file of the comparison results of each sample treated with chemical reduction (hereinafter referred to as the treatment group) and without chemical treatment (hereinafter referred to as the control group) for subsequent analysis.

[0062] S12. Use the pileup mode of JACUSA software to count the frequency of occurrence of the four bases A, C, G, and T at each site in the bam file.

[0063] S13. Calculate the base error rate of each site according to the frequency of the base at that site. The calculation formula is as follows:

[0064] Correct base number: the frequency of cytosine occurrence in the site;

[0065] Number of incorrect bases: the sum of the frequencies of the four bases at the site minus the number of correct bases;

[0066] Base error rate = number of incorrect bases / sum of the frequencies of occurrence of the four bases.

[0067] S14. For each site, statistical tests were performed based on the number of correct bases and incorrect bases in the Treat sample and the Control sample. The specific statistical test method is as follows:

[0068] a. If each sample has multiple replicates, the generalized linear mixed-effects model GLMM (Generalized Linear Mixed-effects Model) is used for statistical testing; GLMM combines the characteristics of the generalized linear model GLM (Generalized Linear Models) and the linear mixed model LMM (Linear Mixed Model), and can handle both fixed effects and random effects. Through GLMM, data with a replication hierarchy can be better processed, especially longitudinal data (here refers to repeated samples). In the case of ac4C modification site identification, the fixed effect part of GLMM refers to cytosine information (whether it is cytosine) and grouping information (Treat or Control), the former is the dependent variable, and the latter is the independent variable. The random effect part refers to the sample information, and each sample has a random intercept.

[0069] b. If each sample has only one replicate, use the chi-square test for statistical testing. The chi-square test is a commonly used nonparametric statistical method that is suitable for analyzing categorical data. By calculating the chi-square statistic, degrees of freedom, and P value, it can be determined whether there is a significant relationship between categorical variables.

[0070] S15. Filter the tested data to obtain potential ac4C modification sites with statistical significance, write the filtered data into an output file, and complete the identification of ac4C modification sites.

[0071] S2. Based on the ac4C modification site data, a difference analysis is performed between the treatment group and the control group to obtain the differential ac4C modification sites. The specific steps are as follows:

[0072] S21, merging the data of the treatment group and the control group. There may be missing values ​​in the merged data, so all the missing values ​​are filled in the merged data;

[0073] S22, standardizing the filled data according to the scaling factor;

[0074] S23. For the standardized data, different methods are used to perform difference tests on the treatment group and the control group, depending on whether there are repeated samples. The specific statistical test methods are as follows:

[0075] a. If each sample has multiple replicates, the generalized linear mixed effects model GLMM is used for difference test; in the case of differential ac4C modification site identification, the fixed effect part in GLMM refers to cytosine information (whether it is cytosine) and group information (treatment group or control group), the former is the dependent variable, and the latter is the independent variable. The random effect part refers to the sample information, and each sample has a random intercept.

[0076] b. If there is only one replicate for each sample, Fisher's exact test was used for difference test.

[0077] S24, filtering the tested data to obtain potential differential ac4C modification sites with statistical significance, writing the filtered data into an output file, and completing the identification of differential ac4C modification sites.

[0078] S3, processing and analyzing the ac4C modification site data (obtained in step S1), and drawing the distribution characteristics of the modification sites; annotating and analyzing the differential ac4C modification site data (obtained in step S2), and obtaining biological process information related to ac4C modification, the specific steps are as follows:

[0079] S31. Analyze the distribution characteristics of ac4C modification sites on chromosomes and gene elements. The specific analysis steps are as follows:

[0080] S311, ac4C modified sites distribution characteristics on chromosomes. First, according to the length information of the chromosome, it is divided into continuous intervals of the same length. Finally, the number of sites with ac4C modification in each interval is counted, and the statistical results are plotted to draw a distribution characteristic diagram.

[0081] S312, ac4C modification site distribution characteristics on gene elements. First, the genome annotation gff file of the species was parsed to obtain the basic information of each element (UTR5, CDS, UTR3), and then the distribution of ac4C modification sites on the elements was counted and plotted.

[0082] S32. Gene annotation was performed based on the differential information of ac4C modification sites to obtain genes related to ac4C modification; then, enrichment analysis was performed using this gene information, including GO (Gene Ontology) enrichment analysis and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis, to identify biological processes related to ac4C modification. Through this process, the role of ac4C modification in cell function and metabolic pathways can be revealed.

[0083] The present invention provides a simple method for identifying ac4C modification sites for chemical ac4C-seq data, improves the resolution of ac4C modification site identification to single base, and provides a method for identifying ac4C modification sites with differences between groups. For experiments with repeated samples, a generalized linear mixed effects model (GLMM) is used to test the data, taking into account both the fixed effects of bases and groupings and the random effects between different samples. In addition, the distribution characteristics of the identified ac4C modification sites are analyzed and the role of ac4C modification in cell function and metabolic pathways is discovered through enrichment analysis. Finally, the method for identifying ac4C modification sites is written in Julia language, with fast calculation speed and high efficiency.

[0084] Specific examples of the present invention are as follows:

[0085] T1. Parse the reference genome GFF file, extract the exon region location information and construct an interval tree. Figure 2 As shown, specifically including:

[0086] T1.1. Extract the exon region location information from the reference genome GFF file and save it in bed6 file format as the region that needs to be retained for subsequent identification of SNV sites.

[0087] T1.2, read the exon area location information file generated in T1.1, such as Fig. 9 As shown, an interval tree containing all exon regions of the reference genome is constructed based on the chromosome name in the first column, the interval start site in the second column, and the interval end site in the third column.

[0088] T2. Count the base type frequency of each site in the bam files of the Treat sample and the Control sample. Figure 3 As shown, specifically including:

[0089] T2.1, Bam file records the position and sequence information of each read aligned to the reference genome. For cDNA sequences, each site has only four bases: A, T, C, and G. Therefore, the number of reads aligned to the four types of bases at each site of the reference genome can be obtained through the bam file, and the original base type of each site can be obtained according to the fasta sequence file of the reference genome. Use the pileup mode of the software JACUSA_v2.0.0, input the bam file and the reference genome fasta file of each sample treatment group and control group respectively, select FR-SECONDSTRAND as the library type, and use the default settings for the rest of the parameters. The result file contains at least 9 columns of information, of which the first six columns are the chromosome name where the site is located, the starting position, the ending position, the data type, the number of reads aligned to the site, and the chain information where the site is located. Starting from the seventh column, the number of A, C, G, and T bases matched at the site by each biological replicate of the sample is shown in turn, and the number of bases of each type is separated by commas. For example, sample A has two biological replicates A1 and A2, and the value of the seventh column is "3,411,2,57", indicating that A1 matched 3 A bases, 411 C bases, 2 G bases, and 57 T bases at the site, respectively. Similarly, the eighth column indicates the number of A, C, G, and T bases matched by A2 at the site. The second to last column of the file stores some descriptive information of the site, and the last column stores the base type of the site on the reference genome.

[0090] T2.2. Save the statistical data of each sample as a pileup file for subsequent analysis.

[0091] T3, read and process the pileup data of each sample. Figure 4 As shown, the specific steps are as follows:

[0092] T3.1, read the pileup files of the Treat and Control samples generated in T2 and store them in Julia DataFrame format. The first three columns of the DataFrame record the chromosome name, start position, and end position of the site in turn. This information is compared with the interval tree constructed in T1, and only the sites on the exon region of the reference genome are retained for subsequent analysis.

[0093] T3.2. According to the positive and negative strand information of each line of the site in the file generated in T3.1 and the base type of the site on the reference genome, the actual base information based on the positive strand is obtained and added to the end of the line. The specific logic is as follows: when the site is on the positive strand, the actual base is the base type of the site on the reference genome; when the site is on the negative strand, the actual base is the complementary base of the base type of the site on the reference genome.

[0094] T3.3. Based on the actual base information of the site obtained in T3.2 and the frequency information of each base at the site, the single nucleotide variation information of the site can be calculated and added to the end of the line, such as Fig.10 As shown. For example: the frequency information of each base is "3,411,2,57", which corresponds to A, C, G, T, and the actual base information is C. Arrange the base frequencies in descending order, remove the frequency information corresponding to the actual base, take the maximum value of the remaining value, and the base corresponding to the maximum value is the single nucleotide variation information of the site. For the convenience of recording and viewing, the base is spliced ​​after the actual base when recording. In this example, the single nucleotide variation information of the site is CT.

[0095] T3.4. Calculate the number of correct bases, the number of incorrect bases and the base error rate based on the frequency information of each base, and add these three types of information to the end of the line. For example, the frequency information of each base is "3,411,2,57", which corresponds to A, C, G, T and the actual base information is C. It can be obtained that the number of correct bases is 411, the number of incorrect bases is 62, and the base error rate is 13.1%.

[0096] T3.5. Combine and integrate the data of each sample output file obtained in T3.4 according to the first six columns. For the integrated data, calculate the median base error rate of each site in the Treat sample and the Control sample, and add the median difference information to the end of the line. Then filter the data of the updated file according to four conditions: 1. The maximum number of incorrect bases in each sample in the Treat sample must be greater than minpileup (default value 3); 2. The minimum base error rate of each sample in the Treat sample must be greater than minmis (default value 0.02); 3. The maximum base error rate of each sample in the Control sample must be less than maxmis (default value 0.05); 4. The median difference in base error rate between the Treat sample and the Control sample must be greater than mindiff (default value 0.02).

[0097] T4. Statistical test to obtain statistically significant modification sites. Figure 5 As shown, the specific steps are as follows:

[0098] T4.1. When there are no repeated samples, use the correct base numbers and incorrect base numbers of the Treat sample and the Control sample at each site to construct a 2×2 matrix. Each row corresponds to the information of each sample. The first column shows the correct base number information of each sample, and the second column shows the incorrect base number information of each sample. For example, the correct base number and incorrect base number of the Treat1 sample are 20 and 5, respectively, and the correct base number and incorrect base number of the Control1 sample are 10 and 2, respectively. Then the first row of the constructed 2×2 matrix is ​​(20, 5), and the second row is (10, 2). Then use the chi-square test to test the constructed matrix, and add the pvalue of the test result to the end of the row.

[0099] T4.2. When there are repeated samples, an n×3 matrix (n is related to the number of samples and bases) is constructed based on the correct base numbers and incorrect base numbers of the Treat samples and Control samples at each site. The first column shows the grouping information (Treat group or Control group), the second column shows the base information (0 represents the correct base, 1 represents the incorrect base), and the third column shows the sample name information (e.g. Treat1, Treat2, Control1, Control2). It is necessary to include all the base information of the sample. For example, when the correct base numbers and incorrect base numbers of the sample Treat1 are 2 and 1 respectively, the generated data format is as follows: (Treat, 0, Treat1), (Treat, 0, Treat1), (Treat, 1, Treat1). The data of all samples are formatted in turn, and then the converted data are tested using the generalized linear mixed effects model (GLMM), and the pvalue of the test result is added to the end of the row.

[0100] T4.3. Filter the data generated by T4.1 or T4.2 according to the pvalue (default pvalue < 0.05). The sites corresponding to each row obtained by screening are potential ac4C modification sites with statistical significance. Write the filtered data into the output file to complete the identification of ac4C modification sites.

[0101] T5. Perform differential analysis on the ac4C modification sites in the treatment group and the control group. Figure 6 As shown, the specific steps are as follows:

[0102] T5.1 merges, fills missing values, and standardizes the data generated by T4.3. The data of the treatment group and the control group are merged according to the first six columns, and two merging modes are supported: intersection merge and union merge. If there are 0 values ​​or missing values ​​in the merged base error rate data, the 0 values ​​and missing values ​​are modified and filled. The default value of modification and filling is 0.01. Finally, the base error rate data of all samples are scaled, and the default value of scaling is 10000.

[0103] T5.2. When there are no repeated samples, extract the base error rate of the treated group Treat sample (recorded as ttmisrate), the base error rate of the treated group Control sample (recorded as tcmisrate), the base error rate of the control group Treat sample (recorded as ctmisrate) and the base error rate of the control group Control sample (recorded as ccmisrate) in each row of the processed data. The base error rates here all refer to the scaled values. Then use the Fisher Exact Test to calculate, the calculation formula is FisherExactTest(ttmisrate, tcmisrate, ctmisrate, ccmisrate), and get the pvalue of the result. Further calculate the log2fc value, the calculation formula is log2((ttmisrate / tcmisrate) / (ctmisrate / ccmisrate)). Add the pvalue and log2fc values ​​to the end of the row.

[0104] T5.3. When there are repeated samples, extract the base error rate of the treated group Treat sample (recorded as ttmisrate-i, i is the i-th repeated sample), the base error rate of the treated group Control sample (recorded as tcmisrate-i), the base error rate of the control group Treat sample (recorded as ctmisrate-i) and the base error rate of the control group Control sample (recorded as ccmisrate-i) in each row of the processed data. The base error rates here all refer to the scaled values. Based on this information, construct an n×3 matrix (n is related to the number of samples and the base error rate after treatment). The first column shows the group information (treated group or control group), the second column shows the base error rate information (Control sample is recorded as 0, Treat sample is recorded as 1), and the third column shows the sample information. For example, in the treated group (W), the Treat samples are Treat1 and Treat2, and the Control samples are Control1 and Control2. Assuming that the base error rate of the Treat1 sample is 1 and the base error rate of the Control1 sample is 2, the constructed data format is as follows: (W, 1, W1), (W, 0, W1), (W, 0, W1). Treat2 and Control2 samples and samples in the control group are processed in the same way, and then the processed data are tested using a generalized linear mixed effects model (GLMM) to obtain the pvalue of the result. The log2fc value is further calculated using the formula: log2(sum(ttmisrate-i / tcmisrate-i) / sum(ctmisrate-i / ccmisrate-i)). Add the pvalue and log2fc value to the end of the line.

[0105] T5.4: Filter the data generated by T5.3 or T5.4 according to the pvalue (default pvalue < 0.05) and log2fc (log2fc > 1.0) values. The sites corresponding to each row obtained by screening are potential differential ac4C modification sites with statistical significance. Write the filtered data into the output file to complete the identification of differential ac4C modification sites.

[0106] T6. Analyze the distribution characteristics of the identified ac4C modification sites. Figure 7 As shown, the specific steps are as follows:

[0107] T6.1. Analyze the distribution characteristics of ac4C modification sites on chromosomes. According to the length information of chromosomes, each chromosome is divided into several intervals according to the specified interval length (the default value is 1000000). For example, if the total length of chr1 chromosome is 248956422, the result of chr1 chromosome segmentation according to the default interval length is: [1, 1000000], [1000001, 2000000], ..., [248000001, 248956422].

[0108] T6.2. According to the location information of the identified ac4C modification sites, the number of ac4C modification sites in each interval of each chromosome is counted. Finally, according to the chromosome length information and the number of ac4C modification sites in each interval of the chromosome, a distribution characteristic map of ac4C modification sites on the chromosome is drawn.

[0109] T6.3, analyze the distribution characteristics of ac4C modification sites on gene elements. First, extract the information of all specified gene elements (UTR5, CDS, UTR3) according to the gff file. Different from the T6.1 analysis method, here we use the position information of ac4C modification sites to determine which gene element each ac4C modification site is distributed on (UTR5, CDS, UTR3 or none of them).

[0110] T6.4. Determine the relative position of the ac4C modification site on the element. For example, suppose there are three ac4C modification site information: [chr1, 500], [chr1, 1300], [chr1, 2000], and it is known that 450-700 on the chr1 chromosome is the UTR5 region, 1200-1600 is the CDS region, and 1900-2100 is the UTR3 region. Considering the upstream and downstream relationship of the gene element (UTR5 is upstream of CDS and UTR3 is downstream of CDS), these three modification sites are respectively on the UTR5, CDS, and TUR3 elements, and the relative positions on their respective corresponding gene elements are 0.2, 1.25, and 2.5. Finally, a density distribution map is drawn based on the annotated gene element information and relative position information of each ac4C modification site.

[0111] T7. Annotate the genes of the identified differentially modified ac4C sites and perform enrichment analysis. Figure 8 As shown, the specific steps are as follows:

[0112] T7.1. Gene annotation of differentially modified sites was performed based on the gff file information and the location information of the identified differentially modified sites. Then, enrichment analysis was performed based on the annotated gene information, including GO enrichment and KEGG enrichment. The GO (Gene Ontology) database describes gene functions consistently, develops a controllable vocabulary, and is species-free. Currently, three independent ontology vocabularies have been established: biological process, cellular component, and molecular function. The KEGG (Kyoto Encyclopedia of Genes and Genomes) database is a bioinformatics database established in 1995 by the Kanehisa Laboratory of the Bioinformatics Center of Kyoto University, Japan. The core of these are the KEGG PATHWAY and KEGG ORTHOLOGY databases. In the KEGG PATHWAY database, biological metabolic pathways are divided into six categories: Cellular Processes, Environmental Information Processing, Genetic Information Processing, Human Diseases, Metabolism, and Organismal Systems.

[0113] T7.2. Use hypergeometric test for enrichment analysis. The hypergeometric test pvalue is calculated as follows:

[0114]

[0115] Among them, N is the number of genes with pathway / GO among all genes; M is the number of genes annotated as a specific pathway / GOterm among all genes; n is the number of genes associated with differential ac4C modification sites in N; m is the number of genes associated with differential ac4C modification sites annotated as a specific pathway / GOterm.

[0116] T7.3. Based on the results calculated in T7.2, use the pvalue (default value is 0.05) to filter and obtain statistically significant pathways / GO terms, and draw bar charts and bubble charts.

[0117] It should be noted that, in this article, the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, device, article or method including a series of elements includes not only those elements, but also includes other elements not explicitly listed, or also includes elements inherent to such process, device, article or method. In the absence of further restrictions, an element defined by the sentence "includes a ..." does not exclude the presence of other identical elements in the process, device, article or method including the element.

[0118] The above description is only a preferred embodiment of the present invention, and does not limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made by using the contents of the present invention specification and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A method for identifying ac4C modification sites, characterized in that: include: Based on the data from ac4C-seq high-throughput sequencing, the treated and controlled groups were obtained, and the ac4C modification sites were identified at single-base resolution. Based on the ac4C modification site data, performing difference analysis between the treatment group and the control group to obtain differential ac4C modification sites; Processing and analyzing the ac4C modification site data, and drawing the distribution characteristics of the modification sites; The differential ac4C modification site data are annotated and analyzed to obtain biological process information related to ac4C modification.

2. The method for identifying ac4C modification sites according to claim 1, characterized in that: The step of obtaining the treatment group and the control group based on the ac4C-seq high-throughput sequencing data and identifying the ac4C modification site at a single-base resolution comprises: The data from the ac4C-seq high-throughput sequencing were quality controlled, the set sequences were removed and the PCR bias was corrected. The quality-controlled data were compared with the reference genome of the corresponding species using STAR software to obtain the bam file of the comparison results of each sample after and without chemical reduction treatment; the comparison results after chemical reduction treatment were used as the treatment group, and the comparison results without chemical treatment were used as the control group; The pileup mode of JACUSA software was used to count the frequency of occurrence of four bases, A, C, G, and T, at each site in the bam file; According to the frequency of bases at each site, the base error rate of the site is calculated; For each site, statistical tests were performed based on the number of correct bases and incorrect bases in the Treat sample and the Control sample; The tested data were filtered to obtain statistically significant potential ac4C modification sites, and the filtered data were written into the output file to complete the identification of ac4C modification sites.

3. The method for identifying ac4C modification sites according to claim 2, characterized in that: In the step of calculating the base error rate of each site according to the frequency of the base at the site, the calculation formula is: Correct base number: the frequency of cytosine occurrence in the site; Number of incorrect bases: the sum of the frequencies of the four bases at the site minus the number of correct bases; Base error rate = number of incorrect bases / sum of the frequencies of occurrence of the four bases.

4. The method for identifying ac4C modification sites according to claim 3, characterized in that: In the step of performing statistical testing on each site based on the number of correct bases and the number of incorrect bases in the Treat sample and the Control sample, the statistical testing method is: If each sample had multiple replicates, the generalized linear mixed effects model GLMM was used for statistical testing; If there was only one replicate for each sample, the chi-square test was used for statistical testing.

5. The method for identifying ac4C modification sites according to claim 1, characterized in that: The step of performing difference analysis between the treatment group and the control group based on the ac4C modification site data to obtain differential ac4C modification sites comprises: Merging the data of the treatment group and the control group, and filling all missing values ​​in the merged data; Normalize the padded data according to the scaling factor; For the standardized data, different methods are used to perform difference tests on the treatment group and the control group according to whether there are repeated samples; The tested data were filtered to obtain potential differential ac4C modification sites with statistical significance, and the filtered data were written into the output file to complete the identification of differential ac4C modification sites.

6. The method for identifying ac4C modification sites according to claim 5, characterized in that: In the step of performing a difference test on the treatment group and the control group using different methods according to whether there are repeated samples for the standardized data, the test method is: If each sample had multiple replicates, the generalized linear mixed effects model GLMM was used for difference testing; If there was only one replicate for each sample, Fisher's exact test was used for difference testing.

7. The method for identifying ac4C modification sites according to any one of claims 4 and 6, characterized in that: The generalized linear mixed effects model GLMM handles both fixed effects and random effects; The fixed effect part refers to cytosine information and grouping information, the former is the dependent variable and the latter is the independent variable; wherein the grouping includes Treat and Control; The random effect part refers to the sample information, and each sample has a random intercept.

8. The method for identifying ac4C modification sites according to claim 1, characterized in that: The ac4C modification site data are processed and analyzed to draw the distribution characteristics of the modification sites; The step of annotating and analyzing the differential ac4C modification site data to obtain biological process information related to ac4C modification comprises: Analyze the distribution characteristics of ac4C modification sites on chromosomes and gene elements; Gene annotation was performed based on the differential information of ac4C modification sites to obtain genes related to ac4C modification; The genes are used to perform enrichment analysis to identify biological processes related to ac4C modification; wherein the enrichment analysis includes GO enrichment analysis and KEGG enrichment analysis.

9. The method for identifying ac4C modification sites according to claim 8, characterized in that: The step of analyzing the distribution characteristics of ac4C modification sites on chromosomes and gene elements comprises: According to the length information of the chromosome, it is divided into continuous intervals of the same length. Finally, the number of sites modified by ac4C in each interval is counted, and the statistical results are drawn into a distribution characteristic diagram to analyze the distribution characteristics of ac4C modified sites on the chromosome. The genome annotation gff file of the species was parsed to obtain the basic information of each element, and then the distribution of ac4C modification sites on the elements was counted and mapped to analyze the distribution characteristics of ac4C modification sites on gene elements; wherein the elements include UTR5, CDS, and UTR3.

Images