A method and device for intergenic allelic expression analysis and a storage medium

Abstract

The application discloses a dikaryon inter-nuclear allelic expression analysis method and device and a storage medium, relates to the field of bioinformatics transcriptome analysis of dikaryophyte fungi, and the method comprises the following steps: performing alignment analysis on sequencing data of dikaryophyte fungi and two nuclear reference genomes respectively; performing nuclear source classification on each read in the sequencing data by using an SVM model; taking all gene sequences of nucleus 1 as a query set, taking gene sequences of nucleus 2 as a target set, and using a BLAST tool to construct a similarity score matrix between the gene sequences; combining the similarity score with gene annotation information, analyzing and establishing an inter-nuclear allelic relationship by using a LightGBM model; generating a counts file, performing FPKM standardization processing, combining the expression after the standardization with an inter-nuclear allelic relationship index for analysis, and obtaining an analysis result. Through the application, inter-nuclear allelic analysis of dikaryophyte fungi is realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of bioinformatics transcriptome analysis of dikaryotic fungi, and particularly relates to a method and device for analyzing inter-nuclear allelic expression and a storage medium. BACKGROUND

[0002] Fungi are a kind of eukaryotic, spore-producing, and non-photosynthetic eukaryotic organisms, including molds, yeasts, mushroom fungi, and other mushroom species known to humans. By 2024, about 150,000 fungal species have been discovered and recorded by the scientific community, including yeasts, molds, and various mushroom species. However, dikaryotic fungi, such as the Basidiomycota phylum of fungi, have only been noticed by people in the past decade. Dikaryons can generate monokaryons through sexual and asexual reproduction, and compatible monokaryons can form dikaryons with lock-like unions through plasmogamy. Current research on edible fungal monokaryons has expanded the understanding of the biological properties of monokaryons and dikaryons. Dikaryotic cells exist as a haploid dikaryon, which has many biological properties different from diploids. Many researchers have conducted studies on the biological properties of dikaryons, such as hierarchical order, fitness, and stability, which have brought many insights into the maintenance of fungal species, hybrid breeding innovation, and the construction of genetic research paradigms. The analysis of the incrementality, fertility, and totipotency of dikaryons in fungal species can provide insights into the long-term maintenance of edible fungal species. The dominance and complementary relationship between the nuclei of dikaryotic cells may have important implications for hybrid breeding innovation. The different expression levels of inter-nuclear alleles in dikaryotic cells provide new ideas for studying edible fungi from the perspective of dikaryons. Research on edible fungal dikaryons may form a more discipline-specific research paradigm for edible fungi.

[0003] In the current research progress of dikaryotic fungi, researchers have proved that there is a certain order between the two nuclei of dikaryotic fungi, which affects the biased segregation phenomenon of cell development and the expression mechanism of inter-nuclear alleles in dikaryotic fungi, that is, the alleles existing on the two sets of chromosomes in dikaryotic cells exhibit 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 dikaryotic fungi as monokaryons, and to use the reference genome of one of the nuclei or the corresponding reference genome of fungi on NCBI to align the reads, ignoring the contribution of the two nuclei in dikaryotic fungi to gene expression and their order, resulting in incomplete analysis of dikaryotic fungi and the inability to discover the biological laws behind dikaryons. In order to more completely analyze the transcriptome of dikaryotic fungi and explore the biological laws between the two nuclei of dikaryotic fungi, the present application provides a unique bioinformatics analysis method for dikaryotic fungi. SUMMARY

[0004] The purpose of the present application is to provide a method and device for analyzing inter-nuclear allele expression and a storage medium, in order to study the contribution of each nucleus to gene expression and the order of biological laws, and to solve the incomplete analysis and the lack of inter-nuclear allele analysis method of dikaryotic fungi in the background art.

[0005] To achieve the above object, the present application adopts the following technical solution:

[0006] According to one aspect of the present application, a method for analyzing inter-nucleus allelic expression is provided, comprising: obtaining sequencing data of a dikaryotic fungus, and performing quality control on the sequencing data; performing alignment analysis on the sequencing data after quality control with reference to two nuclear reference genomes of the dikaryotic fungus, respectively, to generate two sam files, wherein the two nuclear reference genomes of the dikaryotic fungus are nucleus 1 reference genome and nucleus 2 reference genome, respectively, and the two sam files record alignment information of each read in the sequencing data on the two nuclear reference genomes, respectively; performing data fusion on the two sam files obtained after alignment to generate a summary sam file containing all alignment information; performing nucleus origin classification on each read in the sequencing data according to the alignment information in the summary sam file using a trained SVM model, and classifying reads from nucleus 1 into a nucleusl.sam file and reads from nucleus 2 into a nucleus2.sam file; extracting the unique ID, chromosome name, start position and end position of each gene according to the genome annotation file of the two nuclei of the dikaryotic fungus; extracting the sequence of the corresponding gene from the nucleus 1 reference genome and the nucleus 2 reference genome according to the position information provided in the genome annotation file; using the gene sequences of all genes in the nucleus 1 genome annotation file as a query set and the gene sequences of all genes in the nucleus 2 genome annotation file as a target set, performing pairwise sequence alignment of the genes in the nucleus 1 genome annotation file and the genes in the nucleus 2 genome annotation file and their corresponding base sequences using the BLAST tool to obtain similarity scores between gene sequences and construct a score matrix; combining the obtained similarity scores with the gene annotation information to analyze and establish the inter-nucleus allelic relationship between the two nuclei by a LightGBM model; converting the nucleusl.sam and nucleus2.sam files after nucleus origin classification into computer-readable nucleusl.bam and nucleus2.bam files; performing quantitative expression analysis on the nucleusl.bam and nucleus2.bam files using the corresponding genome annotation files of the nuclei, respectively, to obtain expression information of each gene and generate counts files; performing FPKM standardization processing on the generated counts files, combining the standardized expression of nucleusl and nucleusl2 with the index of the inter-nucleus allelic relationship, analyzing the standardized expression data of the genes in nucleus 1 and their inter-nucleus alleles in nucleus 2, comparing the expression differences between the two inter-nucleus alleles, or judging the nuclear order relationship in the dikaryotic fungus according to the standardized expression of the inter-nucleus alleles, and obtaining the order relationship of the expression of different genes between the nuclei.

[0007] Optionally, the quality control of the sequencing data comprises: quality control of the obtained sequencing data of the dikaryotic fungi by bioinformatics quality control technology; filtering of reads with base quality score lower than a preset quality score threshold and read length lower than a preset length threshold in the sequencing data, wherein the base quality score is calculated by the formula Q = -10 log 10 P, wherein 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 comprises: read alignment position, alignment score, alignment quality, number of base mismatches, read coverage length, and insert length.

[0009] Optionally, the quality-controlled sequencing data is respectively aligned with the two nuclear reference genomes of the dikaryotic fungi to generate two sam files, and the two sam files obtained after alignment are fused to generate a summary sam file containing all alignment information, which comprises: aligning the quality-controlled sequencing data with the nuclear 1 reference genome to obtain alignment file 1 sam file; aligning the quality-controlled sequencing data with the nuclear 2 reference genome to obtain alignment file 2 sam file; extracting alignment information fields from the two sam files of alignment file 1 and alignment file 2, and placing the extracted field information into a summary sam file, wherein the alignment information fields comprise mandatory fields and optional fields, the mandatory fields comprise QNAME, RNAME, POS, MAPQ, CIGAR, MRNM and PNEXT, and the optional fields are other fields that are helpful to SVM model classification in addition to the mandatory fields, wherein QNAME represents the name of the read, RNAME represents the name of the aligned reference sequence, POS represents the starting position of the alignment to the reference sequence, MAPQ represents the alignment quality score reflecting the reliability of the alignment position, CIGAR represents the compressed 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 on the reference genome in the paired alignment.

[0010] Optionally, before performing the nucleus origin classification on each read in the sequencing data according to the alignment information in the summary sam file using the trained SVM model, the method further comprises training the SVM model by classifying the alignment information of the reads in the sequencing data to determine whether each read is from nucleus 1 or nucleus 2, thereby realizing the classification of the nucleus origin, wherein the training process comprises: 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 rule is:

[0011] ;

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

[0013] MAPQ1 and MAPQ2 represent the alignment quality of reads aligned to reference genome 1 and 2 respectively; Δ represents the alignment score difference, Δ = AS1-AS2, AS1 is the alignment score aligned to reference genome 1, and AS2 is the alignment score aligned to reference genome 2; T Δ represents the significance threshold of the alignment score difference, used to determine whether the score has a significant difference; T MAPQ represents the significance threshold of the alignment quality difference, used to determine whether the alignment quality has a significant difference; parameter optimization: using cross-validation to determine the kernel function and hyperparameters of the SVM model to optimize the classification performance of the model; training model: training the SVM model with training data labeled as Category to obtain an SVM model based on input features for classifying the nucleus origin of reads.

[0014] Optionally, according to the alignment information in the summary sam file, the nucleus source classification of each read in the sequencing data is performed, and the reads from the nucleus 1 are classified into the nucleusl.sam file, and the reads from the nucleus 2 are classified into the nucleus2.sam file, and the method further comprises: according to the alignment information of the reads on the two-nucleus reference genome, a confidence value of each read belonging to the nucleus 1 or the nucleus 2 is calculated; when the confidence value of the read belonging to the nucleus 1 is greater than that of the nucleus 2 and is higher than a certain threshold value, the read is classified into the nucleus 1 and written into the nucleusl.sam file; when the confidence value of the read belonging to the nucleus 2 is greater than that of the nucleus 1 and is higher than a certain threshold value, the read is classified into the nucleus 2 and written into the nucleus2.sam file; when the confidence values of the read belonging to the nucleus 1 and the nucleus 2 are both lower than a preset threshold value, the read is discarded.

[0015] Optionally, according to the alignment information in the summary sam file, the nucleus source classification of each read in the sequencing data is performed, and the reads from the nucleus 1 are classified into the nucleusl.sam file, and the reads from the nucleus 2 are classified into the nucleus2.sam file, and the method further comprises: if the read can only find a unique alignment position on one reference genome, the classification result is marked according to the reference genome, and specifically, if the read is uniquely aligned on the nucleus 1 reference genome, the read is classified into the nucleus 1 and written into the nucleusl.sam file; if the read is uniquely aligned on the nucleus 2 reference genome, the read is classified into the nucleus 2 and written into the nucleus2.sam file.

[0016] Optionally, after the inter-nucleus allelic relationship between the two nuclei is analyzed and established by the LightGBM model by combining the obtained similarity score with the gene annotation information, the method further comprises: for a gene that cannot establish the inter-nucleus allelic relationship by the LightGBM model, an association with NA is established, and the NA is used to represent no inter-nucleus allele.

[0017] According to another aspect of the present application, an inter-nucleus allelic expression analysis device is also provided, which comprises a memory, a processor and a computer program stored on the memory and executable on the processor, and the computer program implements the steps of the method as described above when executed by the processor.

[0018] According to another aspect of the present application, a storage medium is also provided, and the storage medium stores a computer program, and the computer program implements the steps of the method as described above when executed by a processor.

[0019] Compared with the traditional single-core analysis mode, the integrity of the dikaryon fungus analysis is improved, and the data processing and classification are performed through the SVM (Support Vector Machine) model and the LightGBM model two machine learning modes, so that the analysis result is more accurate and reliable, and the expression mode between the inter-nuclear alleles of the dikaryon fungus can be clearly and obviously indicated through the method in the application, and more biological rules between the inter-nuclear allele groups are explored. BRIEF DESCRIPTION OF DRAWINGS

[0020] The accompanying drawings, which are included to provide a further understanding of the application and are incorporated in and constitute a part of this application, illustrate embodiments of the application and together with the description serve to explain the application. In the drawings:

[0021] Figure 1 is a flowchart of an inter-nuclear allele expression analysis method in embodiment 1 of the application;

[0022] Figure 2 is a flowchart of an inter-nuclear allele expression analysis method in embodiment 2 of the application;

[0023] Figure 3 is a schematic diagram of the percentage of genes with established allele relationships in embodiment 2 of the application;

[0024] Figure 4 is a Wean diagram of the number of genes with judged nuclear order relationships in embodiment 2 of the application. DETAILED DESCRIPTION

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

[0026] Secondly, the "one embodiment" or "embodiment" referred to herein means that the specific features, structures or characteristics can be included in at least one implementation of the present application. The "in one embodiment" appearing in different places in the specification does not mean the same embodiment, nor is it an embodiment that is separate or alternative to other embodiments.

[0027] It should be noted that the terms "first", "second", etc. in the specification and claims of the present application and in the above drawings are used to distinguish similar objects, and do not necessarily have to be used to describe a specific order or sequence. It should be understood that the data thus used can be interchanged under appropriate circumstances, so that the embodiments of the present application described herein can be implemented in an order other than that illustrated or described herein. In addition, the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, product or device including a series of steps or units does not have to be limited to only those steps or units clearly listed, but can include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0028] Embodiment 1

[0029] In this embodiment, a method for analyzing inter-nuclear allelic expression is provided.

[0030] Reference Figure 1 , Figure 1 is a flowchart of a method for analyzing inter-nuclear allelic expression in an embodiment of the present application, as shown in Figure 1 , the method comprises the following steps:

[0031] S1, obtaining sequencing data of a dikaryotic fungus, and performing quality control on the sequencing data;

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

[0033] Then, the obtained dikaryotic fungus sequencing data are subjected to quality control by bioinformatics quality control technology, such as fastp, fastqC and other bioinformatics software, and low-quality reads are filtered to obtain high-quality sequencing data for analysis.

[0034] The specific quality control filtering strategy includes: according to the formula Q = -10·log 10P, where P represents the probability of the base being misidentified; filtering reads with base quality scores lower than a preset quality score threshold, which is usually set according to the requirement of data quality, and common quality score thresholds are Q20 (i.e. 1% error rate) or Q30 (i.e. 0.1% error rate), if the quality score of some bases in the reads is lower than the set quality score threshold, then these reads will be considered as low quality and filtered out. In addition, shorter reads contain less information and are more susceptible to sequencing errors, therefore, in the embodiments of the present application, reads with read lengths lower than a preset length threshold are also filtered, and the read length threshold can be set according to specific analysis requirements, and the embodiments of the present application do not make specific limitations.

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

[0036] S2, the sequencing data after quality control is respectively compared and analyzed with the two nuclear reference genomes of the dikaryophytic fungus to generate two sam files;

[0037] The two nuclear reference genomes of the dikaryophytic fungus are nucleus 1 reference genome (nucleus 1) and nucleus 2 reference genome (nucleus 2), and these reference genomes are known genome sequences, which are used for comparison with sequencing data.

[0038] The high-quality sequencing data after quality control is compared and analyzed with the two nuclear reference genomes of the dikaryophytic fungus using bioinformatics alignment tools such as hisat2, and the comparison process will generate two sam (Sequence Alignment / Map) files, each file corresponding to the comparison results of one reference genome, wherein the comparison file 1. sam file contains the comparison information of sequencing data with nucleus 1 reference genome, and the comparison file 2. sam file contains the comparison information of sequencing data with nucleus 2 reference genome. The comparison information includes read comparison position, comparison score, comparison quality, number of base mismatches, read coverage length, insert length, etc. Among them, the comparison position is the specific position of the sequencing reads on the reference genome; the comparison score is the matching degree score between the reads and the reference genome; the comparison quality represents the reliability or accuracy of the comparison result; the number of base mismatches refers to the number of bases that do not match between the reads and the reference genome; the read coverage length refers to the length of the continuous region successfully compared by the reads on the reference genome; and the insert length refers to the length of the fragment that is more than the reference genome in the reads.

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

[0040] In the two sam files, find and extract necessary alignment information fields, including but not limited to QNAME, RNAME, POS, MAPQ, CIGAR, MRNM, PNEXT, and optional field information that helps SVM classification. The extracted field information in the two sam files is merged into a new sam file to obtain a summary sam file containing all alignment information, wherein 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, reflecting the reliability of the alignment position; CIGAR represents the compressed 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 (1 for the starting point) of the partner read on the reference genome in the paired alignment.

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

[0042] The merged sam file retains all reads' QNAME, RNAME, POS, MAPQ, CIGAR, MRNM, PNEXT, and other alignment information, which helps further analyze alignment quality, sequencing accuracy, and genome coverage, providing necessary support for subsequent feature extraction and classification model construction.

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

[0044] Since the two cell nuclei of dikaryotic fungi are from different parents, and the base in the gene sequence will accumulate gene base mutations during the generation development, resulting in base differences in the gene sequence of the two nuclei. According to the principle that the same gene on the two reference genomes still has base differences in the sequence, the SVM (Support Vector Machine) model trained by the fungal data set classifies the source of each read according to the QNAME, RNAME, POS, MAPQ, CIGAR, MRNM, PNEXT and other mandatory field information in the summary sam file and the optional field information (determined according to different alignment conditions) that helps SVM classification, and gives a confidence value for each nuclear source. When the confidence value of reads in nucleus 1 is greater than that in nucleus 2 and higher than a certain threshold, the reads will be classified as nucleus 1 and written into the nucleus1.sam file. When the confidence value of reads in nucleus 2 is greater than that in nucleus 1 and higher than a certain threshold, the reads will be classified as nucleus 2 and written into the nucleus2.sam file. When the confidence values of reads in nucleus 1 and nucleus 2 are both 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, and there is no alignment result on the other reference genome, the classification result can be directly marked according to the reference genome. Specifically, if the read is uniquely aligned on the nucleus 1 reference genome, the read will be classified as 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 as nucleus 2 and written into the nucleus2.sam file.

[0046] Reads from nucleus 1 (nucleus 1) will be classified into a sam file, and reads from nucleus 2 (nucleus 2) will be classified into another sam file. Reads that cannot be classified are filtered.

[0047] Further, before classifying the 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 comprises training the SVM model by classifying the alignment information (including alignment position, alignment score, alignment quality, etc.) of the reads in the sequencing data to determine whether each read is from nucleus 1 or nucleus 2, thereby realizing the classification of the nuclear source. The training process includes:

[0048] (1) Feature normalization: normalize the input features to ensure the same numerical range for different features, such as Δ and |MAPQ1-MAPQ2| are standardized according to z-score;

[0049] (2) Model construction: use SVM algorithm to construct a classification model, the input features are RNAME1, RNAME2, Δ, |MAPQ1-MAPQ2|, define classification rules, including basic classification rules and complex classification rules, wherein,

[0050] Basic classification rule: if a read can only be aligned to one genome, directly mark the classification result as 1 or 2, such reads usually have high reliability and can be used as clear samples in the training set;

[0051] Complex classification rule: if a read is aligned to two genomes, i.e. there are alignment results on two genomes, use the alignment score difference Δ and the alignment quality difference

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

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

[0054] When Δ < -T Δ , it means that the score of aligning to 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 , the classification result is Undecided, indicating uncertain classification; for Undecided category, rules can be designed separately or a third category can be added as an output of SVM classification;

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

[0057]

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

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

[0060] AS2 is the alignment score of reference genome 2; T Δ represents the significance threshold of the difference of alignment score, used to determine whether the score has a significant difference; T MAPQ represents the significance threshold of the difference of alignment quality, used to determine whether the alignment 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) Train the model: train the SVM model with the training data labeled as Category to obtain an SVM model that can classify the reads based on the input features.

[0063] S5, according to the genome annotation file of the two nuclei of the dikaryophyte fungus, extract the unique ID, chromosome name, start position and end position of each gene;

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

[0065] S6, according to the position information provided in the genome annotation file, extract the sequence of the corresponding gene from the reference genome of nucleus 1 and nucleus 2 respectively;

[0066] When constructing the allelic relationship of genes between the two nuclei of the dikaryophyte fungus through the genome annotation files of the two nuclei, record the start position and end position on the corresponding chromosome and ensure that each gene has a unique ID, then use the reference genome file (FASTA format) to locate the corresponding reference genome according to the chromosome name, start position and end position of the gene in the annotation file, and extract the sequence of the corresponding gene from the reference genome of nucleus 1 and nucleus 2 respectively.

[0067] S7, use the gene sequences of all genes in the nucleus 1 genome annotation file as the query set, and the gene sequences of all genes in the nucleus 2 genome annotation file as the target set, use the BLAST tool to perform pairwise sequence alignment of the genes in the nucleus 1 genome annotation file and the genes in the nucleus 2 genome annotation file and their corresponding base sequences, obtain the similarity score between the gene sequences, and construct the score matrix;

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

[0069] Take the gene sequences of all genes in the nucleic acid 1 genome annotation file as the query sequence set, and the gene sequences of all genes in the nucleic acid 2 genome annotation file as the target sequence set. BLAST each pair of gene sequences to obtain the similarity score (such as the alignment score) between the gene sequences. Finally, take all gene IDs of nucleic acid 1 as row names and all gene IDs of nucleic acid 2 as column names to obtain a score matrix. Each element in the score matrix represents the similarity score between a gene in nucleic acid 1 and a gene in nucleic acid 2.

[0070] S8, combine the obtained similarity score with the gene annotation information, and analyze and establish the inter-nucleus allelic relationship between the two nuclei through the LightGBM model;

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

[0072] LightGBM (Light Gradient Boosting Machine) is a distributed gradient boosting model based on decision tree algorithm. This step uses the LightGBM model to further analyze and infer the allelic relationship between genes that have not established a relationship. When the difference between the BLAST similarity and the similarity of all genes in another nucleus is below a certain threshold, it is determined that the inter-nucleus allelic relationship cannot be established. For genes that still cannot be established, NA is established, indicating that there is no inter-nucleus allelic relationship, and it is believed that the gene is only regulated by one nucleus for expression.

[0073] S9, convert the nucleic acid 1.sam and nucleic acid 2.sam files classified by the nucleus source into computer-readable nucleic acid 1.bam and nucleic acid 2.bam files;

[0074] The sam file can be converted into a computer-readable bam file by bioinformatics tools such as Samtools for subsequent analysis.

[0075] S10, use the genome annotation file corresponding to the nucleus to perform quantitative expression analysis on nucleic acid 1.bam and nucleic acid 2.bam respectively to obtain the expression information of each gene and generate a counts file;

[0076] The quantitative expression analysis is performed by comparing the reads and the genome annotation file, and the number of reads aligned to each gene is counted to reflect the transcription level of the gene in nucleus 1 and nucleus 2. The tool for bam file quantitative expression analysis, such as featureCounts, is used to perform quantitative expression analysis on the nucleus1.bam and nucelus2.bam files respectively with the nucleus1.gff(gtf) and nucleus2.gff(gtf) genome annotation, and the counts data of each nucleus is analyzed.

[0077] S11, the generated counts file is subjected to FPKM normalization processing, the normalized expression of nucleus1 and nucleus2 is combined with the inter-nuclear allele relationship index, the normalized expression data of the gene in nucleus 1 and its inter-nuclear allele in nucleus 2 are analyzed, and the expression difference between the two inter-nuclear alleles is compared, or the order relationship of the nuclei in the dikaryophyte is judged according to the normalized expression of the inter-nuclear allele, and the order relationship of the expression of different genes between the nuclei is obtained.

[0078] The quantitative files of nucleus1 and nucleus2 are subjected to FPKM (Fragments Per Kilobase of transcript per Million mapped reads) data normalization to eliminate the influence of the number difference between reads and the library during sequencing, the normalized expression of nucleus1 and nucleus2 FPKM is combined with the corresponding inter-nuclear allele relationship index, the normalized expression data of the gene in nucleus 1 and its inter-nuclear allele in nucleus 2 are compared, and the expression difference between the two inter-nuclear alleles is analyzed and compared by combining statistical methods such as paired t-test.

[0079] Based on the expression difference between the two inter-nuclear alleles, or according to the normalized expression of the inter-nuclear allele, the order relationship of the nuclei in the dikaryophyte is judged, and based on the analysis result, the order relationship of the expression of different genes between the nuclei is obtained, such as which gene is expressed higher in one nucleus and which gene is expressed higher in the other nucleus, and the “dominant” or “subordinate” status of the nucleus is inferred accordingly.

[0080] The inter-nuclear allele expression analysis method of the dikaryophyte fungus of the present application perfects the integrity of the dikaryophyte fungus analysis, and the data processing and classification are performed by two machine learning modes of SVM (Support Vector Machine) model and LightGBM model, so that the analysis result is more accurate and reliable. The expression mode between the inter-nuclear alleles of the dikaryophyte fungus can be clearly and clearly represented by the method in the present application, and more biological rules between the inter-nuclear allele groups are explored.

[0081] Example 2

[0082] The technical solutions of the present application will be further described in detail below in combination with one specific embodiment and the accompanying drawings. It should be understood that the following embodiments are only used to explain the present application and do not limit the present application.

[0083] The implementation scenario in this embodiment is Auricularia heimuer. The inter-nucleus allelic expression analysis method of the present application is used to perform double-nucleus inter-nucleus allelic transcriptome analysis on Auricularia heimuer.

[0084] Auricularia heimuer (scientific name Auricularia heimuer), also known as Yun'er, Auricularia, Auricularia, Auricularia, Auricularia, Auricularia, Auricularia, Auricularia, is a kind of edible fungus of Auriculariaceae and Auricularia. Auricularia is a kind of bi-nuclear fungus. The life history of Auricularia includes the germination of basidiospores into single-nucleus mycelium or the formation of sickle-shaped conidia, and then the germination of conidia into single-nucleus mycelium. Single-nucleus mycelium of different sexes combines to form double-nucleus mycelium, which continuously proliferates through lock-like combination. The double-nucleus mycelium that reaches the physiological maturity stage forms primordium of fruiting body on the surface of the substrate and develops into fruiting body.

[0085] In this embodiment, the reference genome of nucleus1 has a genome size of about 41.12 Mb, a GC content of 52.68%, and an N50 of 1168003 bp, and the reference genome of nucleus2 has a genome size of about 40.52 Mb, a GC content of 54.32%, and an N50 of 1302490 bp. The genome size refers to the DNA content of a genome; the GC content refers to the proportion of guanine (Guanine) and cytosine (Cytosine) in the genome; and the N50 refers to the length of the shortest fragment that exceeds half (50%) of the total length of all fragments when all assembled fragments are sorted by length from long to short.

[0086] As shown in Figure 2 , the inter-nucleus allelic expression analysis method in the embodiment of the present application includes the following steps:

[0087] Step 1: Use fastp (0.23.4) to perform quality control analysis and processing on the original sequencing data. The high-quality sequencing data after quality control is aligned to the reference genomes of double nuclei (nucleus1 and nucleus2) using Hisat2 (version 2.2.1), and the global alignment rate of nucleus1 is 97.72%, and the global alignment rate of nucleus2 is 92.53%;

[0088] Step 2, according to the above compared data, filter the unaligned data, keep the valid reads, the total reads after quality control by fastp tool are 8.23M, the remaining reads in nucleus1 after filtering are 8.035M, and the remaining reads in nucleus2 are 7.703M;

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

[0090] Step 4, using the trained SVM classification model, classify each read in the sam file after step 3 fusion, write the reads information of nucleus1 into nucleus1.sam file, and write the reads information of nucleus2 into 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 nucleus1 and nucleus2 is 5704710 and 6696834 respectively;

[0091] Step 5, using samtools view(1.9) to convert sam file to binary bam file;

[0092] Step 6, according to the genome annotation file of the two nuclei of the dikaryophyte fungus, extract the unique ID, chromosome name, start position and end position of each gene, and cut the gene sequence of the reference genome of nucleus1 and nucleus2 respectively;

[0093] Step 7, using BLAST tool to pairwise align the genes in the genome annotation file from nucleus1 with the genes in the genome annotation from nucleus2, obtaining a 15420*14903 scoring matrix for evaluating the similarity between gene sequences, and constructing the scoring matrix;

[0094] Step 8, give the 15420*14903 scoring matrix score and the annotation information of each nucleus gene in its genome annotation file to the LightGBM machine learning model which has been preliminarily trained locally to establish the relationship between alleles between the two nuclei, and discard the reads that cannot establish the relationship between alleles between the two nuclei; refer to Figure 3 , Figure 3 The Percentage of Genes with Established Allele Relationships in this embodiment is shown in the schematic diagram as follows: Figure 3As shown, the proportion of genes with allele relationships is much higher than that of genes without allele relationships.

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

[0096] Step 10: Based on the processing results of featureCounts, the expression levels of each gene are standardized using FPKM (Fragments Per Kilobase of Transcript per Million Mapped Reads) 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 bases; N = total number of all valid fragments; L = length of the target gene (in bases); M = total length of all mapped fragments (in bases).

[0097] Step 11: Based on the FPKM-normalized expression data and combined with the internuclear allele relationship index, analyze the allele expression differences between the two nuclei. The following statistical methods are used to compare the allele expression between the nuclei to obtain the analysis results:

[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 normalized expression level of the target gene in nucleus 2; A represents the set judgment threshold; Mean(FPKMN ucleus1 FPKMN ucleus2 () indicates 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 Venn diagram shows the number of gene expression dominance in nucleus 1 and nucleus 2. Nucleus1Dominant indicates the gene set whose expression in nucleus 1 is significantly higher than that in nucleus 2, Nucleus2Dominant indicates the gene set whose expression in nucleus 2 is significantly higher than that in nucleus 1, and Undecided indicates that the gene set whose expression in nucleus 1 and nucleus 2 is not significantly different and whose dominance cannot be determined.

[0101] In the embodiment, an inter-nuclei allelic expression analysis device is also provided for implementing the above-described embodiments and preferred embodiments, which have been described and will not be repeated here.

[0102] Alternatively, specific examples in the embodiment can refer to the examples described in the above-described embodiments and optional embodiments, which will not be repeated here.

[0103] The embodiment of the present application further provides a storage medium, and the storage medium stores a computer program.

[0104] Specific examples in the embodiment can refer to the examples described in the above-described embodiments and optional embodiments, which will not be repeated here.

[0105] Alternatively, in the embodiment, the storage medium can include but is not limited to a U disk, a read-only memory (ROM), a random access memory (RAM), a mobile hard disk, a magnetic disk or an optical disk and various storage media that can store computer programs.

[0106] The above-mentioned serial numbers of the embodiments of the present application are only for description, and do not represent the advantages and disadvantages of the embodiments.

[0107] In the above-described embodiments of the present application, the description of each embodiment has its own focus, and the parts not described in detail in a certain embodiment can be referred to the relevant description of other embodiments.

[0108] The above-mentioned only is the preferred embodiment of the present application, it should be pointed out, for the ordinary skilled in the art, without departing from the principles of the present application, can make a number of improvements and refinements, these improvements and refinements should also be considered as the protection scope of the present application.

Claims

1. A method for analyzing internuclear allele expression, characterized in that, include: Sequencing data of binucleate fungi were obtained, and the sequencing data were subjected to quality control. The quality-controlled sequencing data were compared and analyzed with the two nuclear reference genomes of binuclear fungi to generate two SAM files. The two nuclear reference genomes of binuclear fungi are nuclear reference genome 1 and nuclear reference genome 2, respectively. The two SAM files record the alignment 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, based on the alignment information in the summarized SAM file, each read in the sequencing data is classified by its nuclear origin. Reads from nuclear 1 are assigned to the nucleus1.sam file, and reads from nuclear 2 are assigned to the nucleus2.sam file. Based on the genome annotation files of the two nuclei of binucleate fungi, extract the unique ID, chromosome name, start position, and end position of each gene; Based on the location information provided in the genome annotation file, the sequences of the corresponding genes were extracted from the nuclear 1 reference genome and the nuclear 2 reference genome, respectively. The gene sequences of all genes in the nucleus 1 genome annotation file are used as the query set, and the gene sequences of all genes in the nucleus 2 genome annotation file are used as the target set. The BLAST tool is used to perform pairwise sequence alignment between the genes in the nucleus 1 genome annotation file and the genes and their corresponding base sequences in the nucleus 2 genome annotation file to obtain similarity scores between gene sequences and construct a score matrix. The obtained similarity scores were combined with gene annotation information, and the internuclear allele relationship between the two nuclei was established by analyzing the LightGBM model. Convert the kernel source-classified nucleusl.sam and nucleus2.sam files into computer-readable nucleus1.bam and nucleus2.bam files; Quantitative expression analysis was performed on nucleus1.bam and nucleus2.bam using the genome annotation files corresponding to the nucleus, respectively, to obtain the expression level information of each gene and generate a counts file; The generated counts file was normalized using FPKM. The normalized expression levels of nucleus1 and nucleusl2 were combined with the internuclear allele relationship index to analyze the normalized expression levels of genes in nucleus 1 and their internuclear alleles in nucleus 2. The expression differences between the two internuclear alleles were compared, or the nuclear order relationship in binuclear fungi was determined based on the normalized expression levels of internuclear alleles, thus obtaining the order relationship of different genes expressed in the nuclei.

2. The method for analyzing internuclear allele expression according to claim 1, characterized in that, The quality control of the sequencing data includes: Quality control was performed on the acquired binucleosome fungal sequencing data using bioinformatics quality control techniques. Reads in the sequencing data whose base quality score is lower than a preset quality score threshold and whose read length is lower than a preset length threshold are filtered out. The base quality score is determined by the formula... Q = -10 log 10 P is calculated as follows: Q represents the base quality score, and P represents the probability that a base is misidentified.

3. The method for analyzing internuclear allele expression according to claim 1, characterized in that, The alignment information for 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 quality-controlled sequencing data were aligned with two nuclear reference genomes of binucleate fungi, generating two SAM files. These two SAM files were then fused to produce a single summary SAM file containing all alignment information, including: The quality-controlled sequencing data was compared with the nucleus 1 reference genome to obtain the alignment file 1sam file; The quality-controlled sequencing data was compared with the nucleus 2 reference genome to obtain the alignment file 2sam. Alignment information fields were extracted from two SAM files, alignment file 1 and alignment file 2, and the extracted field information was placed into a summary SAM file. The alignment information fields include required fields and optional fields. The required fields include QNAME, RNAME, POS, MAPQ, CIGAR, MRNM, and PNEXT. The optional fields are other fields that are helpful for SVM model classification besides the required fields. Specifically, QNAME represents the name of the read, RNAME represents the name of the reference sequence to be aligned, POS represents the start position of the alignment to the reference sequence, MAPQ represents the alignment quality score, reflecting the confidence of the alignment position, CIGAR represents the alignment compression information, indicating how the query sequence is aligned with the reference sequence, MRNM represents the name of the reference sequence that corresponds to the partner read in the alignment, and PNEXT represents the start position of the partner read in the alignment on the reference genome.

5. The method for analyzing internuclear allele expression according to claim 1, characterized in that, Before classifying the nuclear origin of each read in the sequencing data using the trained SVM model based on the alignment information in the summarized SAM file, the method further includes training the SVM model. By classifying the alignment information of reads in the sequencing data, the SVM model is trained to determine whether each read originates from nucleus 1 or nucleus 2, thereby achieving nuclear origin classification. The training process includes: Feature standardization: Normalize the input features to ensure that the numerical range of different features is consistent; Model Construction: A classification model is constructed using the SVM algorithm, with input features RNAME1, RNAME2, Δ, and |MAP. Q1 -MAP Q2 The classification rules are as follows: ； Here, Category represents the classification result of the reads, 1 represents from nucleus 1, 2 represents from nucleus 2, and Undecided represents no clear classification; RNAME1 represents the sequence name of the reads aligned to reference genome 1; RNAME2 represents the sequence name of the reads aligned to reference genome 2; MAPQ1 and MAPQ2 represent the alignment quality of reads to reference genomes 1 and 2, respectively; Δ represents the difference in alignment scores, Δ = AS1 - AS2, where AS1 is the alignment score to reference genome 1 and AS2 is the alignment score to reference genome 2; T Δ T represents the significance threshold for the difference in scores, used to determine whether the scores have a significant difference; MAPQ The significance threshold for comparison quality differences is used to determine whether there are significant differences in comparison quality. Parameter optimization: Cross-validation is used to determine the kernel function and hyperparameters of the SVM model in order to optimize the model's classification performance; Training the model: Train the SVM model using training data labeled Category to obtain an SVM model that classifies the kernel source of reads based on input features.

6. The method for analyzing internuclear allele expression according to claim 1, characterized in that, Based on the alignment information in the summarized SAM file, each read in the sequencing data is classified by its nuclear origin. Reads from nucleus1 are assigned to the nucleus1.sam file, and reads from nucleus2 are assigned to the nucleus2.sam file. Based on the alignment information of reads on two nuclear reference genomes, calculate the confidence value for each read to belong to nuclear 1 or nuclear 2; When the confidence value of reads belonging to core 1 is greater than that of core 2 and exceeds a certain threshold, the reads will be classified as core 1 and written to the nucleus1.sam file; When the confidence value of reads belonging to core 2 is greater than that of core 1 and exceeds a certain threshold, the reads will be classified as core 2 and written to the nucleus2.sam file; When the confidence level of reads belonging to both core 1 and core 2 is lower than the preset threshold, the reads will be discarded.

7. The method for analyzing internuclear allele expression according to claim 1, characterized in that, Based on the alignment information in the summarized SAM file, each read in the sequencing data is classified by its nuclear origin. Reads from nucleus1 are assigned to the nucleus1.sam file, and reads from nucleus2 are assigned to the nucleus2.sam file. This also includes: If reads can only find a unique alignment location on a single reference genome, then the classification is based on the markers of that reference genome, specifically as follows: If a read is uniquely aligned to the nucleus1 reference genome, the read will be classified as nucleus1 and written to the nucleus1.sam file; If a read is uniquely aligned to the nucleus2 reference genome, the read 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 gene annotation information and establishing the internuclear allele relationship between the two nuclei using the LightGBM model, the method further includes: For genes whose internuclear allele relationships cannot be established using the LightGBM model, an association is established with NA, where NA represents the absence of internuclear alleles.

9. An apparatus for analyzing internuclear allele expression, characterized in that, It includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the computer program, when executed by the processor, implements the steps of the method as described in any one of claims 1 to 8.

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

Images