A method and system for analyzing m6a high-throughput sequencing data

Abstract

The application discloses a kind of m6A high-throughput sequencing data bioinformation analysis methods, belong to bioinformatics technical field.The analysis method includes original data output statistics and pre-processing, data pollution filtration, sequencing data and genome alignment, whole genome level Peak Calling and annotation, difference Peak analysis between samples, gene sequence motif analysis, gene and lncRNA expression abundance analysis, differential gene expression profile analysis, differential lncRNA expression profile analysis.Using the analysis method provided by the application, comprehensive detailed analysis data can be obtained, and m6A methylation modification information of animals, plants and clinical samples can be obtained from whole genome level using high-throughput sequencing means, to provide help for subsequent scientific research.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and specifically relates to an analysis method and system for m6A high-throughput sequencing data. Background Art

[0002] RNA modification, a significant biological phenomenon following DNA and histone modifications, can influence gene expression and protein synthesis without altering the RNA sequence. Over 170 modifications are known to occur on RNA, among which m6A is one of the most abundant RNA base modifications and has emerged in recent years as a novel regulatory mechanism for controlling gene expression in eukaryotic cells. This methyl modification, occurring at the sixth position of the nitrogen adenine (A) base, is widely found in messenger RNA (mRNA), long noncoding RNA (lncRNA), microRNA (miRNA), circular RNA (circRNA), and ribosomal RNA (rRNA). As a reversible epigenetic modification, m6A influences the fate of modified RNA molecules and plays a crucial role in nearly all important biological processes.

[0003] LncRNAs are a group of endogenous RNA molecules that are not translated into proteins but have specialized functions in regulating gene expression. They can be divided into lncRNAs (those longer than 200 nt) and small ncRNAs (those less than 200 nt). In addition to protein-coding mRNAs, m6A modifications have been found in ncRNAs, including lncRNAs, circRNAs, and miRNAs, and have been shown to be important for the expression and function of these RNAs.

[0004] In recent years, m6A modification has expanded from basic scientific research to a variety of downstream applications. Research in nearly every field of medicine and agriculture, including the development of targeted tumor cells, immunotherapy inhibitors, and increasing grain yields, requires researchers to initially detect potential RNA modification information in samples. High-throughput sequencing methods can detect modification information across all RNA in a sample. The m6A-seq method, based on the m6A antibody IP principle, is currently the most common experimental approach, generating a large amount of m6A-seq data daily. Therefore, a bioinformatics analysis method that can simultaneously mine m6A modifications in both mRNA and lncRNA is particularly important. Summary of the Invention

[0005] In order to solve the above technical problems, the present invention aims to provide a bioinformatics analysis method for m6A-seq high-throughput sequencing that can simultaneously analyze mRNA and lncRNA. To this end, the technical solution adopted by the present invention is as follows:

[0006] The first aspect of the present invention provides a method for analyzing m6A-seq high-throughput sequencing data, comprising the following steps:

[0007] S1, m6A-seq high-throughput sequencing data is processed to filter out unqualified sequences, wherein the m6A-seq sequencing data comes from at least one pair of paired IP samples and input samples. IP refers to the abbreviation of immunoprecipitation in m6A sequencing, and input is a positive control in m6A experiments. Input samples are not subjected to IP experiments, which can help determine the IP effect of m6A. By comparing input and IP, the degree of m6A modification and the experimental effect can be determined;

[0008] S2, filtering out contaminating sequences using a database including ribosomal RNA sequences and / or mycoplasma RNA sequences;

[0009] S3, aligning the filtered sequences with the reference genome of the species from which the samples were derived, and obtaining reads that can be aligned to the reference genome;

[0010] S4, genome-wide peak calling and annotation: perform genome-wide peak identification on reads that can be aligned to the reference genome, obtain peak positions on the genome, peak lengths, and read count information, annotate peak positions on the genome, and analyze the distribution of reads at transcription start and end positions;

[0011] S5, differential peak analysis: Take the peak of the intersection of IP sample and input sample for differential peak analysis to obtain the peak with differences. The related genes corresponding to the peak with differences are analyzed by GO and KEGG functional enrichment using the negative binomial distribution method;

[0012] S6, Peak sequence motif analysis: Perform motif analysis on the single sample Peak obtained in step S4 and the Peak with differences obtained in step S5;

[0013] S7, mRNA and lncRNA expression abundance analysis: Quantitative analysis of mRNA and lncRNA in the sample was performed to obtain the original counts and normalized FPKM value of each mRNA and lncRNA in the sample;

[0014] S8. Differential mRNA and lncRNA expression profiling analysis: Differential analysis was performed on each mRNA and lncRNA in the IP sample and the input sample, and GO and KEGG functional enrichment analysis was performed on all differentially associated genes using the negative binomial distribution method.

[0015] In some embodiments of the present application, in step S1, the unqualified sequences include: ① reads with adapters; ② reads containing more than 5% of N bases; ③ reads with more than 20% of bases having a quality value Q≤10; and ④ reads with all bases being A or N.

[0016] In some embodiments of the present application, in step S4, the method for identifying the Peak is as follows:

[0017] S41, scanning the exon region to obtain a candidate Peak region;

[0018] S42, assuming that a gene is located in an exon window w=[w1, w2], calculating the number of reads of the IP and Input samples aligned to the window w=[w1, w2]:

[0019] X 1,w ~ Poisson(n1λ 1,w ),

[0020] X 0,w ~ Poisson(n0λ 0,w ),

[0021] wherein, X 1,w and X 0,w are the number of reads of the IP and Input samples aligned to the window w=[w1, w2], n1 and n0 are the total number of reads of the IP and Input samples, λ 1,w and λ 0,w represent the probability of the sequencing reads of the IP and Input samples falling into the window w, respectively;

[0022] S43, using C-test to compare the average value ratio c=λ1,w / λ0,w of the two Poisson distributions X1,w~Poisson(n1λ1,w) and X0,w~Poisson(n0λ0,w) of the candidate Peak region, and calculating the p-value;

[0023] S44, using Fisher test to calculate the p-value of the candidate Peak region, and the region with a p-value less than 0.05 is the Peak.

[0024] In some embodiments of the present application, in step S7, the FPKM value is calculated according to the following formula:

[0025]

[0026] Among them, total exon fragments indicates the total number of exon fragment reads, mappedreads indicates the number of millions of reads that can be mapped to the genome, and exon length indicates the kb length of the exon.

[0027] The analysis method according to claim 1, characterized in that in step S8, the GO and KEGG functional enrichment analysis is performed using the negative binomial distribution method, and the formula is:

[0028]

[0029] Among them, pvalue of enriched GO / KEGG represents the GO / KEGG functional enrichment analysis value, TB represents all background genes, TS represents all differentially expressed genes, B represents all genes with GO or KEGG annotations, and S represents all differentially expressed genes with GO or KEGG annotations.

[0030] A second aspect of the present invention provides an analysis system for m6A-seq sequencing data, comprising the following modules:

[0031] Data input module, used to accept m6A-seq high-throughput sequencing data input and perform preprocessing;

[0032] The alignment module and the data input module are used to align the preprocessed sequence with the reference genome of the sample source species to obtain reads that can be aligned to the reference genome;

[0033] The Peak Calling and Annotation module is connected to the alignment module and is used to perform genome-wide peak identification on reads that can be aligned to the reference genome, obtain peak positions on the genome, peak lengths, and read quantity information, annotate peak positions on the genome, and analyze the distribution of reads at transcription start and end positions;

[0034] The differential peak analysis module is connected to the Peak Calling and Annotation module to perform differential peak analysis on the peaks at the intersection of the IP sample and the Input sample to obtain differential peaks, and to perform GO and KEGG functional enrichment analysis on the related genes corresponding to the differential peaks using the negative binomial distribution method;

[0035] Peak sequence motif analysis module, connected to the difference peak analysis module, for performing motif analysis on the obtained single sample peaks and peaks with differences;

[0036] The mRNA and lncRNA expression abundance analysis module is connected to the comparison module and is used to quantitatively analyze the mRNA and lncRNA in the sample to obtain the original value Counts and normalized expression quantity FPKM value of each mRNA and lncRNA in the sample;

[0037] The differential mRNA and lncRNA expression profile analysis module is connected to the mRNA and lncRNA expression abundance analysis module to perform differential analysis on each mRNA and lncRNA in the IP sample and the Input sample, and perform GO and KEGG functional enrichment analysis on all differentially related genes using the negative binomial distribution method.

[0038] In some embodiments of the present invention, preprocessing in the data input module refers to filtering out unqualified sequences and / or contaminating sequences. The unqualified sequences include: ① reads with adapters; ② reads containing more than 5% N bases; ③ reads with a quality value Q ≤ 10 where the number of bases accounts for more than 20% of the entire read; and ④ reads where the bases are all A or N. Contaminating sequences are filtered out using a database including ribosomal RNA sequences and / or mycoplasma RNA sequences.

[0039] In some embodiments of the present invention, in the Peak Calling and Annotation module, the Peak identification method is:

[0040] S41, scan the exon region to obtain the candidate Peak region;

[0041] S42, assuming the gene is located in the exon window w = [w1, w2], calculate the number of reads that the IP and Input map to the window w = [w1, w2]:

[0042] X1,w~Poisson(n1λ1,w),

[0043] X0,w~Poisson(n0λ0,w),

[0044] Where X1,w and X0,w are the number of reads of the IP sample and the Input sample aligned to the window w = [w1,w2], n1 and n0 are the total number of reads of the IP sample and the Input sample, respectively, and λ1,w and λ0,w are the probabilities that the sequencing reads of the IP sample and the Input sample fall into the window w, respectively.

[0045] S43, using C-test to compare the ratio c = λ1,w / λ0,w of X1,w ~ Poisson(n1λ1,w) and X0,w ~ Poisson(n0λ0,w) of the candidate Peak region, and calculate the p-value;

[0046] S44, using Fisher test to calculate the p-value of the candidate Peak region, and the region with a p-value less than 0.05 is the Peak.

[0047] In some embodiments of the present application, in the mRNA and lncRNA expression abundance analysis module, the FPKM value calculation formula is as follows:

[0048]

[0049] Wherein, total exon Fragments represents the total number of exon fragment reads, Mapped reads represents the number of reads that can be aligned to the genome, and exon length represents the kb length of the exon.

[0050] In some embodiments of the present application, in the differential mRNA and lncRNA expression profile analysis module, the GO and KEGG function enrichment analysis is performed using the negative binomial distribution method, and the formula is as follows:

[0051]

[0052] Wherein, pvalue of enriched GO / KEGG represents the GO / KEGG function enrichment analysis value, TB represents all background genes, TS represents all differential expression genes, B represents all genes with GO or KEGG annotation, and S represents all differential expression genes with GO or KEGG annotation.

[0053] Advantages of the present application

[0054] Compared with the prior art, the present application has the following advantages:

[0055] (1) Using the method or system of the present application, the m6A modification information and expression information on mRNA and lncRNA can be detected simultaneously from high-throughput sequencing data obtained by an experimental method, and the data can be fully mined from four dimensions

[0056] (2) The analysis method and system of the present application only scan the exon region and do not scan the intron region during Peak Calling, and only report the Peak information of annotated mRNA and lncRNA, so that the information is more targeted.

[0057] (3) The analysis method and system of the present invention outputs the Peak counts value, which will make the Peak abundance value more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 The figure shows the bioinformatics analysis flow chart of the m6A-seq high-throughput sequencing data of lncRNA of the present invention.

[0059] Figure 2 A pie chart showing the proportion of exons and introns in the genome of reads in an embodiment of the present invention is shown.

[0060] Figure 3 The figure shows the distribution of Peak relative gene positions in the embodiment of the present invention.

[0061] Figure 4 A pie chart showing the proportion of Peak in different gene functional elements in an embodiment of the present invention is shown.

[0062] Figure 5 The Peak corresponding motif diagram in an embodiment of the present invention is shown.

[0063] Figure 6 Shown are box plots of gene expression in an embodiment of the present invention.

[0064] Figure 7 Shown is a gene expression density map in an embodiment of the present invention.

[0065] Figure 8 A heat map of differentially expressed genes in an embodiment of the present invention is shown.

[0066] Figure 9 A volcano plot of differentially expressed genes in an embodiment of the present invention is shown.

[0067] Figure 10 GO term functional enrichment histogram in the embodiment of the present invention.

[0068] Figure 11 KEGG pathway enrichment bubble diagram in the embodiment of the present invention. DETAILED DESCRIPTION

[0069] Unless otherwise indicated, all parts and percentages expressed herein are based upon weight and all tests and measurements are conducted synchronously with the filing date of this application. The contents of any patent, patent application, or publication referenced herein are hereby incorporated by reference in their entirety for the teachings relevant to the sentence and / or paragraph in which the reference is presented. Equivalent same family patents to the extent available to the public are also incorporated by reference. In the event that the definition of a term in the prior art is inconsistent with the definition provided herein, the definition provided herein shall control.

[0070] Numerical ranges expressed in the present application as comprising one or more endpoints are approximate ranges which include all values falling within such ranges, unless otherwise indicated. The range of values includes all values within the range, including the values defining the range. For example, a range of "1 to 10" is intended to include all whole integers, fractions, and decimals within the given range. Any numerical range recited herein is intended to include all sub-ranges of the same numerical precision subsumed therewithin. For example, a range of "1 to 10" is intended to include a range of "2 to 5," "3 to 7," "4 to 6," etc. Any reference to a numerical range is intended to include all values subsumed therein.

[0071] The terms "comprising," "including," "containing," and variations thereof, do not exclude the presence of other components, steps or processes, and are used synonymously with the term "including." For purposes of the subject application, the term "consisting essentially of to the extent it is used, means including also any additional component, step or process that does not materially affect the basic and novel characteristics of the composition(s) or method(s) described herein. The term "consisting of excludes any component, step or process not specifically listed. The term "or" as used herein is used in the inclusive sense, i.e., "and / or," unless otherwise stated.

[0072] In order to make the technical problems solved by the present application, the technical solutions and the beneficial effects clearer, the present application will be further described in detail below with reference to the embodiments.

[0073] Example

[0074] The following examples are provided to illustrate preferred embodiments of the present invention. Those skilled in the art will appreciate that the techniques disclosed in the following examples represent techniques discovered by the inventors that can be used to practice the present invention and, therefore, can be considered preferred embodiments of the present invention. However, those skilled in the art will appreciate from this disclosure that many modifications may be made to the specific embodiments disclosed herein while still achieving the same or similar results without departing from the spirit or scope of the present invention.

[0075] Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs, and the disclosure and materials they cite are hereby incorporated by reference.

[0076] Those skilled in the art will recognize, or be able to ascertain using no more than routine experimentation, many technical equivalents to the specific embodiments of the invention described herein. Such equivalents are intended to be encompassed by the claims.

[0077] The experimental methods in the following examples, unless otherwise specified, are all conventional methods. The instruments and equipment used in the following examples, unless otherwise specified, are all conventional laboratory instruments and equipment; the experimental materials used in the following examples, unless otherwise specified, are all purchased from conventional biochemical reagent stores.

[0078] Example 1: Bioinformatics analysis of lncRNA m6A-seq high-throughput sequencing data

[0079] Combine Figure 1 , the detailed analysis process is as follows:

[0080] 1. Raw data output statistics and preprocessing

[0081] The raw data generated by the paired-end sequencing method using the Illumina Novaseq 6000 or HiSeq series sequencers is in the fastq format with adapter sequences. This data format contains base quality and specific base information. These raw data require preliminary preprocessing and use fastq software to filter out unqualified sequences. The specific steps include: ① reads with adapters; ② removal of reads containing N bases (N indicates base information that the sequencer cannot recognize) greater than 5%; ③ removal of low-quality reads (the number of bases with a quality value Q ≤ 10 accounting for more than 20% of the total read); ④ removal of reads with bases that are all A or N.

[0082] Then, the Q20%, Q30%, GC%, effective data volume, effective read number and other data of the sequencing data were counted, as shown in Table 1.

[0083] Table 1 Sequencing sequence quality control statistics

[0084] Sample_ID Raw_Reads Raw_Bases Valid_Reads Valid_Bases Valid_Ratio% Q20％ Q30％ GC% nc1_input 109333796 16.51G 105863002 13.93G 84.38 98.09 94.36 49.25 nc1_IP 121467832 18.34G 114893386 14.95G 81.51 97.98 94.11 49.26

[0085] Note:

[0086] Sample_ID: sample name

[0087] Raw_Reads: The number of reads of the original offline data

[0088] Raw_Bases: The number of bases in the original data.

[0089] Valid_Reads: Number of reads of valid data

[0090] Valid_Bases: The number of bases of valid data

[0091] Valid_Ratio%: The proportion of valid reads

[0092] Q20%: The proportion of bases with a quality value ≥ 20 (sequencing error rate less than 0.01)

[0093] Q30%: The proportion of bases with a quality value ≥ 30 (sequencing error rate less than 0.001)

[0094] GC%: the proportion of GC content

[0095] 2. Contaminated sequence filtering

[0096] (1) Download of rRNA and mycoplasma RNA sequences

[0097] All 5s / 5.8s / 18s / 28s rRNA sequences were downloaded from the NCBI NR database and the Ensembl database. These sequences cover species such as humans, animals, and plants, including Homo sapiens, mouse (Mus musculus), rat (Rattus norvegicus), zebrafish (Danio rerio), Drosophila melanogaster, Arabidopsis thaliala, rice (Oryza sativa), tobacco (Nicotiana tabacum), and wheat (Triticum aestivum). All common mycoplasma RNA sequences were downloaded from the NCBI NR database using the keyword "mycoplasma." All RNA sequences were in fasta format.

[0098] (2) Create In-House ribosome and mycoplasma sequence database

[0099] After downloading all ribosomal RNA sequences from NCBI database and Ensembl database, merge them into a large sequence file for subsequent contamination data filtering.

[0100] (3) Contaminated sequence filtering

[0101] Using bowtie software, aligning with the ribosome and mycoplasma sequence database with default parameters, base mismatch is 3, filtering the relevant contaminated sequences in the sequencing data.

[0102] 3. Sequencing data and genome alignment

[0103] (1) Download the sequence file (fasta or fa format) and annotation file (gtf format) of the reference genome.

[0104] (2) Use extract_exons.py and extract_exons.py in HISAT2 software to convert the gtf file of the reference genome into exon file and ss file, representing exon information and alternative splicing information respectively.

[0105] (3) Use hisat2-build instruction in HISAT2 software to construct genome index file (ht2 format) through genome fa file, exon file and ss file, which is used for later genome alignment.

[0106] (4) Use hisat2 instruction of HISAT2 software to align each read in IP and Input CleanData file (fastq format) with reference genome index file, and the results are shown in Table 2.

[0107] Table 2 Reference genome alignment reads statistics table

[0108]

[0109] Note:

[0110] Sample: sequencing library name

[0111] Valid reads: data after quality control

[0112] Mapped reads: number of reads that can be aligned to the genome

[0113] Unique Mapped reads: number of reads that can only be uniquely aligned to one position of the genome

[0114] Multi mapped reads: The number of reads that can be mapped to multiple locations in the genome

[0115] Reads map to sense strand: The number of reads mapped to the sense strand on the genome

[0116] Reads map to antisense strand: Statistics of reads mapped to the negative sense strand on the genome

[0117] Non-splice reads: Reads can be aligned end-to-end to the genomic region, and statistics are for the entire segment alignment.

[0118] Splice reads: Reads cannot be aligned end-to-end to the genomic region, and are segmented alignment statistics.

[0119] Output the reads that can be mapped to the genome as a Mapped.bam file.

[0120] (5) Draw a pie chart of the alignment rates between IP samples and Input genomes.

[0121] The results are as follows Figure 2 As shown, the results show that the reads are compared to the exon region (exon), intron region (intron) and intergenic region (intergenic) of the genome, and more than 90% of the reads are located in the exon region.

[0122] 4. Peak Calling and Annotation at the Whole-Genome Level

[0123] (1) Peak Calling analysis was performed using exomePeak2 software. The main input files included the bam files of the IP sample and the input sample and the genome annotation gtf file. The main parameters included "glm_type = c("Poisson", "NB", "DESeq2")", "peak_calling_mode = exon", and "fragment_length = 100". This step only scans the exon region and does not scan the intron region, so the Peak Calling mode is only exon. DESeq2 is used by default for differential analysis. Different scan lengths will affect the number of peaks. The shorter the length, the more peaks there are. The default scan length is 100. Due to the heterogeneity of the input sample, the background is estimated by using the position tag density method.

[0124] X1,w and X 0,w is the number of reads of the IP sample and the Input sample aligned to the window w = [w1, w2]. Therefore, the gene is located in the exon window w = [w1, w2]. The reads follow the Poisson distribution. Calculate the two Poisson means, which are X 1,w ~Poisson(n1λ 1,w ) and X 0,w ~Poisson(n0λ 0,w ). X 1,w and X 0,w is the number of reads from IP and Input aligned to the window w = [w1, w2], n1 and n0 are the total number of reads for each sample, and λ 1,w and λ 0,w Represent the probability that the IP and Input sample sequencing reads fall into window w respectively.

[0125] The C-test is then used to compare the mean ratio (c = λ1,w / λ0,w) of the two Poisson equations for the candidate peak region (IP sample X1,w~Poisson(n1λ1,w)) and the input sample X0,w~Poisson(n0λ0,w), and the p-value is calculated. The p-value of the candidate peak region is then calculated using the Fisher test. By default, a region with a p-value less than 0.05 is considered a peak. The output format is bed format, and key information includes chromosome number and region, peak length, gene ID, IP and input sample counts, log2 fold change (Fold Change) value, p-value, and FDR value.

[0126] (2) After reading the bed file of the Peak using ChIPseeker software, the distribution information of the Peak on the chromosome and the distribution information at the transcription start site TSS were statistically analyzed, and the gene functional elements corresponding to the Peak were annotated, including the Promoter region, 5'UTR region, 3'UTR region, and CDS region. The results were output in the form of a statistical table, as shown in Table 3:

[0127] Table 3 Peak calling result output table

[0128]

[0129]

[0130] Note:

[0131] chr: The chromosome number where the mRNA or lncRNA is located. For example, chr1 means that the gene corresponding to the transcript is located on chromosome 1. start: The starting site of the mRNA or lncRNA peak on the chromosome.

[0132] end: The end site of the mRNA or lncRNA Peak on the chromosome

[0133] width: The length of the Peak range, that is, the size of the peak area

[0134] geneId: The ID number corresponding to the gene, generally derived from authoritative databases such as NCBI or ensembl

[0135] gene_biotype: mRNA or lncRNA corresponding to the properties of the encoded protein, such as protein coding, lincRNA, etc.

[0136] geneName: abbreviation or abbreviation of mRNA or lncRNA, also called gene symbol

[0137] GO: full name gene ontology, which means annotation of genes, also called GO functional annotation

[0138] KEGG: small k refers to the name of the KEGG metabolic pathway, which is used to indicate a specific biological pathway

[0139] ReadsCount.Input: The number of reads in the Peak area of ​​the Input sample

[0140] ReadsCount.IP: The number of reads in the Peak region of the IP sample

[0141] log2fc: The multiple of the peak in the IP sample to the peak in the input, when processed with log2

[0142] pval: p-value of the peak and then log10 processing

[0143] FDR: The value after correcting the p-value of the peak using the FDR value

[0144] Annotation: Peak genomic location functional element annotation, i.e., whether the Peak is in exon or 3'UTR, etc.

[0145] (3) Draw the enrichment peak of m6A at the starting site TSS of each sample

[0146] like Figure 3 As shown in the figure, after normalization of all gene distribution reads, more than 50% of the fragments were enriched in the 5'UTR end and the transcription start site TSS at the start of CDS, and more than 80% of the fragments were enriched in the CDS end and the transcription termination site TES at the start of 3'UTR. It is consistent with the characteristics of m6A experiment.

[0147] 5. Differential Peak Analysis

[0148] (1) The exomePeak2 software used in differential Peak analysis has the same analysis parameters as Peak Calling, and the output file is in bed format. The specific content includes the location of differential Peak on the chromosome, the ID number of gene or IncRNA, p-value value, log2FC, Peak gene functional element annotation, Peak fold enrichment, etc. The default p-value < 0.05 and the absolute value of log2FC > 1 are differential Peak. The results are shown in Table 4:

[0149] Table 4 Differential Peak Information Statistics

[0150]

[0151]

[0152] Note:

[0153] Chr: The number of the chromosome where the mRNA or IncRNA is located, such as chr1 means that the gene corresponding to the transcript is located on chromosome 1

[0154] end: The end of the Peak peak on the chromosome of mRNA or IncRNA

[0155] width: The range length of Peak, that is, the size of the peak area

[0156] geneId: The ID number corresponding to mRNA or IncRNA, generally from authoritative databases such as NCBI or ensembl

[0157] gene_biotype: The attribute corresponding to mRNA or IncRNA, such as protein coding, IncRNA, etc.

[0158] geneName: The abbreviation or abbreviation of mRNA or IncRNA, also called gene symbol

[0159] ReadsCount.Input.nc: The number of reads in the Peak region where the Input sample of the nc group is located

[0160] ReadsCount.IP.nc: The number of reads in the Peak region where the IP sample of the nc group is located

[0161] ReadsCount.Input.t1: The number of reads in the Peak region where the Input sample of group t1 is located

[0162] ReadsCount.IP.t1: The number of reads in the Peak region where the IP sample in group t1 is located

[0163] log2fc: The difference in peak values ​​between the nc combination and t1 groups (two groups) is processed using log2

[0164] pval: p-value obtained by Peak difference analysis between nc group and t1 group (two groups) fdr: the value obtained by correcting the p-value obtained by Peak difference analysis between nc group and t1 group (two groups) using fdr value

[0165] Annotation: Peak genomic location functional element annotation, i.e., whether the Peak is in exon or 3'UTR, etc.

[0166] (2) Draw a pie chart showing the distribution of differential peaks in gene functional elements

[0167] like Figure 4 As shown in the figure, over 60% of the peaks are located in the 3'UTR region, also known as UTR3. The second peak, over 13%, is located in UTR5, also known as the 5'UTR. Overall, over 70% of the peaks are located in the UTR region, consistent with the characteristics of m6A experiments.

[0168] 6.Peak sequence motif analysis

[0169] Use the findmotifsGenome module in HOMER. Output files include a Peak bed file and genomic position annotation information files (such as tss, tts, aug, stop, and repeats.rna files). Key parameters include RegionSize = 100 and Motif length = 5. Output the Homer de novo Motif Results webpage.

[0170] The results are as follows Figure 5As shown in the figure, the motif analysis structure is very typical m6A characteristic motif RGACH. Among them, R is G base or A base, and H is A base, C base or U base. And pvalue < 0.05.

[0171] 7. mRNA and lncRNA expression abundance analysis

[0172] (1) Use StringTie software, input file is single sample Input bam file and genome gtf file, and finally output merged.gtf file. According to the transcript_type definition information in gtf file, such as protein_coding, lincRNA, ncRNA and other annotations, first annotate each RNA as mRNA and lncRNA.

[0173] (2) Extract the Counts value and normalized expression FPKM value of each mRNA and lncRNA in merged.gtf file, and then output the expression profile matrix. The FPKM value calculation formula is as follows:

[0174]

[0175] Among them, total exon Fragments represents the total number of exon reads, Mapped reads represents the number of reads that can be aligned to the genome (million), and exon length represents the length of exon (kb).

[0176] The results are shown in Table 5:

[0177] Table 5 mRNA and lncRNA overall expression level

[0178]

[0179] Note:

[0180] gene_id: mRNA or lncRNA corresponding ID number, generally from authoritative databases such as NCBI or ensembl

[0181] gene_name: abbreviation or abbreviation of mRNA or lncRNA, also called gene

[0182] GO: full name gene ontology, that is, gene annotation, also called GO function annotation

[0183] KEGG: small k, refers to the name of KEGG metabolic pathway, used to represent a specific biological path

[0184] KO_Entry: Big K, refers to KEGG entries with similar results or functions in all species

[0185] EC: Generally refers to a specific enzyme code. For example, EC2.7.1.1 refers to serine / threonine kinase

[0186] Description: Detailed description of the mRNA or lncRNA, including its function, etc.

[0187] FPKM.nc1_input: FPKM normalized expression level of mRNA or lncRNA of sample cn1

[0188] FPKM.t1_input: FPKM normalized expression level of mRNA or lncRNA in sample t1

[0189] (3) Based on the FPKM values ​​of genes and lncRNAs, draw gene expression box plots and gene expression density maps.

[0190] Gene expression box plots Figure 6 As shown in the figure, after log10 normalization of the expression levels of all mRNA and lncRNA in the sample FPKM, the overall expression level of the nc1 sample is slightly higher than that of the t1 sample.

[0191] Gene expression density map Figure 7 As shown in the figure, after log10 processing of the FPKM expression of all mRNA and lncRNA in the sample, it becomes the relative expression. Therefore, a regular distribution diagram of the relative expression size can be found, in which the high and low expression distribution patterns of both nc1 and t1 are similar.

[0192] 8. Differential mRNA and lncRNA Expression Profiling

[0193] (1) Using the edgeR software, differential analysis was performed based on the counts of each mRNA and lncRNA in the samples between groups, and the p-value was calculated. The default p-value < 0.05 was considered differentially expressed mRNA and lncRNA. The results are shown in Table 6:

[0194] Table 6. Differential analysis of mRNA and lncRNA

[0195]

[0196] Note:

[0197] gene_id: The ID number corresponding to the mRNA or lncRNA, generally derived from authoritative databases such as NCBI or ensembl

[0198] gene_name: the abbreviation or abbreviation of mRNA or lncRNA, also called gene GO: the full name of gene ontology, which is the annotation of genes, also called GO functional annotation

[0199] KEGG: small k refers to the name of the KEGG metabolic pathway, which is used to indicate a specific biological pathway

[0200] KO_Entry: Big K, refers to KEGG entries with similar results or functions in all species

[0201] EC: Generally refers to a specific enzyme code. For example, EC2.7.1.1 refers to serine / threonine kinase

[0202] Description: Detailed description of the mRNA or lncRNA, including its function, etc.

[0203] FPKM.nc1_input: FPKM normalized expression level of mRNA or lncRNA of sample cn1

[0204] FPKM.t1_input: FPKM normalized expression level of mRNA or lncRNA in sample t1

[0205] fc: fold change

[0206] pval: p-value obtained after difference analysis between the two groups qval: q value obtained by FDR correction of p-value

[0207] (2) Draw differential gene heat maps and volcano maps based on the expression levels and p-values ​​of differentially expressed genes

[0208] Heat map of differentially expressed genes Figure 8 As shown in the figure, the relative expression levels of the two sample groups were compared after performing log10 processing on the FPKM expression levels. The relative expression levels of each gene in the two sample groups are correlated with color. The closer the color of a gene is to dark blue, the lower the relative expression level, while the closer the color is to dark red, the higher the relative expression level.

[0209] Volcano plot of differentially expressed genes Figure 9 As shown in the figure, the pvalue and fold difference obtained from the difference analysis of a gene between the two groups of samples are processed by -log10 and log2, respectively. Each point represents a gene, and is displayed in a dual quadrant plot on the X and Y axes. The small point in the upper left corner outside the dotted line indicates a significant decrease in gene expression, while the small point in the upper right corner outside the dotted line indicates a significant increase in gene expression.

[0210] (3) For all differentially correlated genes, GO and KEGG functional enrichment analysis was performed using the negative binomial distribution method based on their corresponding GO and KEGG annotation information. The formula is:

[0211]

[0212] Among them, pvalue of enriched GO / KEGG represents the GO / KEGG functional enrichment analysis value, TB represents all background genes, TS represents all differentially expressed genes, B represents all genes with GO or KEGG annotations, and S represents all differentially expressed genes with GO or KEGG annotations.

[0213] (4) If the p-value is less than 0.05, the GO or KEGG entry is considered to be significantly enriched. Finally, a bar chart and bubble chart of the GO and KEGG enrichment analysis are drawn.

[0214] The histogram of GO enrichment analysis is as follows Figure 10 As shown in the figure, it can be seen that the higher the number of genes annotated with a certain gene function, the higher the corresponding column, which means that the enrichment degree of this gene function is higher.

[0215] The bubble diagram of KEGG enrichment analysis is as follows Figure 11 As shown in the figure, the corresponding gene pathway is also called pathway. The more genes are enriched, the larger the bubble is. The smaller the pvalue of the corresponding pathway enrichment, the closer the color is to blue, which means that the degree of enrichment of the pathway is higher.

[0216] All documents mentioned in this application are incorporated herein by reference, just as if each document were incorporated herein by reference individually. It should also be understood that after reading the above teachings of the present invention, those skilled in the art may make various changes or modifications to the present invention, and that such equivalents also fall within the scope of the claims appended hereto.

Claims

1. A method for analyzing m6A-seq high-throughput sequencing data, characterized in that: The following steps are involved: S1, processing the m6A-seq high-throughput sequencing data to filter out unqualified sequences, wherein the m6A-seq sequencing data comes from at least one pair of paired IP samples and input samples; S2, filtering out contaminating sequences using a database including ribosomal RNA sequences and / or mycoplasma RNA sequences; S3, aligning the filtered sequences with the reference genome of the species from which the samples were derived, and obtaining reads that can be aligned to the reference genome; S4, genome-wide peak calling and annotation: perform genome-wide peak identification on reads that can be aligned to the reference genome, obtain peak positions on the genome, peak lengths, and read count information, annotate peak positions on the genome, and analyze the distribution of reads at transcription start and end positions; S5, differential peak analysis: Take the peak of the intersection of IP sample and input sample for differential peak analysis to obtain the peak with differences. The related genes corresponding to the peak with differences are analyzed by GO and KEGG functional enrichment using the negative binomial distribution method; S6, Peak sequence motif analysis: Perform motif analysis on the single sample Peak obtained in step S4 and the Peak with differences obtained in step S5; S7, mRNA and lncRNA expression abundance analysis: Quantitative analysis of mRNA and lncRNA in the sample was performed to obtain the original counts and normalized FPKM value of each mRNA and lncRNA in the sample; S8, Differential mRNA and lncRNA expression profiling analysis: Differential analysis was performed on each mRNA and lncRNA in the IP sample and the Input sample, and GO and KEGG functional enrichment analysis was performed on all differentially associated genes using the negative binomial distribution method.

2. The analysis method according to claim 1, characterized in that In step S1, the unqualified sequences include: ① reads with adapters; ② reads containing more than 5% N bases; ③ reads with a quality value Q≤10 that accounts for more than 20% of the total reads; and ④ reads with all bases being A or N.

3. The analysis method according to claim 1, characterized in that In step S4, the Peak identification method is: S41, scan the exon region to obtain the candidate Peak region; S42, assuming the gene is located in the exon window w = [w1, w2], calculate the number of reads that the IP and Input map to the window w = [w1, w2]: X 1,w ~Poisson(n1λ 1,w ), X 0,w ~Poisson(n0λ 0,w ), Among them, X 1,w and X 0,w is the number of reads of the IP sample and the Input sample aligned to the window w=[w1,w2], n1 and n0 are the total number of reads of the IP sample and the Input sample, respectively, 1,w and λ 0,w Represent the probability that the sequencing reads of IP samples and Input samples fall into window w respectively; S43, use C-test to analyze the X of the candidate Peak region 1,w ~Poisson(n1λ 1,w ) and X 0,w ~Poisson(n0λ 0,w ) ratio c=λ 1,w / λ 0,w Make comparisons and calculate p-value; S44, use Fisher test to calculate the p-value of the candidate Peak region. The region smaller than the p-value is the Peak, where p-value < 0.

05.

4. The analysis method according to claim 1, characterized in that In step S7, the FPKM value calculation formula is as follows: Among them, total exon fragments represents the total number of exon fragment reads, It indicates the number of millions of reads that can be mapped to the genome, and exon length indicates the kb length of the exon.

5. The analysis method according to claim 1, characterized in that In step S8, the negative binomial distribution method is used to perform GO and KEGG functional enrichment analysis, and the formula is: ， Among them, pvalue of enriched GO / KEGG represents the GO / KEGG functional enrichment analysis value, TB represents all background genes, TS represents all differentially expressed genes, B represents all genes with GO or KEGG annotations, and S represents all differentially expressed genes with GO or KEGG annotations.

6. An analysis system for m6A-seq sequencing data, characterized in that: Includes the following modules: Data input module, used to accept m6A-seq high-throughput sequencing data input and perform preprocessing; The alignment module and the data input module are used to align the preprocessed sequence with the reference genome of the sample source species to obtain reads that can be aligned to the reference genome; The Peak Calling and Annotation module is connected to the alignment module and is used to perform genome-wide peak identification on reads that can be aligned to the reference genome, obtain peak positions on the genome, peak lengths, and read quantity information, annotate peak positions on the genome, and analyze the distribution of reads at transcription start and end positions; The differential peak analysis module is connected to the Peak Calling and Annotation module to perform differential peak analysis on the peaks at the intersection of the IP sample and the Input sample to obtain differential peaks, and to perform GO and KEGG functional enrichment analysis on the related genes corresponding to the differential peaks using the negative binomial distribution method; Peak sequence motif analysis module, connected to the difference peak analysis module, for performing motif analysis on the obtained single sample peaks and peaks with differences; The mRNA and lncRNA expression abundance analysis module is connected to the comparison module and is used to quantitatively analyze the mRNA and lncRNA in the sample to obtain the original value Counts and normalized expression quantity FPKM value of each mRNA and lncRNA in the sample; The differential mRNA and lncRNA expression profile analysis module is connected to the mRNA and lncRNA expression abundance analysis module to perform differential analysis on each mRNA and lncRNA in the IP sample and the Input sample, and perform GO and KEGG functional enrichment analysis on all differentially related genes using the negative binomial distribution method.

7. The system according to claim 6, characterized in that Preprocessing in the data input module refers to filtering out unqualified sequences and / or contaminating sequences. The unqualified sequences include: ① reads with adapters; ② reads containing more than 5% N bases; ③ reads with a quality value Q≤10 that accounts for more than 20% of the total reads; and ④ reads with all bases being A or N. Contaminating sequences are filtered out using a database including ribosomal RNA sequences and / or mycoplasma RNA sequences.

8. The system according to claim 6, wherein: In the Peak Calling and Annotation module, the Peak identification method is: S41, scan the exon region to obtain the candidate Peak region; S42, assuming the gene is located in the exon window w = [w1, w2], calculate the number of reads that the IP and Input map to the window w = [w1, w2]: X 1,w ~Poisson(n1λ 1,w ), X 0,w ~Poisson(n0λ 0,w ), Among them, X 1,w and X 0,w is the number of reads of the IP sample and the Input sample aligned to the window w=[w1,w2], n1 and n0 are the total number of reads of the IP sample and the Input sample, respectively, 1,w and λ 0,w Represent the probability that the sequencing reads of IP samples and Input samples fall into window w respectively; S43, use C-test to analyze the X of the candidate Peak region 1,w ~Poisson(n1λ 1,w ) and X 0,w ~Poisson(n0λ 0,w ) ratio c=λ 1,w / λ 0,w Make comparisons and calculate p-value; S44, use Fisher test to calculate the p-value of the candidate Peak region. The region smaller than the p-value is the Peak, where p-value < 0.

05.

9. The system according to claim 6, wherein: In the mRNA and lncRNA expression abundance analysis module, the FPKM value calculation formula is as follows: Among them, total exon fragments represents the total number of exon fragment reads, It indicates the number of millions of reads that can be mapped to the genome, and exon length indicates the kb length of the exon.

10. The system according to claim 6, wherein: In the differential mRNA and lncRNA expression profile analysis module, the negative binomial distribution method is used to perform GO and KEGG functional enrichment analysis, and the formula is: ， Among them, pvalue of enriched GO / KEGG represents the GO / KEGG functional enrichment analysis value, TB represents all background genes, TS represents all differentially expressed genes, B represents all genes with GO or KEGG annotations, and S represents all differentially expressed genes with GO or KEGG annotations.

Images