Internuclear allele expression analysis method and device and storage medium

Abstract

The invention discloses a binuclear internuclear allele expression analysis method and device and a storage medium, and relates to the technical field of bioinformatics transcriptome analysis of binuclear fungi. The method comprises the following steps: respectively carrying out comparative analysis on sequencing data of binuclear fungi and two nuclear reference genomes; performing kernel source classification on each read in the sequencing data by using an SVM model; taking all gene sequences of the kernel 1 as a query set, taking a gene sequence of the kernel 2 as a target set, and constructing a similarity score matrix among the gene sequences by using a BLAST tool; combining the similarity score with gene annotation information, and analyzing and establishing an internuclear allele relationship through a LightGBM model; and generating a counts file, performing FPKM standardization processing, and analyzing the standardized expression quantity in combination with the internuclear allele relationship index to obtain an analysis result. According to the invention, the analysis of the interkaryotic alleles of the binuclear fungi is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics transcriptome analysis of dikaryotic fungi, and in particular to an inter-nuclear allele expression analysis method, device and storage medium. Background Art

[0002] Fungi are eukaryotic, spore-forming, chloroplast-free eukaryotic organisms, including molds, yeasts, mushrooms, and other mushrooms known to humans. As of 2024, about 150,000 fungal species have been discovered and recorded by the scientific community, including yeasts, molds, and various mushrooms. However, dikaryotic fungi have only been noticed in the past decade, such as the Basidiomycetes of fungi. Dikaryons can generate monokaryons through sexual and asexual methods, and compatible monokaryons form dikaryons with lock-like union through plasmogamy. The current research results on the monokaryons of edible fungi have expanded the understanding of the biological properties of monokaryons and dikaryons. As a haploid dikaryotic form, dikaryotic cells have many biological characteristics different from diploids. Many researchers have conducted some research on the hierarchical order, adaptability, stability and other biological properties of dikaryons. The relevant research results have brought many inspirations for fungal species maintenance, hybrid breeding innovation, and the construction of genetic research paradigms. Analysis of the incrementality, fertility, and totipotency of dikaryons as fungi can provide inspiration for solving the problem of maintaining the long-term stability of edible fungi species. The dominant and recessive relationships and complementary relationships between dikaryon nuclei may have important guiding significance for hybrid breeding innovation. Internuclear alleles have different expression levels in dikaryon cells, providing new ideas for studying scientific issues of edible fungi from the perspective of dikaryons. Research on dikaryons of edible fungi may form a more disciplinary research paradigm for edible fungi.

[0003] At present, in the research progress of dikaryon fungi, researchers have proved that there is a certain order relationship between the two nuclei in dikaryon fungi, which affects the phenomenon of partial separation of cell development and the expression mechanism of inter-nuclear alleles in dikaryon fungi, that is, the alleles existing on the two sets of chromosomes in dikaryon cells show different expression patterns, and one of the nuclei may be the dominant nucleus. These studies have confirmed the existence of such biological laws in some papers, but there are still deficiencies in the study of dikaryons. The traditional transcriptome analysis method is to regard dikaryon fungi as monokaryons, and use the reference genome of one of the nuclei or the reference genome corresponding to the fungus on NCBI to compare its reads, ignoring the contribution of the two nuclei in gene expression and their order relationship between dikaryon fungi, resulting in incomplete analysis of dikaryon fungi and the inability to discover the biological laws behind the binucleus. In order to more completely analyze the transcriptome of dikaryon fungi and explore the biological laws between the binuclei of dikaryon fungi, the present invention provides a unique bioinformatics analysis method for dikaryon fungi. Summary of the invention

[0004] The purpose of the present invention is to provide a method, device and storage medium for analyzing internuclear allele expression, so as to study the contribution of each nucleus to gene expression and biological laws such as order, and to solve the problems of incomplete existing analysis and the lack of internuclear allele analysis methods of binuclear fungi mentioned in the above background technology.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] According to one aspect of the present invention, a method for analyzing inter-nuclear allele expression is provided, comprising: obtaining sequencing data of a dikaryon fungus and performing quality control on the sequencing data; comparing the sequencing data after quality control with two nuclear reference genomes of the dikaryon fungus to generate two sam files, wherein the two nuclear reference genomes of the dikaryon fungus are respectively the nucleus 1 reference genome and the nucleus 2 reference genome, and the two sam files respectively record the comparison information of each read in the sequencing data on the two nuclear reference genomes; fusing the two sam files obtained after the comparison to generate a summary sam file containing all the comparison information; using a trained SVM model, according to the comparison information in the summary sam file, classifying the nuclear origin of each read in the sequencing data, classifying the reads from nucleus 1 into the nucleus1.sam file, and the reads from nucleus 2 into the nucleus2.sam file; according to the genome annotation files of the two nuclei of the dikaryon fungus, extracting the unique ID, chromosome name, start position and end position of each gene; according to the position information provided in the genome annotation file , extract the sequences of corresponding genes from the nuclear 1 reference genome and the nuclear 2 reference genome respectively; use the gene sequences of all genes in the nuclear 1 genome annotation file as the query set, and the gene sequences of all genes in the nuclear 2 genome annotation file as the target set, use the BLAST tool to perform pairwise sequence comparisons on the genes in the nuclear 1 genome annotation file and the genes in the nuclear 2 genome annotation file and their corresponding base sequences, obtain the similarity scores between the gene sequences, and construct a score matrix; combine the obtained similarity scores with the gene annotation information, and establish the inter-nuclear allele relationship between the two nuclei through LightGBM model analysis; convert the nucleus1.sam and nucleus2.sam files after nuclear source classification into computer-readable nucleus1.bam and nucleus2.bam files; use the genome annotation files corresponding to the nuclei to perform quantitative expression analysis on nucleus1.bam and nucleus2.bam respectively, obtain the expression information of each gene, and generate counts files; perform FPKM normalization on the generated counts files, and convert nucleus1 and nucleusl2 The FPKM-normalized expression was combined with the internuclear allele relationship index to analyze the standardized expression data of genes in core 1 and their internuclear alleles in core 2, to compare the expression differences of the two internuclear alleles, or to judge the nuclear order relationship in binuclear fungi based on the standardized expression of the internuclear alleles, and to obtain the order relationship of different genes expressed between nuclei.

[0007] Optionally, the quality control of the sequencing data includes: performing quality control on the obtained binuclear fungus sequencing data by bioinformatics quality control technology; filtering reads in the sequencing data whose base quality scores are lower than a preset quality score threshold and whose read lengths are lower than a preset length threshold, wherein the base quality scores are calculated by the formula Q=-10·log 10 P is calculated, where Q represents the base quality score and P represents the probability of the base being misidentified.

[0008] Optionally, the alignment information of each read in the sequencing data on the two nuclear reference genomes includes: read alignment position, alignment score, alignment quality, number of base mismatches, read coverage length, and insert length.

[0009] Optionally, the sequencing data after quality control is compared and analyzed with the two nuclear reference genomes of the dikaryotic fungus respectively to generate two sam files, and the two sam files obtained after the comparison are merged to generate a summary sam file containing all the comparison information, including: comparing the sequencing data after quality control with the nuclear 1 reference genome to obtain the comparison file 1 sam file; comparing the sequencing data after quality control with the nuclear 2 reference genome to obtain the comparison file 2 sam file; extracting the comparison information fields from the two sam files of comparison file 1 and comparison file 2, and putting the extracted field information into a summary sam file, wherein the comparison information fields include required fields and optional fields The required fields include QNAME, RNAME, POS, MAPQ, CIGAR, MRNM and PNEXT, and the optional fields are other fields that are helpful for SVM model classification except the required fields, wherein QNAME represents the name of the reads, RNAME represents the name of the reference sequence to be aligned, POS represents the starting position of the alignment to the reference sequence, MAPQ represents the alignment quality score, reflecting the credibility of the alignment position, CIGAR represents the compression information of the alignment, indicating how the query sequence is aligned with the reference sequence, MRNM represents the reference sequence name of the partner read in the paired alignment, and PNEXT represents the starting position of the partner read in the paired alignment on the reference genome.

[0010] Optionally, before classifying the nuclear origin of each read in the sequencing data using the trained SVM model according to the alignment information in the summary sam file, the method further includes training the SVM model, and by classifying the alignment information of the reads in the sequencing data, the SVM model is trained to determine whether each read is derived from nucleus 1 or nucleus 2, thereby achieving nuclear origin classification, wherein the training process includes: feature standardization: normalizing the input features to ensure that the numerical ranges of different features are consistent; model construction: using the SVM algorithm to construct a classification model, the input features are RNAME1, RNAME2, Δ, |MAPQ1-MAPQ2|, and the classification rules are:

[0011]

[0012] Among them, Category indicates the classification result of reads, 1 represents reads from nucleus 1, 2 represents reads from nucleus 2, and Undecided represents no clear classification; RNAME1 represents the sequence name of reads aligned on reference genome 1; RNAME2 represents the sequence name of reads aligned on reference genome 2;

[0013] MAPQ1 and MAPQ2 represent the alignment quality of reads to reference genomes 1 and 2, respectively; Δ represents the alignment score difference, Δ=AS1-AS2, where AS1 is the alignment score to reference genome 1 and AS2 is the alignment score to reference genome 2; T Δ Indicates the significance threshold of the comparison score difference, which is used to determine whether the score has a significant difference; T MAPQ The significance threshold of the alignment quality difference is used to determine whether the alignment quality has a significant difference; parameter optimization: cross-validation is used to determine the kernel function and hyperparameters of the SVM model to optimize the classification performance of the model; training model: the SVM model is trained with training data labeled Category to obtain an SVM model that classifies reads by nuclear origin based on the input features.

[0014] Optionally, according to the alignment information in the summary sam file, each read in the sequencing data is classified by nuclear origin, and the reads from nucleus 1 are classified into the nucleus1.sam file, and the reads from nucleus 2 are classified into the nucleus2.sam file, including: according to the alignment information of the reads on the two nuclear reference genomes, the confidence value of each read belonging to nucleus 1 or nucleus 2 is calculated; when the confidence value of the reads belonging to nucleus 1 is greater than that of nucleus 2 and higher than a certain threshold, the reads will be classified into nucleus 1 and written into the nucleus1.sam file; when the confidence value of the reads belonging to nucleus 2 is greater than that of nucleus 1 and higher than a certain threshold, the reads will be classified into nucleus 2 and written into the nucleus2.sam file; when the confidence of the reads belonging to both nucleus 1 and nucleus 2 is lower than a preset threshold, the reads will be discarded.

[0015] Optionally, based on the alignment information in the summary sam file, each read in the sequencing data is classified by nuclear origin, and the reads from nucleus 1 are classified into the nucleus1.sam file, and the reads from nucleus 2 are classified into the nucleus2.sam file, which also includes: if the reads can only find a unique alignment position on one reference genome, the classification result is marked according to the reference genome, specifically: if the reads are uniquely aligned on the nucleus 1 reference genome, the reads will be classified into nucleus 1 and written into the nucleus1.sam file; if the reads are uniquely aligned on the nucleus 2 reference genome, the reads will be classified into nucleus 2 and written into the nucleus2.sam file.

[0016] Optionally, after combining the obtained similarity score with the gene annotation information and establishing the inter-nuclear allele relationship between the two nuclei through LightGBM model analysis, the method further includes: for genes for which the inter-nuclear allele relationship cannot be established through the LightGBM model, associating them with NA, wherein NA is used to indicate the absence of inter-nuclear alleles.

[0017] According to another aspect of the present invention, there is also provided an inter-nuclear allele expression analysis device, comprising a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the computer program implements the steps of the method described above when executed by the processor.

[0018] According to another aspect of the present invention, a storage medium is provided, on which a computer program is stored. When the computer program is executed by a processor, the steps of the method described above are implemented.

[0019] Compared with the traditional single-nuclear analysis method, the present invention improves the integrity of binuclear fungal analysis, and processes and classifies data through two machine learning methods, the SVM (support vector machine) model and the LightGBM model, making the analysis results more accurate and reliable. The method of the present invention can clearly show the expression mode of inter-nuclear alleles of binuclear fungi, and explore more biological laws between inter-nuclear allele groups. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] The drawings described herein are used to provide a further understanding of the present invention and constitute a part of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0021] Figure 1 is a schematic flow chart of a method for analyzing internuclear allele expression in Example 1 of the present invention;

[0022] Figure 2 is a schematic flow chart of a method for analyzing internuclear allele expression in Example 2 of the present invention;

[0023] Figure 3 Schematic diagram of the percentage of genes for which allelic relationships have been established in Example 2 of the present invention;

[0024] Figure 4 This is a Venn diagram of the number of genes whose gene-nuclear order relationship has been determined in Example 2 of the present invention. DETAILED DESCRIPTION

[0025] In order to enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only embodiments of a part of the present application, not all embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in the field without creative work should fall within the scope of protection of the present application. It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict.

[0026] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The term "in one embodiment" that appears in different places in this specification does not necessarily refer to the same embodiment, nor does it refer to a separate or selective embodiment that is mutually exclusive with other embodiments.

[0027] It should be noted that the terms "first", "second", etc. in the specification and claims of the present application and the above-mentioned drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequence. It should be understood that the data used in this way can be interchangeable where appropriate, so that the embodiments of the present application described herein can be implemented in an order other than those illustrated or described herein. In addition, the terms "including" and "having" and any of their variations are intended to cover non-exclusive inclusions, for example, a process, method, product or device comprising a series of steps or units is not necessarily limited to those steps or units clearly listed, but may include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0028] Example 1

[0029] In this embodiment, a method for analyzing internuclear allele expression is provided.

[0030] Reference Figure 1 , Figure 1 is a schematic flow chart of a method for analyzing internuclear allele expression in an embodiment of the present invention, such as Figure 1 As shown, the method comprises the following steps:

[0031] S1, obtaining sequencing data of dikaryon fungi, and performing quality control on the sequencing data;

[0032] The genome of dikaryon fungi is sequenced by high-throughput sequencing technology to obtain sequencing data of dikaryon fungi. These data are usually stored in FASTQ format and contain original sequencing reads.

[0033] Then, the obtained binuclear fungal sequencing data were quality controlled through bioinformatics quality control technology, such as fastp, fastqC and other bioinformatics software, and low-quality reads were filtered to obtain high-quality sequencing data to be analyzed.

[0034] Specific quality control filtering strategies include: according to the formula Q = -10·log 10P calculates the base quality score, where P represents the probability of the base being misidentified; reads with base quality scores lower than a preset quality score threshold are filtered. This threshold is usually set according to the data quality requirements. Common quality score thresholds are Q20 (i.e., the error rate is 1%) or Q30 (i.e., the error rate is 0.1%). If the quality scores of some bases in the reads are lower than the set quality score threshold, then these reads will be considered low quality and filtered out. In addition, shorter reads contain less information and are more susceptible to sequencing errors. Therefore, in the embodiment of the present invention, reads with a read length lower than a preset length threshold are also filtered. The read length threshold can be set according to specific analysis requirements, and the embodiment of the present invention does not make specific restrictions.

[0035] The remaining reads after quality control will be used as high-quality sequencing data for subsequent internuclear allele expression analysis.

[0036] S2, the sequencing data after quality control were compared with the two nuclear reference genomes of dikaryotic fungi to generate two sam files;

[0037] The two nuclear reference genomes of dikaryotic fungi are nucleus 1 reference genome (nucleus1) and nucleus 2 reference genome (nucleus2). These reference genomes are known genome sequences and are used for comparison with sequencing data.

[0038] The high-quality sequencing data after quality control were aligned with the reference genomes of the two nuclei of the dikaryotic fungi using bioinformatics alignment tools such as hisat2. The alignment process will generate two sam (SequenceAlignment / Map) files, each of which corresponds to the alignment results of a reference genome. Among them, the alignment file 1.sam file contains the alignment information of the sequencing data with the reference genome of nucleus 1, and the alignment file 2.sam file contains the alignment information of the sequencing data with the reference genome of nucleus 2. The alignment information includes reads alignment position, alignment score, alignment quality, number of base mismatches, reads coverage length, insert fragment length, etc. Among them, the alignment position is the specific position of the sequencing reads on the reference genome; the alignment score is the match degree score between the reads and the reference genome; the alignment quality indicates the reliability or accuracy of the alignment result; the number of base mismatches refers to the number of bases that do not match between the reads and the reference genome; the reads coverage length refers to the length of the continuous region where the reads are successfully aligned on the reference genome; the insert fragment length refers to the length of the fragment in the reads that is longer than the reference genome.

[0039] S3, merge the two sam files obtained after comparison to generate a summary sam file containing all comparison information;

[0040] Find and extract the necessary alignment information fields in the two sam files, including but not limited to the required fields such as QNAME, RNAME, POS, MAPQ, CIGAR, MRNM, PNEXT and optional fields that are helpful for SVM classification. The optional fields can be customized according to actual conditions or needs. Merge the field information extracted from the two sam files into a new sam file to obtain a summary sam file containing all alignment information, where QNAME represents the name of the reads (query name), which is used to uniquely identify each sequencing read; RNAME represents the name of the reference sequence to which the reads are aligned (reference name); POS represents the starting position (position) of the reads aligned to the reference sequence; MAPQ represents the mapping quality score, which reflects the credibility of the alignment position; CIGAR represents the compression information of the alignment, indicating how the query sequence is aligned with the reference sequence; MRNM represents the reference sequence name of the partner read in the pairwise alignment, and PNEXT represents the starting position of the partner read in the pairwise alignment on the reference genome (1 is the starting point).

[0041] By extracting required and optional fields, we filter out field information irrelevant to classification, such as ZS (RNA chain information during alignment), YS (score of suboptimal alignment), YT (alignment type information), CL (command line information for recording SAM file generation), and other redundant fields. These fields usually do not contain information useful for SVM classification and may increase the complexity of data processing. This reduces the redundancy of useless information and the interference of useless information on the model under millions of reads, thereby improving the accuracy of the model.

[0042] The merged sam file retains the alignment information of all reads, such as QNAME, RNAME, POS, MAPQ, CIGAR, MRNM, PNEXT, etc. This information helps to further analyze the alignment quality, sequencing accuracy and genome coverage, and provides necessary support for subsequent feature extraction and classification model construction.

[0043] S4, using the trained SVM model, classify the nuclear origin of each read in the sequencing data according to the alignment information in the summary sam file, and classify the reads from nucleus 1 into the nucleus1.sam file, and the reads from nucleus 2 into the nucleus2.sam file;

[0044] Because the two cells of the dikaryotic fungus The nucleus comes from different parents, and the bases in the gene sequence will accumulate gene base mutations during the generation and development process, resulting in base differences in the gene sequences of the two nuclei. Based on the principle that the same genes on the two reference genomes still have base differences in sequence, the SVM (Support Vector Machine) model trained with the fungal dataset is used to classify the nuclear origin of each read based on the required field information such as QNAME, RNAME, POS, MAPQ, CIGAR, MRNM, PNEXT of the reads in the summary sam file and the optional field information that is helpful for SVM classification (determined according to different alignment situations), and a confidence value is given to each nuclear source. When the confidence value of the reads in nucleus 1 is greater than that of nucleus 2 and higher than a certain threshold, the reads will be classified into nucleus 1 and written into the nucleus1.sam file. When the confidence value of the reads in nucleus 2 is greater than that of nucleus 1 and higher than a certain threshold, the reads will be classified into nucleus 2 and written into the nucleus2.sam file. When the confidence of the reads in nucleus 1 and nucleus 2 is lower than the preset threshold, the reads will be discarded.

[0045] In addition, if a read can only find a unique alignment position on one reference genome, that is, the alignment result is unique on the reference genome, but there is no alignment result on the other reference genome, then the classification result can be directly based on the reference genome marker. Specifically, if the read is uniquely aligned on the nucleus 1 reference genome, the read will be classified into nucleus 1 and written into the nucleus1.sam file; if the read is uniquely aligned on the nucleus 2 reference genome, the read will be classified into nucleus 2 and written into the nucleus2.sam file.

[0046] Reads from nucleus1 will be classified into one sam file, and reads from nucleus2 will be classified into another sam file. Reads that cannot be classified will be filtered.

[0047] Furthermore, before classifying the nuclear source of each read in the sequencing data using the trained SVM model according to the alignment information in the summary sam file, the method further includes training the SVM model, and by classifying the alignment information of the reads in the sequencing data (including alignment position, alignment score, alignment quality and other features), the SVM model is trained to determine whether each read is derived from core 1 or core 2, thereby achieving nuclear source classification, wherein the training process includes:

[0048] (1) Feature normalization: Normalize the input features to ensure that the numerical ranges of different features are consistent. For example, Δ and |MAPQ1-MAPQ2| are normalized according to the z-score.

[0049] (2) Model construction: Use the SVM algorithm to build a classification model. The input features are RNAME1, RNAME2, Δ, |MAPQ1-MAPQ2|. Define classification rules, including basic classification rules and complex classification rules.

[0050] Basic classification rules: If a read can only be aligned on one genome, the classification result is directly marked as 1 or 2. In this case, the reads are usually highly reliable and can be used as clear samples in the training set;

[0051] Complex classification rules: If a read is aligned to two genomes at the same time, that is, there are alignment results on both genomes, the alignment score difference Δ and alignment quality difference are used

[0052] |MAPQ1-MAPQ2| conducts further determination, specifically:

[0053] When Δ>T Δ , indicating that the score of the genome 1 is significantly higher, and the classification result is 1, indicating that it comes from nucleus 1;

[0054] When Δ<-T Δ , indicating that the score of genome 2 is significantly higher, and the classification result is 2, indicating that it comes from nucleus 2;

[0055] When |Δ|≤T Δ And |MAPQ1-MAPQ2|≤T MAPQ When , the classification result is Undecided, indicating uncertain classification; for the Undecided category, you can design a separate rule to exclude it, or add a third category as a classification output of SVM;

[0056] The formula of the above classification rule is as follows:

[0057]

[0058] Among them, Category indicates the classification result of reads, 1 represents reads from nucleus 1, 2 represents reads from nucleus 2, and Undecided represents no clear classification; RNAME1 represents the sequence name of reads aligned on reference genome 1; RNAME2 represents the sequence name of reads aligned on reference genome 2;

[0059] MAPQ1 and MAPQ2 represent the alignment quality of reads to reference genomes 1 and 2, respectively; Δ represents the alignment score difference, which is used to measure the relative tightness of the alignment of reads between the two nuclei, Δ = AS1-AS2, where AS1 is the alignment score to reference genome 1,

[0060] AS2 is the alignment score to reference genome 2; T Δ Indicates the significance threshold of the comparison score difference, which is used to determine whether the score has a significant difference; T MAPQ The significance threshold of the comparison quality difference is used to determine whether the comparison quality has a significant difference;

[0061] (3) Parameter optimization: Use cross-validation to determine the kernel function (such as linear kernel, RBF kernel, etc.) and hyperparameters (such as C value and gamma value) of the SVM model to optimize the classification performance of the model;

[0062] (4) Training model: Use the training data labeled Category to train the SVM model to obtain an SVM model that can classify the nuclear origin of reads based on the input features.

[0063] S5, extract the unique ID, chromosome name, start position, and end position of each gene based on the genome annotation files of the two nuclei of the dikaryotic fungus;

[0064] The genome annotation file is a file that interprets and annotates the genome sequence and contains detailed annotation information of the genome, such as the location, structure, and function of the gene.

[0065] S6, extracting the sequences of corresponding genes from the core 1 reference genome and the core 2 reference genome respectively according to the position information provided in the genome annotation file;

[0066] When constructing the allelic relationship of genes between the two nuclei of dikaryotic fungi through the genome annotation files of the two nuclei, by recording the starting and ending positions on the corresponding chromosomes and ensuring that each gene has a unique ID, the reference genome file (FASTA format) is used to locate the chromosome name, starting and ending positions of the genes in the annotation file to the reference genome of the corresponding nucleus, and the sequences of the corresponding genes are extracted from the reference genomes of nuclei 1 and 2, respectively.

[0067] S7, taking the gene sequences of all genes in the core 1 genome annotation file as the query set and the gene sequences of all genes in the core 2 genome annotation file as the target set, using the BLAST tool to perform pairwise sequence comparisons of the genes in the core 1 genome annotation file and the genes in the core 2 genome annotation file and their corresponding base sequences, obtaining similarity scores between the gene sequences, and constructing a score matrix;

[0068] BLAST (Basic Local Alignment Search Tool) is a commonly used tool software in bioinformatics. It can compare the input nucleic acid or protein sequence with the known sequences in the database to obtain information such as sequence similarity, thereby determining the source or evolutionary relationship of the sequence.

[0069] The gene sequences of all genes in the core 1 genome annotation file are used as the query sequence set, and the gene sequences of all genes in the core 2 genome annotation file are used as the target sequence set. BLAST is performed on each pair of gene sequences to obtain the similarity scores (such as alignment scores) between the gene sequences. Finally, a score matrix is ​​obtained with all the gene IDs of core 1 as row names and all the gene IDs of core 2 as column names. Each element in the score matrix represents the similarity score between a gene in core 1 and a gene in core 2.

[0070] S8, the obtained similarity scores were combined with the gene annotation information, and the internuclear allele relationship between the two nuclei was established through LightGBM model analysis;

[0071] Two genetically different cell nuclei carry different allelic variants of the same gene, and the relationship between them is an allelic relationship.

[0072] LightGBM (Light Gradient Boosting Machine) is a distributed gradient boosting model based on a decision tree algorithm. This step uses the LightGBM model to further analyze and infer the allele relationship between these unrelated genes. When the difference between the BLAST similarity and the similarity of all genes in another nucleus is lower than a certain threshold, it is judged that the inter-nuclear allele cannot be established. The genes that still cannot be analyzed will be established with NA, indicating that there is no inter-nuclear allele, and it is considered that the gene is regulated by only one nucleus.

[0073] S9, converting the nucleus1.sam and nucleus2.sam files after nuclear source classification into computer-readable nucleus1.bam and nucleus2.bam files;

[0074] sam files can be converted into computer-readable bam files through bioinformatics tools, such as Samtools, for subsequent analysis.

[0075] S10, using the genome annotation file corresponding to the nucleus to perform quantitative expression analysis on nucleus1.bam and nucleus2.bam respectively, obtain the expression information of each gene, and generate a counts file;

[0076] Quantitative expression analysis is to compare reads and genome annotation files, and count the number of reads of each gene to reflect the transcription level of the gene in nucleus 1 and nucleus 2. Through tools for quantitative expression analysis of bam files, such as featureCounts, the nucleus1.bam and nucelus2.bam files are respectively compared with the nucleus1.gff (gtf) and nucleus2.gff (gtf) genome annotations for quantitative expression analysis, and the counts data of each nucleus are analyzed.

[0077] S11, perform FPKM standardization on the generated counts file, combine the FPKM-standardized expression levels of nucleus1 and nucleus2 with the inter-nuclear allele relationship index, analyze the standardized expression data of the genes in nucleus 1 and their inter-nuclear alleles in nucleus 2, compare the expression differences of the two inter-nuclear alleles, or judge the nuclear order relationship in binuclear fungi based on the standardized expression levels of the inter-nuclear alleles, and obtain the order relationship of different genes expressed between nuclei.

[0078] The quantitative files of nucleus1 and nucleus2 were normalized by FPKM (Fragments Per Kilobase of transcript per Million mapped reads) data to eliminate the quantitative differences between reads during sequencing and the influence of the library. The expression levels of nucleus1 and nucleus2 after FPKM normalization were combined with the corresponding internuclear allele relationship index. By comparing the standardized expression data of genes in nucleus 1 and their internuclear alleles in nucleus 2, combined with statistical methods such as paired t-test, the expression differences of the two internuclear alleles can be analyzed and compared.

[0079] The nuclear order relationship in dikaryotic fungi is determined based on the difference in expression of alleles between two nuclei, or based on the standardized expression of alleles between nuclei. Based on the analysis results, the order relationship of expression of different genes between nuclei is obtained, such as which genes are more highly expressed in one nucleus and which are more highly expressed in the other nucleus, and the "dominant" or "subordinate" status of the nucleus is inferred accordingly.

[0080] The method for analyzing the inter-nuclear allele expression of binary fungi of the present invention improves the integrity of the binary fungi analysis, and processes and classifies data through two machine learning methods, the SVM (support vector machine) model and the LightGBM model, making the analysis results more accurate and reliable. The method of the present invention can clearly show the expression mode of the inter-nuclear alleles of binary fungi, and explore more biological laws between the inter-nuclear allele groups.

[0081] Example 2

[0082] The technical solution of the present invention is further described in detail below in conjunction with a specific embodiment and the accompanying drawings. It should be understood that the following embodiment is only used to explain the present invention, and is not used to limit the present invention.

[0083] The implementation scenario in this embodiment is black fungus, and the internuclear allele expression analysis method of the present invention is used to perform binuclear internuclear allele transcriptome analysis on black fungus.

[0084] Black fungus (scientific name Auricularia heimuer), also known as cloud ear, fungus, wood ear, light fungus, wood core, wood fungus, tree chicken, It is an edible fungus of the genus Auricularia in the family Auriculariaceae. Auricularia is a binuclear fungus. The life history of Auricularia includes basidiospores germinating into monokaryotic hyphae or forming sickle-shaped conidia, and then germinating from conidia into monokaryotic hyphae. Monokaryotic hyphae of different sexes combine to form binucleate hyphae, and continue to proliferate through lock-like union. The binucleate hyphae that have reached the physiological maturity stage form fruiting body primordia on the surface of the substrate and develop into fruiting bodies.

[0085] This example uses reference genomes with a nucleus1 genome size of about 41.12Mb, a GC content of 52.68%, an N50 of 1168003bp and a nucleus2 genome size of about 40.52Mb, a GC content of 54.32%, and an N50 of 1302490bp. Among them, genome size refers to the DNA content in a genome; GC content refers to the proportion of guanine and cytosine in the genome; N50 refers to the length of the shortest fragment that exceeds half (50%) of the total length of all fragments after all assembled fragments are sorted from long to short.

[0086] like Figure 2 As shown, the method for analyzing internuclear allele expression in the embodiment of the present invention comprises the following steps:

[0087] Step 1: fastp (0.23.4) was used to perform quality control analysis on the raw sequencing data. The high-quality sequencing data after quality control were aligned to the reference genomes of the two nuclei (nucleus1 and nucleus2) using Hisat2 (version 2.2.1). The global alignment rates of nucleus1 and nucleus2 were 97.72% and 92.53% respectively.

[0088] Step 2: According to the above aligned data, the unaligned data were filtered and the valid reads were retained. The total reads after quality control by the fastp tool were 8.23M. After filtering, the remaining reads in nucleus1 were 8.035M, and the remaining reads in nucleus2 were 7.703M.

[0089] Step 3, merge the obtained alignment file 1.sam with the alignment file 2.sam, and remove the redundant field information to obtain a sam file containing 15737980 reads alignment information;

[0090] Step 4: Use the trained SVM classification model to classify each read in the sam file after fusion in step 3, write the read information of the nucleus source of nucleus 1 into the nucleus1.sam file, and write the read information of the nucleus source of nucleus 2 into the nucleus2.sam file; the reads that cannot be classified are discarded; the number of reads successfully classified is 12401544, of which the number of reads of nucleus 1 and nucleus 2 are 5704710 and 6696834 respectively;

[0091] Step 5, use samtools view (1.9) to convert the sam file into a binary bam file;

[0092] Step 6, extracting the unique ID, chromosome name, start position and end position of each gene according to the genome annotation files of the two cores of the dikaryotic fungus, and performing gene sequence cutting on the reference genomes of core 1 and core 2 respectively;

[0093] Step 7: Use BLAST tool to compare the base sequences of genes in the genome annotation file from core 1 with those in the genome annotation file from core 2, and obtain a 15420*14903 scoring matrix, which is used to evaluate the similarity between gene sequences and construct a scoring matrix;

[0094] Step 8: Give the 15420*14903 scoring matrix score and the annotation information of each nuclear gene in its genome annotation file to the LightGBM machine learning model that has been preliminarily trained locally to establish the association between nuclear alleles, and discard the reads that cannot establish the association between nuclear alleles; refer to Figure 3 , Figure 3 Schematic diagram of the percentage of genes with established allele relationships between nuclei in this embodiment. Figure 3As shown, the proportion of genes with allelic relationships far exceeds that of genes without allelic relationships;

[0095] Step 9, measure the expression of nucleus1.bam and nucleus2.bam using featureCounts (v2.0.2), obtain nucleus1.counts and nucleus2.counts files, and obtain the expression information of each gene;

[0096] Step 10: According to the processing result of featureCounts, FPKM (Fragments Per Kilobase of transcript per Million mapped reads) is used to standardize the expression level of each gene to eliminate the differences caused by sequencing depth and read length. The standardization formula is: FPKM = (C / N) / (L / M)*10^6, where C = fragment count on the base; N = the total number of all valid fragments; L = the length of the target gene (in bases); M = the total length of all mapped fragments (in bases);

[0097] Step 11, based on the FPKM-normalized expression data, combined with the inter-nuclear allele relationship index, analyze the difference in allele expression between the two nuclei, and use the following statistical method to compare the inter-nuclear allele expression to obtain the analysis result:

[0098] |FPKMNucleus1-FPKMNucleus2|>A·Mean(FPKMNucleus1,FPKMNucleus2);

[0099] In the formula, FPKMNucleus1 represents the normalized expression level of the target gene in nucleus 1;

[0100] FPKMNucleus2 represents the standardized expression of the target gene in nucleus 2; A represents the set judgment threshold; Mean (FPKMN ucleus1 ,FPKMN ucleus2 ) represents the overall expression level of the target gene between the two nuclei; when A is 0.3, the judgment result is as follows Figure 4 As shown, Figure 4 The quantitative Venn diagram of the dominance of gene expression in nucleus 1 and nucleus 2 is displayed. Nucleus1Dominant means that the expression of genes in nucleus 1 is significantly higher than that of the gene set in nucleus 2. Nucleus2 Dominant means that the expression of genes in nucleus 2 is significantly higher than that of the gene set in nucleus 1. Undecided means that the expression difference of genes in nucleus 1 and nucleus 2 is not obvious, and the dominant gene set cannot be determined.

[0101] In this embodiment, a device for analyzing internuclear allele expression is also provided, which is used to implement the above embodiments and preferred implementation modes, and those that have been explained will not be repeated here.

[0102] Optionally, the specific examples in this embodiment may refer to the examples described in the above embodiments and optional implementation modes, and this embodiment will not be described in detail here.

[0103] An embodiment of the present invention further provides a storage medium, in which a computer program is stored, wherein the computer program is configured to execute the steps of any of the above method embodiments when running.

[0104] For specific examples in this embodiment, reference may be made to the examples described in the above embodiments and optional implementation modes, and this embodiment will not be described in detail here.

[0105] Optionally, in this embodiment, the above-mentioned storage medium may include but is not limited to: a USB flash drive, a read-only memory (ROM), a random access memory (RAM), a mobile hard disk, a magnetic disk or an optical disk, and other media that can store computer programs.

[0106] The serial numbers of the above embodiments of the present invention are only for description and do not represent the advantages or disadvantages of the embodiments.

[0107] In the above embodiments of the present invention, the description of each embodiment has its own emphasis. For parts that are not described in detail in a certain embodiment, reference can be made to the relevant descriptions of other embodiments.

[0108] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for analyzing internuclear allele expression, characterized in that: include: Acquiring sequencing data of dikaryotic fungi, and performing quality control on the sequencing data; The sequencing data after quality control are compared and analyzed with two nuclear reference genomes of the dikaryon fungus to generate two sam files, wherein the two nuclear reference genomes of the dikaryon fungus are respectively the nuclear 1 reference genome and the nuclear 2 reference genome, and the two sam files respectively record the comparison information of each read in the sequencing data on the two nuclear reference genomes; The two sam files obtained after comparison are merged to generate a summary sam file containing all comparison information; Using the trained SVM model, each read in the sequencing data is classified by nuclear origin according to the alignment information in the summary sam file, and the reads from nucleus 1 are classified into the nucleus1.sam file, and the reads from nucleus 2 are classified into the nucleus2.sam file; Based on the genome annotation files of the two nuclei of dikaryotic fungi, the unique ID, chromosome name, start position, and end position of each gene were extracted; According to the position information provided in the genome annotation file, the sequences of the corresponding genes were extracted from the core 1 reference genome and the core 2 reference genome respectively; The gene sequences of all genes in the core 1 genome annotation file are used as the query set, and the gene sequences of all genes in the core 2 genome annotation file are used as the target set. The BLAST tool is used to perform pairwise sequence comparisons on the genes in the core 1 genome annotation file and the genes in the core 2 genome annotation file and their corresponding base sequences, and the similarity scores between the gene sequences are obtained to construct a score matrix. The obtained similarity scores were combined with the gene annotation information, and the internuclear allele relationship between the two nuclei was established through LightGBM model analysis; Convert the nucleus1.sam and nucleus2.sam files after nuclear source classification into computer-readable nucleus1.bam and nucleus2.bam files; Use the genome annotation files corresponding to the nucleus to perform quantitative expression analysis on nucleus1.bam and nucleus2.bam, obtain the expression information of each gene, and generate counts files; The generated counts file was subjected to FPKM standardization. The FPKM-normalized expression levels of nucleus1 and nucleus12 were combined with the internuclear allele relationship index to analyze the standardized expression data of genes in nucleus 1 and their internuclear alleles in nucleus 2. The expression differences of the two internuclear alleles were compared, or the nuclear order relationship in binuclear fungi was determined based on the standardized expression levels of the internuclear alleles, and the order relationship of different genes expressed between nuclei was obtained.

2. The method for analyzing internuclear allele expression according to claim 1, characterized in that: The quality control of the sequencing data comprises: The quality of the obtained binkaryon fungal sequencing data was controlled by bioinformatics quality control technology; The reads in the sequencing data whose base quality scores are lower than the preset quality score threshold and whose read lengths are lower than the preset length threshold are filtered, wherein the base quality scores are calculated by the formula Q = -10·log 10 P is calculated, where Q represents the base quality score and P represents the probability of the base being misidentified.

3. The method for analyzing internuclear allele expression according to claim 1, characterized in that: The alignment information of each read in the sequencing data on the two nuclear reference genomes includes: read alignment position, alignment score, alignment quality, number of base mismatches, read coverage length, and insert length.

4. The method for analyzing internuclear allele expression according to claim 1, characterized in that: The sequencing data after quality control were compared with the two nuclear reference genomes of dikaryotic fungi to generate two sam files. The two sam files obtained after comparison were merged to generate a summary sam file containing all the comparison information, including: The sequencing data after quality control are aligned with the nuclear 1 reference genome to obtain the alignment file 1sam file; The sequencing data after quality control are aligned with the nuclear 2 reference genome to obtain the alignment file 2sam file; Extract alignment information fields from alignment file 1 and alignment file 2, and put the extracted field information into a summarized sam file, wherein the alignment information fields include required fields and optional fields, the required fields include QNAME, RNAME, POS, MAPQ, CIGAR, MRNM and PNEXT, and the optional fields are other fields that are helpful for SVM model classification except the required fields, wherein QNAME represents the name of the reads, RNAME represents the name of the reference sequence for alignment, POS represents the starting position of the alignment to the reference sequence, MAPQ represents the alignment quality score, reflecting the credibility of the alignment position, CIGAR represents the compression information of the alignment, indicating how the query sequence is aligned with the reference sequence, MRNM represents the reference sequence name of the partner read in the pairwise alignment, and PNEXT represents the starting position of the partner read in the pairwise alignment on the reference genome.

5. The method for analyzing internuclear allele expression according to claim 1, characterized in that: Before classifying the nuclear source of each read in the sequencing data using the trained SVM model according to the alignment information in the summary sam file, the method further includes training the SVM model, and by classifying the alignment information of the reads in the sequencing data, the SVM model is trained to determine whether each read is derived from nucleus 1 or nucleus 2, thereby achieving nuclear source classification, wherein the training process includes: Feature normalization: Normalize the input features to ensure that the value ranges of different features are consistent; Model construction: Use the SVM algorithm to build a classification model, and the input features are RNAME1, RNAME2, Δ, |MAP Q1 -MAP Q2 |, the classification rules are: Among them, Category indicates the classification result of reads, 1 represents reads from nucleus 1, 2 represents reads from nucleus 2, and Undecided represents no clear classification; RNAME1 represents the sequence name of reads aligned on reference genome 1; RNAME2 represents the sequence name of reads aligned on reference genome 2; MAPQ1 and MAPQ2 represent the alignment quality of reads to reference genomes 1 and 2, respectively; Δ represents the alignment score difference, Δ=AS1-AS2, where AS1 is the alignment score to reference genome 1 and AS2 is the alignment score to reference genome 2; T Δ Indicates the significance threshold of the comparison score difference, which is used to determine whether the score has a significant difference; T MAPQ The significance threshold of the comparison quality difference is used to determine whether the comparison quality has a significant difference; Parameter optimization: Use cross-validation to determine the kernel function and hyperparameters of the SVM model to optimize the classification performance of the model; Training model: Use the training data labeled Category to train the SVM model to obtain an SVM model that classifies the nuclear origin of reads based on the input features.

6. The method for analyzing internuclear allele expression according to claim 1, characterized in that: According to the alignment information in the summary sam file, each read in the sequencing data is classified by nuclear source, and the reads from nucleus 1 are classified into the nucleus1.sam file, and the reads from nucleus 2 are classified into the nucleus2.sam file, including: Based on the alignment information of reads on the two nuclear reference genomes, the confidence value of each read belonging to core 1 or core 2 is calculated; When the confidence value that a read belongs to nucleus 1 is greater than that of nucleus 2 and is above a certain threshold, the read will be classified as nucleus 1 and written into the nucleus1.sam file; When the confidence value that the reads belong to nucleus 2 is greater than that of nucleus 1 and above a certain threshold, the reads will be classified as nucleus 2 and written into the nucleus2.sam file; When the confidence that a read belongs to both core 1 and core 2 is lower than the preset threshold, the read will be discarded.

7. The method for analyzing internuclear allele expression according to claim 1, characterized in that: According to the alignment information in the summary sam file, each read in the sequencing data is classified by nuclear source, and the reads from nucleus 1 are classified into the nucleus1.sam file, and the reads from nucleus 2 are classified into the nucleus2.sam file. If reads can only find a unique alignment position on a reference genome, the classification results are based on the reference genome markers, specifically: If the reads are uniquely aligned to the nucleus 1 reference genome, the reads will be grouped into nucleus 1 and written into the nucleus1.sam file; If the reads are uniquely aligned to the nucleus2 reference genome, the reads will be classified as nucleus2 and written to the nucleus2.sam file.

8. The method for analyzing internuclear allele expression according to claim 1, characterized in that: After combining the obtained similarity score with the gene annotation information and establishing the inter-nuclear allele relationship between the two nuclei through LightGBM model analysis, the method further includes: For genes for which the inter-nuclear allele relationship cannot be established by the LightGBM model, an association is established with NA, which is used to indicate the absence of inter-nuclear alleles.

9. A method and device for analyzing internuclear allele expression, characterized in that: The method comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the computer program implements the steps of the method according to any one of claims 1 to 8 when executed by the processor.

10. A storage medium, characterized in that: The storage medium stores a computer program, which implements the steps of the method according to any one of claims 1 to 8 when executed by a processor.

Images