Crop whole genome phenotype prediction method and system fused with environmental indicator gene

Abstract

The invention relates to the technical field of bioinformatics, and provides a crop whole genome phenotype prediction method and system fused with an environmental indicator gene, and the method comprises the following steps: collecting re-sequencing data, and carrying out genome comparison to obtain variation site data; performing whole genome association analysis by using the variation site data to obtain phenotype association site information; carrying out gene expression quantity measurement on samples of the crop population material in different environments to obtain gene expression quantity data; performing differential expression analysis on the gene expression quantity data to screen environmental indicator genes to obtain an environmental indicator gene set; constructing a phenotype prediction model of double-branch fusion; and predicting a to-be-predicted material through the phenotype prediction model to obtain phenotype prediction results for different environments. According to the method, environmental factors are incorporated into the whole genome selection model, so that the phenotype prediction precision in different environments is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a crop whole genome phenotype prediction method and system integrating environmental indicator genes. Background Art

[0002] Whole-genome selection is a highly efficient breeding technique that uses genome-wide molecular marker information to predict individual breeding values. With the recent advancements in sequencing technology and the reduction in its cost, it has been widely used in animal and plant breeding. Currently, mainstream whole-genome selection models rely solely on genotypic data, predicting phenotypes solely at the genetic level without considering the impact of environmental factors on the final phenotype. This results in a discrepancy between the predicted phenotypic values ​​and the actual phenotypic values ​​in field environments.

[0003] Currently, a common solution is to incorporate real-world environmental indicators of the planting site (such as temperature, precipitation, and sunshine) into the prediction model to improve prediction accuracy. However, the environmental factors in real fields are numerous and complex, the observation costs are high, and existing environmental indicators cannot fully reflect the impact of the environment on individuals. In addition, the level of mRNA in an organism is an important indicator of an individual's life activities. When an individual perceives changes in the surrounding environment, it will change the expression levels of some of its own genes through various pathways and methods to achieve the purpose of adapting to the environment. However, due to the nature of plants being unable to move on their own, their response to the environment is more dependent on the regulation of their own gene expression levels. Summary of the Invention

[0004] The present invention provides a crop whole genome phenotype prediction method and system integrating environmental indicator genes, to address the deficiencies of the prior art.

[0005] The present invention provides a crop genome-wide phenotype prediction method integrating environmental indicator genes, comprising:

[0006] S1: Collect resequencing data and perform genome alignment on the resequencing data to obtain variant site data;

[0007] S2: Perform genome-wide association analysis using the variant site data to obtain phenotype-associated site information;

[0008] S3: Gene expression data were obtained by measuring the gene expression of samples of crop population materials under different environments based on transcriptome sequencing technology;

[0009] S4: performing differential expression analysis on the gene expression data to screen environmental indicator genes and obtain an environmental indicator gene set;

[0010] S5: constructing a double-branch fusion phenotype prediction model using the environmental indicator gene set, the gene expression data, and the phenotype-associated site information;

[0011] S6: Predicting the material to be predicted using the phenotypic prediction model to obtain phenotypic prediction results for different environments.

[0012] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, step S1 further comprises:

[0013] S11: Collect resequencing data of crop population materials to obtain original sequencing data;

[0014] S12: performing quality control on the original sequencing data to obtain quality control sequencing data;

[0015] S13: Align the quality control sequencing data to the reference genome to obtain variant site data.

[0016] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, step S12 specifically includes:

[0017] Discard reads whose base quality values ​​at the beginning of the reads are lower than the preset threshold;

[0018] Discard reads whose base quality values ​​at the end of the reads are lower than the preset threshold;

[0019] After clipping, reads whose length is less than a preset threshold are discarded.

[0020] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, step S13 further comprises:

[0021] S131: Align the quality control sequencing data to the reference genome using the BWA-MEM algorithm to obtain preliminary alignment results;

[0022] S132: Screening the preliminary comparison results according to comparison quality and distance standards to obtain valid comparison results;

[0023] S133: Perform variation detection on the valid comparison results to obtain a single sample variation identification file for a single material;

[0024] S134: Merge the single sample variation identification files of all materials to obtain a population variation information file;

[0025] S135: Filtering the population variation information file according to minor allele frequency and data missing rate standards to obtain variation site data.

[0026] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, step S2 further comprises:

[0027] S21: performing principal component analysis on the variant site data to obtain population structure information;

[0028] S22: Calculating a kinship matrix based on the group structure information to obtain group kinship information;

[0029] S23: performing association analysis on the phenotypic data and the variant site data based on the population structure information and the population kinship information to obtain phenotype-associated site information.

[0030] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, step S3 further comprises:

[0031] S31: Sampling the same crop population materials in different environments and extracting RNA to obtain sample RNA;

[0032] S32: performing high-throughput sequencing on the sample RNA to obtain raw transcriptome sequencing data;

[0033] S33: Perform quality control on the original transcriptome sequencing data to obtain transcriptome quality control sequencing data;

[0034] S34: aligning the transcriptome quality control sequencing data to a reference genome, and extracting unique aligned reads;

[0035] S35: Calculating expression data for the uniquely aligned reads to obtain gene expression data.

[0036] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, step S4 further comprises:

[0037] S41: Perform differential expression analysis on gene expression data to obtain differentially expressed gene sets between different environments;

[0038] S42: Screening genes with significant differences from the differentially expressed gene set to obtain candidate environmental indicator genes;

[0039] S43: analyzing the expression stability of the candidate environmental indicator genes in different environments, screening genes with stable differential expression using a preset threshold, and obtaining a stable differential expression gene set;

[0040] S44: constructing an environmental prediction model based on the XGBoost algorithm, ranking the importance of the stable differentially expressed gene set, screening the stable differentially expressed genes according to the ranking results and the model prediction accuracy, and obtaining an environmental indicator gene set.

[0041] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, step S5 further comprises:

[0042] S51: performing structural conversion on the phenotype-associated site information according to a preset encoding method to obtain genotype encoding data;

[0043] S52: normalizing the expression levels of multiple genes in the environmental indicator gene set to obtain environmental indicator characteristic data;

[0044] S53: using a convolutional neural network to learn the genotype encoding data to obtain a genotype feature representation;

[0045] S54: Learning the environmental indicator feature data through a fully connected neural network to obtain an environmental feature representation;

[0046] S55: Based on the feature fusion layer, the genotype feature representation and the environmental feature representation are integrated to obtain a fused feature representation;

[0047] S56: Outputting the phenotype prediction result through the joint decision layer according to the fusion feature representation to obtain a phenotype prediction model.

[0048] According to a crop genome-wide phenotype prediction method integrating environmental indicator genes provided by the present invention, the phenotype prediction result in step S6 is a predicted phenotypic value, which represents the phenotypic quantitative value of the material to be predicted under the target environment.

[0049] The present invention also provides a crop genome-wide phenotype prediction system integrating environmental indicator genes, comprising:

[0050] Alignment module: used to collect resequencing data and perform genome alignment on the resequencing data to obtain variant site data;

[0051] Association analysis module: used to perform genome-wide association analysis using the variant site data to obtain phenotype-associated site information;

[0052] Measurement module: used to measure gene expression levels of samples of crop population materials under different environments based on transcriptome sequencing technology to obtain gene expression data;

[0053] Screening module: used to perform differential expression analysis on the gene expression data to screen environmental indicator genes,

[0054] Obtaining an environmental indicator gene set;

[0055] Construction module: for constructing a double-branch fusion phenotype prediction model using the environmental indicator gene set, the gene expression data, and the phenotype-associated site information;

[0056] The prediction module is configured as the phenotypic prediction model obtained by training the building module, and is used to predict the material to be predicted and obtain phenotypic prediction results for different environments.

[0057] The present invention provides a method and system for predicting crop whole-genome phenotypes by integrating environmental indicator genes. By combining resequencing data acquisition and processing, variant site data analysis, transcriptome sequencing technology, differential expression analysis, and dual-branch fusion model construction, the system achieves accurate prediction of crop phenotypes under different environmental conditions. The present invention first obtains high-quality variant site data by performing genome alignment on resequencing data, providing a solid foundation for subsequent genotype analysis; secondly, the present invention uses whole-genome association analysis technology to screen out site information significantly correlated with the target phenotype from massive variant sites, effectively reducing the data dimension and improving the targeted prediction; thirdly, the present invention uses transcriptome sequencing technology to measure the gene expression levels of crop population materials under different environments, capturing the actual impact of environmental factors on gene expression, and also screens out environmental indicator gene sets that can stably indicate environmental differences through differential expression analysis. These genes, as environmental biological indicators, can comprehensively reflect the overall impact of complex environmental factors on plant gene expression regulation; finally, a two-branch fusion phenotype prediction model is constructed based on these environmental indicator gene sets and phenotypic association site information. The genotype data and environmental indicator gene expression data are processed respectively by convolutional neural networks and fully connected neural networks, and integrated in the feature fusion layer to achieve deep learning of the interaction between genotype and environment.

[0058] The dual-branch fusion architecture constructed by the present invention using artificial intelligence algorithms can not only automatically learn complex association patterns in genotype data, but also capture the influence of environmental factors reflected by the expression levels of environmental indicator genes, thereby realizing the fusion processing of genotype and environment at the model level. It has generally broken through the limitation of traditional whole-genome selection models that rely solely on genotype information. Through the deep integration of bioinformatics and artificial intelligence technology, environmental factors are incorporated into the prediction model in a low-cost manner, effectively solving the technical problem of insufficient accuracy of phenotypic prediction under different environmental conditions.

[0059] Overall, the present invention significantly improves the accuracy of phenotypic prediction compared with traditional methods, provides a more accurate phenotypic prediction tool for crop breeding, can accelerate the crop variety selection process, and promote the sustainable development of agricultural production. It has important theoretical value and broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] In order to more clearly illustrate the technical solutions in the present invention or the prior art, a brief introduction will be given below to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0061] Figure 1 A schematic flow chart of a method for predicting crop genome-wide phenotypes by integrating environmental indicator genes provided by the present invention;

[0062] Figure 2 A schematic diagram of the saturation analysis results of the present invention for predicting environmental accuracy using different gene numbers;

[0063] Figure 3 This is a schematic diagram of the accuracy evaluation results of the phenotypic prediction model for the normal environmental phenotype without environmental indicator genes provided by the present invention;

[0064] Figure 4 This is a schematic diagram of the accuracy evaluation results of the phenotypic prediction model for low nitrogen environmental phenotypes without environmental indicator genes provided by the present invention;

[0065] Figure 5 This is a schematic diagram of the accuracy evaluation results of the phenotypic prediction model fused with environmental indicator genes provided by the present invention;

[0066] Figure 6 This is a schematic diagram of the structure of a crop genome-wide phenotype prediction system that integrates environmental indicator genes provided by the present invention. DETAILED DESCRIPTION

[0067] In order to make the purpose, technical solutions and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below in conjunction with the drawings in the present invention. Obviously, the embodiments described are part of the embodiments of the present invention, not all of the embodiments, and they should not be understood as limitations on the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. In the description of the present invention, it should be understood that the terms used are only for descriptive purposes and cannot be understood as indicating or implying relative importance.

[0068] like Figure 1 As shown, the present invention provides a crop genome-wide phenotype prediction method integrating environmental indicator genes, comprising:

[0069] S1: Collect resequencing data and perform genome alignment on the resequencing data to obtain variant site data.

[0070] Wherein, step S1 further includes:

[0071] S11: Collect resequencing data of crop population materials to obtain original sequencing data.

[0072] Raw sequencing data is obtained by whole-genome resequencing of crop populations using high-throughput sequencing technology. This data contains sequence information for every base position in the genome of an individual crop, stored in FASTQ format. Each sequence read contains a base sequence and a corresponding quality score. Specifically, the present invention first downloaded high-depth resequencing data from 208 natural populations from published literature to obtain the raw sequencing data.

[0073] S12: Performing quality control on the original sequencing data to obtain quality control sequencing data.

[0074] Among them, step S12 specifically includes: discarding reads whose base quality values ​​at the beginning of the reads are lower than a preset threshold; discarding reads whose base quality values ​​at the end of the reads are lower than a preset threshold; and discarding reads whose length after clipping is lower than a preset threshold.

[0075] In step S12, the present invention uses the Trimmomatic program to perform quality control on the original resequencing data. The specific parameters are LEADING: 3 (reads with a base quality value lower than the threshold are discarded); TRAILING: 3 (reads with a base quality value lower than the threshold are discarded); SLIDINGWINDOW: 4:15 (sliding window size starting from the 5' end of the reads); MINLEN: 36 (reads with a length less than the threshold after trimming are discarded). Finally, the filtered data is obtained for subsequent analysis.

[0076] S13: Align the quality control sequencing data to the reference genome to obtain variant site data.

[0077] Wherein, step S13 further includes:

[0078] S131: The quality control sequencing data were aligned to the reference genome using the BWA-MEM algorithm to obtain preliminary alignment results.

[0079] Furthermore, the BWA-MEM algorithm is a sequence alignment algorithm based on suffix arrays and FM indexing, designed to quickly and accurately locate short sequence reads to corresponding locations on a reference genome. The algorithm first constructs an FM index for the reference genome. The FM index is a compressed full-text index structure that enables both exact and approximate string matching in linear time. The BWA-MEM algorithm's workflow begins with seed expansion. For each quality control sequencing read, a seed sequence is searched for that exact match to the reference genome. The seed length is typically set to 19 bases. The seed region is then bidirectionally expanded using the Smith-Waterman dynamic programming algorithm to find the optimal local alignment. The algorithm calculates the alignment score by weighting matches, mismatches, and indels differently, assigning positive scores to matches and negative scores to mismatches and indels. Dynamic programming is used to find the highest-scoring alignment path. The BWA-MEM algorithm employs a divide-and-conquer strategy for long reads. When a read cannot be fully aligned, the algorithm splits the read into multiple subsegments for separate alignments, then connects the subsegments to form a complete alignment using chained alignment. The preliminary alignment results are stored in SAM format, including the alignment position, alignment quality, mismatch information and insertion and deletion information of each read segment. The alignment quality value reflects the reliability of the alignment position, the mismatch information records the base position where the read segment is inconsistent with the reference genome, and the insertion and deletion information marks the insertion and deletion variations of the sequence.

[0080] S132: Screening the preliminary comparison results according to the comparison quality and distance standards to obtain valid comparison results.

[0081] Furthermore, alignment quality screening is based on the MAPQ value, a log-likelihood ratio representation of alignment quality that reflects the probability that the alignment position is correct. During the screening process, a MAPQ threshold is set, and only alignments exceeding the threshold are retained. At the same time, the alignment must be a unique optimal match, meaning that the read has only one optimal alignment position on the reference genome. Distance-based screening is performed on paired-end sequencing data. Paired-end sequencing generates two reads from opposite ends of the same DNA fragment. The distance between the alignment positions of the two reads on the reference genome should fall within the expected range of DNA fragment length. Duplicate reads are identified and removed by comparing the alignment positions and sequence information of the reads. The final valid alignment is high-quality alignment data retained after rigorous quality and distance screening. The filtered data has accurate genomic location information and reliable alignment quality.

[0082] S133: Perform variation detection on the valid comparison results to obtain a single sample variation identification file of a single material.

[0083] The variant detection process is implemented using the HaplotypeCaller function of the GATK software package. This tool identifies single nucleotide polymorphisms and indel variants based on local genotype reassembly technology. The variant quality score is represented by a QUAL value, with higher QUAL values ​​indicating higher confidence in the variant. The final single-sample variant identification file is stored in the GVCF format. The GVCF format not only records the information of the detected variant site but also retains the genotype probability distribution at each genomic position, including the confidence level of the reference genotype and variant genotype. The GVCF file contains detailed annotations such as chromosome position, reference base, variant base, genotype quality, depth information, and allele frequency, providing complete single-sample variant information for subsequent population variation analysis.

[0084] S134: Merge the single sample variation identification files of all materials to obtain a population variation information file.

[0085] The population variation information merging process was implemented using Glnexus software, which is used for population genotype merging and quality control of large-scale genomic data. The final population variation information file is stored in VCF format, which contains standardized variant site coordinates, reference and variant allele sequences, genotype information for each sample, and associated quality annotations. Each line of the VCF file represents a variant site and contains the genotype distribution of that site across all 208 wheat accessions. Genotypes are represented by 0 / 0 for the homozygous reference genotype, 0 / 1 for the heterozygous genotype, and 1 / 1 for the homozygous variant genotype. This standardized format enables subsequent population genetics analysis and variant filtering steps to efficiently process large-scale variation data.

[0086] S135: Filtering the population variation information file according to minor allele frequency and data missing rate standards to obtain variation site data.

[0087] In a specific embodiment, the present invention uses the BWA-MEM tool in the BWA software to align the sequenced reads after quality control to the reference genome IWGSCv1.0. To ensure the accuracy of the alignment, the only optimally matched read pair is subsequently retained. Specifically, bamtools is used to remove fragments with a pairing distance greater than 10,000bp or 0bp and a quality lower than 1, and samtools is used to remove PCR duplications to obtain a valid BAM format file; the HaplotypeCaller function of the GATK software package is used, based on the effective alignment results, the single nucleotide polymorphism (SNP) variation information of each material in the population is detected and identified to obtain a GVCF format file; the GVCF format files of all materials are merged using the Glnexus software, and finally a VCF format file containing 208 materials is obtained.

[0088] S2: Perform genome-wide association analysis using the variant site data to obtain phenotype-associated site information.

[0089] Wherein, step S2 further includes:

[0090] S21: Perform principal component analysis on the variant site data to obtain population structure information.

[0091] Principal component analysis (PCA) was performed using the PCA function in PLINK software. PLINK first constructed a genotype matrix from the genotype data for 48,522,876 SNPs, with each row representing a wheat accession and each column representing a SNP. A subsequent algorithm then calculated the covariance matrix of the genotype matrix, in which each element represents a measure of genotypic similarity between two accessions. The eigenvalue decomposition process mathematically decomposes the covariance matrix, yielding 208 eigenvalues ​​and corresponding eigenvectors. The eigenvalues ​​reflect the proportion of genetic variation explained by each principal component, while the eigenvectors represent the score of each accession on that principal component. Principal components are sorted from largest to smallest eigenvalue, with the first principal component explaining the greatest proportion of genetic variation, the second the second, and so on. PLINK software outputs the eigenvectors of the first few principal components as population structure information. Each material has a numerical score on each principal component. These score values ​​constitute the population structure characteristics of the material. The population structure information reflects the genetic background differences and population stratification between different materials. When the principal component scores of two materials are similar, it means that they have similar genetic backgrounds. When the scores differ greatly, it indicates the existence of population stratification.

[0092] S22: Calculate a kinship matrix based on the group structure information to obtain group kinship information.

[0093] The kinship matrix was calculated using GEMMA software, and a genetic correlation matrix between accessions was constructed based on variance-covariance analysis. GEMMA software uses a standardized kinship matrix algorithm. The genotype data for each SNP locus are first centered. The centering process involves subtracting the mean genotype value for each locus from the genotype value across all accessions. The standardization process normalizes the variance of each SNP locus, calculates the genotype variance for each locus, and then divides the centered genotype value by the standard deviation of that locus. The resulting standardized kinship matrix has each element representing the genetic correlation coefficient between accessions. The diagonal elements of the matrix represent the genetic correlation between each accession and itself, which is close to 1. The off-diagonal elements represent the degree of relatedness between different accessions, with larger values ​​indicating closer relatedness. Population relatedness information, which contains a measure of genetic similarity between each pair of accessions, is used to control for the confounding effects of family structure and recessive population stratification in genome-wide association studies.

[0094] S23: performing association analysis on the phenotypic data and the variant site data based on the population structure information and the population kinship information to obtain phenotype-associated site information.

[0095] Genome-wide association analysis was performed using a mixed linear model using GEMMA software, which considers the contributions of both fixed and random effects to phenotypic variation. The resulting phenotypic association loci include the chromosomal location of the significant SNPs, allele information, effect size, significance level, and the proportion of phenotypic variation explained. This information forms the genotype input features for subsequent phenotypic prediction models.

[0096] In a specific embodiment, the ear length phenotype is taken as an example. First, a one-year field experiment is conducted, and then the ear length of the experimental materials at maturity in two environments (normal plots and low nitrogen plots) is measured. The ear length is measured from the first spikelet below the ear to the top of the ear. A total of three replicates are set, and five ears of each material are measured in each replicate. Finally, the BLUP (Best Linear Unbiased Prediction) value of each material is calculated and used as the final phenotype.

[0097] The obtained VCF files were subsequently filtered using PLINK software for minor allele frequency (MAF > 0.05) and missing data rate (<10%), retaining only biallelic SNPs. The filtered dataset contained 48,522,876 high-quality variants. To correct for population structure bias, principal component analysis (PCA) was performed on the genotype dataset using PLINK software, and the first five principal components were used to estimate population structure. GEMMA software was also used to calculate the variance-covariance kinship matrix as a correction for kinship.

[0098] Finally, genome-wide association analysis was performed using GEMMA software for the spike length phenotypes under the two environments. The significance of the association was ≤1*10 -6 The variant sites were considered to be significantly associated with ear length. The associated sites obtained under the two environments were merged to form a set of associated sites (genotype data).

[0099] S3: Based on transcriptome sequencing technology, gene expression levels of samples of crop population materials under different environments are measured to obtain gene expression data.

[0100] Wherein, step S3 further includes:

[0101] S31: Sampling the same crop population materials in different environments and extracting RNA to obtain sample RNA.

[0102] Furthermore, the sampling site of the present invention selects the flag leaf of the main stem of the plant. The flag leaf is the uppermost leaf of the wheat plant, which is directly involved in the nutrient supply and photosynthesis during the development of the ear. Its gene expression status can reflect the response of the plant to environmental conditions. Specifically, the present invention samples 208 materials in normal environment and low nitrogen environment respectively. About 100 mg of fresh flag leaf tissue is collected for each material, and it is immediately put into liquid nitrogen for quick freezing to prevent RNA molecule breakage caused by RNA degrading enzyme activity. RNA extraction is achieved by the Trizol method. The Trizol reagent can dissolve the cell membrane and cell wall at the same time, release the nucleic acid molecules in the cell, and the RNA finally obtained is the purified sample RNA.

[0103] S32: Performing high-throughput sequencing on the sample RNA to obtain raw transcriptome sequencing data.

[0104] Furthermore, the high-throughput transcriptome sequencing is achieved using RNA-seq technology, with a sequencing depth of 6Gbp per sample, which can cover the detection needs of most expressed genes in the wheat genome. The sequencing reaction adopts the principle of sequencing by synthesis. DNA polymerase adds fluorescently labeled deoxynucleotides one by one. After each addition reaction, the fluorescent signal is detected by laser excitation, and the type of added base is identified based on the fluorescent color. The raw transcriptome sequencing data is stored in FASTQ format, and each sample generates tens of millions of sequencing reads of 150 nucleotides in length. Each read represents the sequence information of a cDNA fragment and the corresponding sequencing quality score.

[0105] S33: Perform quality control on the original transcriptome sequencing data to obtain transcriptome quality control sequencing data.

[0106] Specifically, the present invention uses the Trimmomatic program to perform quality control and filtering on the raw transcriptome sequencing reads, with the parameters set to "SLIDING-WINDOW:4:20, MINLEN:40." The processed transcriptome quality control sequencing data has high base accuracy, appropriate length, and clean sequence content, which can accurately reflect the level of gene transcription activity.

[0107] S34: Aligning the transcriptome quality control sequencing data to a reference genome, and extracting unique aligned reads.

[0108] Furthermore, in step S34, the filtered high-quality reads were aligned to the reference genome (IWGSCv1.0) using the default parameters of the STAR software, and the bam file generated after alignment was used by the SAMTOOLS software to extract the unique aligned reads for subsequent calculation of gene expression.

[0109] S35: Calculating expression data for the uniquely aligned reads to obtain gene expression data.

[0110] Gene expression was calculated using featureCounts software, which counts the number of reads mapping to each gene region based on gene annotation information. The featureCounts software reads the reference genome's gene annotation file, which contains information on each gene's chromosomal location, exon coordinates, transcript structure, and gene function. The wheat reference genome (IWGSCv1.0) contains 107,891 annotated genes. Read counting uses a joint mode. When a read's mapping position spans multiple exons, the algorithm assigns the read to the gene containing all mapped exons. When a read's mapping position overlaps with multiple genes, the algorithm uses the maximum overlap principle to determine gene affiliation. The counting strategy is set to a gene-level approach, treating all transcripts of the same gene as a single unit. All reads mapping to any transcript of that gene are counted towards that gene's expression. Raw count data represent the number of reads per gene in each sample. This value directly reflects the level of transcriptional activity, but requires normalization to facilitate comparison across samples. The TPM normalization method first calculates the effective length of each gene. The effective length is equal to the total exon length of the gene minus the sequencing read length plus 1. This calculation compensates for the effect of gene length on read counts. The RPK value (RPK) is calculated by dividing the raw counts by the effective length. RPK represents the number of reads per kilobase, eliminating the influence of gene length differences on expression comparisons. The TPM value is calculated by dividing the RPK value by the sum of the RPK values ​​of all genes and multiplying by 1,000,000. TPM normalization ensures that the sum of all gene expression values ​​in each sample is equal to 1,000,000, eliminating the influence of sequencing depth differences on expression comparisons. Differential expression analysis was performed using the DESeq2 software package. DESeq2 performs statistical tests on raw count data based on the negative binomial distribution model, which effectively models the dispersion and overdispersion of RNA-seq count data. The DESeq2 algorithm first estimates the dispersion parameter for each gene, which reflects the variability in gene expression between biological replicates. Then, a generalized linear model is used to fit gene expression differences under different conditions. Differentially expressed genes were identified based on an expression fold change greater than 2 and a corrected P value less than 0.05. Fold changes were calculated using log2-transformed expression differences, and P values ​​were corrected using the Benjamini-Hochberg method to control the false discovery rate. The resulting gene expression data included TPM-normalized expression levels and differential expression statistics for each gene in each sample. These data formed the foundational dataset for subsequent environmental indicator gene screening.

[0111] S4: performing differential expression analysis on the gene expression data to screen environmental indicator genes and obtain an environmental indicator gene set.

[0112] Wherein, step S4 further includes:

[0113] S41: Perform differential expression analysis on gene expression data to obtain differentially expressed gene sets between different environments.

[0114] In step S41, the gene expression data are first analyzed by statistical methods to screen out genes with significant differential expression. The standard adopted in the present invention is |log2(FoldChange)|>1 and the corrected p-value <0.05. The purpose is to distinguish between genes with significant and non-significant differences and provide a basic data set for subsequent analysis.

[0115] S42: Screening genes with significant differences from the differentially expressed gene set to obtain candidate environmental indicator genes.

[0116] Specifically, a stability analysis was conducted on the differential expression characteristics of genes across the dual-environment system in all experimental materials. A threshold of 80% was used to preliminarily screen genes that were stably expressed significantly high or low in one environment. The top 30 genes were selected as candidate environmental indicator genes. The screened genes were able to combine the significance of differences with the consistency of expression across environments to reduce noise interference.

[0117] S43: Analyze the expression stability of the candidate environmental indicator genes in different environments, and screen genes with stable differential expression using a preset threshold to obtain a stable differential expression gene set.

[0118] In step S43, the present invention evaluates the expression stability of candidate genes across environments. The evaluation parameters include the coefficient of variation and the cross-environmental expression volatility. A set of genes with stable differential expression is screened out through a preset threshold, aiming to enhance the robustness of gene response to the environment and ensure the reliability of the indicator gene.

[0119] S44: constructing an environmental prediction model based on the XGBoost algorithm, ranking the importance of the stable differentially expressed gene set, screening the stable differentially expressed genes according to the ranking results and the model prediction accuracy, and obtaining an environmental indicator gene set.

[0120] In step S44, the present invention constructs an environmental prediction model based on XGBoost, uses the expression levels of 30 candidate environmental indicator genes in each sample as input, predicts the environment in which the sample is located, and thus obtains the importance ranking of the candidate indicator genes. Finally, based on the ranking, the number of genes used in the environmental prediction model is gradually increased to obtain an environmental indicator gene set.

[0121] The present invention uses different gene numbers to predict the accuracy of environmental saturation analysis schematic diagram as shown in Figure 2. Figure 2 The horizontal axis is the number of genes, and the vertical axis is the accuracy. Figure 2As can be seen from the results, the accuracy of environmental prediction using the top 8 genes reached the highest (98.6%), so these 8 genes were ultimately used as environmental indicator genes in subsequent models.

[0122] S5: constructing a dual-branch fusion phenotype prediction model using the environmental indicator gene set, the gene expression data, and the phenotype-associated site information.

[0123] Wherein, step S5 further includes:

[0124] S51: performing structural conversion on the phenotype-associated site information according to a preset encoding method to obtain genotype encoding data.

[0125] In step S51, the present invention converts discrete genotype information (SNP sites) into a numerical feature matrix suitable for deep learning model processing. Genotype encoding needs to consider the balance between genetic nature and computational efficiency, so the present invention selects a coding method suitable for dominant / recessive genetic pattern analysis, binary coding (0 / 1 / 2), to encode the SNP information of each material.

[0126] S52: normalizing the expression levels of multiple genes in the environmental indicator gene set to obtain environmental indicator feature data.

[0127] In step S52, the present invention adopts FPKM (Fragments Per Kilobase Million) or TPM

[0128] (Transcripts Per Million) eliminates differences in sequencing depth, and uses the ComBat algorithm or SVA (Surrogate Variable Analysis) to eliminate systematic errors introduced by experimental operations and instrument fluctuations. In addition, log2(x+1) transformation is performed on highly expressed genes to avoid interference of extreme values ​​on neural network training, thereby achieving the goal of eliminating differences in environmental gene expression across samples and experimental batches.

[0129] S53: using a convolutional neural network to learn the genotype encoding data to obtain a genotype feature representation;

[0130] S54: Learning the environmental indicator feature data through a fully connected neural network to obtain an environmental feature representation;

[0131] S55: Based on the feature fusion layer, the genotype feature representation and the environmental feature representation are integrated to obtain a fused feature representation.

[0132] In steps S53 to S55, the present invention first captures the local linkage disequilibrium (LD) and epistasis effects between SNP sites through CNN to obtain genotype feature representation, then analyzes the nonlinear interaction relationship between the expression levels of environmental indicator genes, quantifies the intensity of environmental stress, and finally integrates the complementary information of genotype and environmental characteristics through a feature fusion layer.

[0133] That is, the model contains two parts of input content: genotype: for the spike length-associated sites obtained through GWAS in the previous step, the SNP information of each material is encoded using the "0-1-2" method; expression level: for the 8 environmental indicator genes obtained in the previous step, the gene expression levels of each material under two environments are standardized as input to indicate the environment in which the current sample is located.

[0134] S56: Outputting the phenotype prediction result through the joint decision layer according to the fusion feature representation to obtain a phenotype prediction model.

[0135] The final model adopts a dual-branch fusion architecture to process genotype data and expression data separately, and outputs the prediction results through a joint decision layer. The dual-branch fusion architecture is divided into: a convolutional neural network branch for genotype; a fully connected network branch for expression; and a feature fusion layer that integrates the features of the two branches.

[0136] S6: Predicting the material to be predicted using the phenotypic prediction model to obtain phenotypic prediction results for different environments.

[0137] The phenotype prediction result in step S6 is a predicted phenotype value, which represents the phenotypic quantitative value of the material to be predicted under the target environment.

[0138] In addition to the above-mentioned predicted phenotypic values, the present invention compares the predicted phenotypic values ​​of each material under two environments output in each cross-validation with the actual phenotypic values, and evaluates the accuracy of the model by calculating the Pearson correlation coefficient.

[0139] Evaluation results such as Figures 3 to 5 As shown, Figure 3 This is a schematic diagram of the evaluation results of the normal environment phenotype. Figure 4 Schematic diagram of the evaluation results of low nitrogen environment phenotype. Figure 5 This is a schematic diagram of prediction accuracy after adding environmental indicator genes. Figures 3 to 5 In the figure, the horizontal axis is the actual ear length, and the vertical axis is the predicted ear length, in cm. Figure 3 and Figure 4 It can be seen that the prediction accuracy without using environmental indicator genes in the two environments is 0.68 and 0.71 respectively. Figure 5It can be seen that the prediction accuracy after adding environmental indicator genes is 0.79, indicating that the integration of environmental indicator genes can significantly improve the accuracy of phenotypic prediction.

[0140] like Figure 6 As shown, the present invention also provides a crop genome-wide phenotype prediction system integrating environmental indicator genes, comprising:

[0141] Alignment module 100: used to collect resequencing data and perform genome alignment on the resequencing data to obtain variant site data;

[0142] Association analysis module 200: used to perform genome-wide association analysis on the variant site data to obtain phenotype-associated site information;

[0143] Determination module 300: used to measure gene expression levels of samples of crop population materials under different environments based on transcriptome sequencing technology to obtain gene expression level data;

[0144] Screening module 400: for performing differential expression analysis on the gene expression data to screen environmental indicator genes and obtain an environmental indicator gene set;

[0145] Construction module 500: for constructing a double-branch fusion phenotype prediction model using the environmental indicator gene set, the gene expression data, and the phenotype-associated site information;

[0146] The prediction module 600 is configured as the phenotypic prediction model obtained by training the construction module 500, and is used to predict the material to be predicted and obtain phenotypic prediction results for different environments.

[0147] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, i.e., they may be located in one location or distributed across multiple network units. Some or all of the modules may be selected based on actual needs to achieve the objectives of the present embodiment. Persons of ordinary skill in the art will be able to understand and implement the present invention without inventive effort.

[0148] Through the description of the above embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus a necessary general hardware platform, or of course, by hardware. Based on this understanding, the essence of the above technical solution or the part that contributes to the existing technology can be embodied in the form of a software product. The computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a magnetic disk, an optical disk, etc., and includes a number of instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or certain parts of the embodiments.

[0149] The present invention provides a whole-genome phenotype prediction method and system for crops that integrate environmental indicator genes. The method obtains resequencing data, performs quality control and alignment of the sequencing data, uses genome-wide association analysis (GWAS) to obtain sites significantly associated with the phenotype, obtains transcriptome data of materials under different environments and calculates gene expression levels, screens environmental indicator genes, constructs a phenotype prediction model, and evaluates the phenotype prediction model. Ultimately, a whole-genome selection model that incorporates overall environmental factors is constructed in a low-cost manner, thereby improving the accuracy of model predictions under different environments.

[0150] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the various embodiments of the present invention.

Claims

1. A method for predicting crop genome-wide phenotypes by integrating environmental indicator genes, characterized in that: include: S1: Collect resequencing data and perform genome alignment on the resequencing data to obtain variant site data; S2: Perform genome-wide association analysis using the variant site data to obtain phenotype-associated site information; S3: Gene expression data were obtained by measuring the gene expression of samples of crop population materials under different environments based on transcriptome sequencing technology; S4: performing differential expression analysis on the gene expression data to screen environmental indicator genes and obtain an environmental indicator gene set; S5: constructing a double-branch fusion phenotype prediction model using the environmental indicator gene set, the gene expression data, and the phenotype-associated site information; S6: Predicting the material to be predicted using the phenotypic prediction model to obtain phenotypic prediction results for different environments.

2. The method for crop genome-wide phenotype prediction integrating environmental indicator genes according to claim 1, characterized in that: Step S1 further comprises: S11: Collect resequencing data of crop population materials to obtain original sequencing data; S12: performing quality control on the original sequencing data to obtain quality control sequencing data; S13: Align the quality control sequencing data to the reference genome to obtain variant site data.

3. The method for crop genome-wide phenotype prediction integrating environmental indicator genes according to claim 2, characterized in that: Step S12 specifically includes: Discard reads whose base quality values ​​at the beginning of the reads are lower than the preset threshold; Discard reads whose base quality values ​​at the end of the reads are lower than the preset threshold; After clipping, reads whose length is less than a preset threshold are discarded.

4. The method for predicting crop genome-wide phenotypes by integrating environmental indicator genes according to claim 2, characterized in that: Step S13 further includes: S131: Align the quality control sequencing data to the reference genome using the BWA-MEM algorithm to obtain preliminary alignment results; S132: Screening the preliminary comparison results according to comparison quality and distance standards to obtain valid comparison results; S133: Perform variation detection on the valid comparison results to obtain a single sample variation identification file for a single material; S134: Merge the single sample variation identification files of all materials to obtain a population variation information file; S135: Filtering the population variation information file according to minor allele frequency and data missing rate standards to obtain variation site data.

5. The method for crop genome-wide phenotype prediction integrating environmental indicator genes according to claim 1, characterized in that: Step S2 further comprises: S21: performing principal component analysis on the variant site data to obtain population structure information; S22: Calculating a kinship matrix based on the group structure information to obtain group kinship information; S23: performing association analysis on the phenotypic data and the variant site data based on the population structure information and the population kinship information to obtain phenotype-associated site information.

6. The method for crop genome-wide phenotype prediction integrating environmental indicator genes according to claim 1, characterized in that: Step S3 further comprises: S31: Sampling the same crop population materials in different environments and extracting RNA to obtain sample RNA; S32: performing high-throughput sequencing on the sample RNA to obtain raw transcriptome sequencing data; S33: Perform quality control on the original transcriptome sequencing data to obtain transcriptome quality control sequencing data; S34: aligning the transcriptome quality control sequencing data to a reference genome, and extracting unique aligned reads; S35: Calculating expression data for the uniquely aligned reads to obtain gene expression data.

7. The method for crop genome-wide phenotype prediction integrating environmental indicator genes according to claim 1, characterized in that: Step S4 further comprises: S41: Perform differential expression analysis on gene expression data to obtain differentially expressed gene sets between different environments; S42: Screening genes with significant differences from the differentially expressed gene set to obtain candidate environmental indicator genes; S43: analyzing the expression stability of the candidate environmental indicator genes in different environments, screening genes with stable differential expression using a preset threshold, and obtaining a stable differential expression gene set; S44: constructing an environmental prediction model based on the XGBoost algorithm, ranking the stable differentially expressed gene set by importance, screening the stable differentially expressed genes according to the ranking results and the model prediction accuracy, and obtaining an environmental indicator gene set.

8. The method for crop genome-wide phenotype prediction integrating environmental indicator genes according to claim 1, characterized in that: Step S5 further comprises: S51: performing structural conversion on the phenotype-associated site information according to a preset encoding method to obtain genotype encoding data; S52: normalizing the expression levels of multiple genes in the environmental indicator gene set to obtain environmental indicator characteristic data; S53: using a convolutional neural network to learn the genotype encoding data to obtain a genotype feature representation; S54: Learning the environmental indicator feature data through a fully connected neural network to obtain an environmental feature representation; S55: Based on the feature fusion layer, the genotype feature representation and the environmental feature representation are integrated to obtain a fused feature representation; S56: Outputting the phenotype prediction result through the joint decision layer according to the fusion feature representation to obtain a phenotype prediction model.

9. The method for crop genome-wide phenotype prediction integrating environmental indicator genes according to claim 1, characterized in that: The phenotype prediction result in step S6 is a predicted phenotype value, which represents the phenotypic quantitative value of the material to be predicted under the target environment.

10. A crop genome-wide phenotype prediction system integrating environmental indicator genes, characterized in that: include: Alignment module: used to collect resequencing data and perform genome alignment on the resequencing data to obtain variant site data; Association analysis module: used to perform genome-wide association analysis using the variant site data to obtain phenotype-associated site information; Measurement module: used to measure gene expression levels of samples of crop population materials under different environments based on transcriptome sequencing technology to obtain gene expression data; Screening module: used to perform differential expression analysis on the gene expression data to screen environmental indicator genes and obtain an environmental indicator gene set; Construction module: for constructing a double-branch fusion phenotype prediction model using the environmental indicator gene set, the gene expression data, and the phenotype-associated site information; The prediction module is configured as the phenotypic prediction model obtained by training the building module, and is used to predict the material to be predicted and obtain phenotypic prediction results for different environments.

Images