Method for identifying genotype-phenotype associated mutation sites of space mutation germplasm resources

Abstract

The application discloses a spaceflight mutagenic germplasm resource genotype-phenotype correlation variation site identification method, and belongs to the technical field of bioinformatics. The environment parameters and germplasm phenotype data collected by a field sensor and a detector in real time are dynamically mapped to a genotype-phenotype correlation physical model, and a simulation platform is constructed in combination with historical observation data; genotype variation and corresponding phenotype response data under different mutagenic treatment conditions are generated through Monte Carlo simulation; a particle swarm optimization algorithm is used to formulate a correlation site screening strategy; parameter configuration in an identification process is adaptively adjusted to maximize variation site identification efficiency; the detection effect of the strategy is hierarchically evaluated; and the optimal mutagenic condition and variation site combination are screened. The application solves the problem that a fixed threshold screening strategy is difficult to balance between positive detection and false positive control, and significantly enhances the repeatability and stability of the identified variation sites in a real field environment.

Description

Technical Field

[0001] This invention discloses a method for identifying genotype-phenotype associated variation sites in space-induced mutagenesis germplasm resources, belonging to the field of bioinformatics technology. Background Technology

[0002] In space-induced mutation breeding research, treating crop germplasm with mutagenic factors such as high-energy particle radiation and microgravity in the space environment can induce random variations across the genome. Screening for favorable variant sites associated with superior phenotypes is a core step in breeding work. Currently, genotype-phenotype association analysis typically relies on static genome sequencing data under laboratory conditions and limited field phenotypic observation records, using conventional statistical association analysis methods for site screening. However, field environmental factors such as temperature, humidity, light, and soil conditions exhibit significant spatiotemporal dynamic variations. These environmental fluctuations strongly influence the phenotypic expression effects of variant sites, and existing methods fail to incorporate real-time environmental parameters into the association analysis framework, resulting in insufficient reproducibility and explanatory power of the identified variant sites under actual field conditions.

[0003] In existing computer-aided breeding technologies, some studies attempt to simulate mutagenesis processes using simulation platforms. For example, Monte Carlo methods are used to generate radiation-induced genomic damage distributions, or crop growth models are constructed to predict phenotypic changes. However, these simulation platforms often treat mutagenesis doses and treatment conditions as fixed input parameters, failing to dynamically couple them with environmental data collected in real time by field sensors. This fails to reflect the modulating effect of different environmental backgrounds on the amplitude of variant expression. On the other hand, for screening strategies for associated loci, existing methods typically rely on fixed significance thresholds or empirical filtering rules, lacking a self-optimization mechanism to dynamically adjust screening parameters based on data characteristics. When faced with massive amounts of simulated variant data and diverse combinations of treatment conditions, fixed threshold strategies struggle to balance positive detection rates with false positive control, often resulting in the omission of key variant loci or false positives of numerous irrelevant loci, limiting further improvements in the validity of variant loci identification.

[0004] Furthermore, existing technologies rely on a simplistic approach to evaluating candidate mutation sites and mutagenesis conditions. They typically employ only global indicators such as overall detection rate or accuracy for a one-time evaluation, failing to consider the stratified impact of different mutagenesis dose ranges on identification effectiveness. They also lack comprehensive consideration of the stability of sites across multiple simulations and the consistency of site effects with environmental changes. This deficiency directly leads to a lack of reliable basis for determining the optimal combination of mutagenesis conditions and mutation sites, making it difficult to stably reproduce the final association analysis results in real field environments. Moreover, existing methods terminate the process after outputting the identification results, lacking a feedback mechanism for verifying and correcting the results using subsequent experimental data. This prevents dynamic adjustment of the optimal combination based on changes in environmental conditions, hindering the continuous improvement of genotype-phenotype association analysis accuracy. To address these issues, this invention proposes a systematic method that combines real-time field sensor data, simulation platform generation, particle swarm optimization screening, stratified effect evaluation, and post-validation correction to effectively overcome the shortcomings of existing technologies. Summary of the Invention

[0005] According to a first aspect of the present invention, the present invention claims protection for a method for identifying genotype-phenotype associated variation sites in space-induced mutagenesis germplasm resources, comprising: S1: Acquire real-time environmental parameters and germplasm phenotypic data from sensors and detectors deployed in the crop field, and dynamically map them to a preset genotype-phenotype association physical model, using a set of differential equations to describe the dynamic coupling relationship between environmental factors and germplasm phenotype. S2, Analyze the real-time environmental parameters, the germplasm phenotypic data and the acquired historical germplasm observation data, construct a simulation platform with input interface for mutagenesis condition parameters and output interface for mutation sites, and use Monte Carlo simulation method to generate genotypic variation and corresponding phenotypic response data under different combinations of mutagenesis doses and treatment conditions; S3. A particle swarm optimization algorithm is used to formulate an associated site screening strategy. The genotype variation and corresponding phenotypic response data output by the simulation platform are used as input to adjust the parameter configuration in the identification process in order to maximize the validity of the variant site identification. S4. After executing the associated site screening strategy, evaluate the detection effect of the strategy on candidate variant sites. The detection effect is measured by comparing the degree of agreement between the preset simulated truth set in the simulation platform and the actual output set of the strategy. S5. Based on the evaluation results of the detection effect, screen and determine the optimal combination of mutagenesis conditions and variant sites to achieve the highest genotype-phenotype association resolution accuracy of crop germplasm under the current environment and genetic background.

[0006] Furthermore, in S1, the real-time environmental parameters are collected wirelessly and synchronously by multiple sets of sensor nodes deployed in different ecological locations in the field. The sensor nodes output environmental temperature, environmental humidity, light intensity, wind speed and direction and soil conductivity parameters at preset time intervals. The detector outputs crop canopy spectral reflectance and leaf surface temperature distribution parameters at preset spatial intervals. The real-time environmental parameters and the germplasm phenotypic data are time-stamp aligned and spatial coordinate registered before being entered into the physical model.

[0007] Further, in S1, the genotype-phenotype association physical model includes an environmental response module and a genetic expression module. The environmental response module receives the real-time environmental parameters and outputs an environmental stress vector. The genetic expression module receives the environmental stress vector and outputs a phenotypic plasticity interval. The dynamic mapping substitutes the real-time environmental parameters at each sampling time into the environmental response module to solve for the transient environmental stress value. Then, the transient environmental stress value is used as a boundary condition to input into the genetic expression module to update the expected phenotypic value at the current time. The residual between the expected phenotypic value and the germplasm phenotypic data measured by the detector is recorded in real time and used to drive the parameter fine-tuning of the physical model.

[0008] Furthermore, in S2, the simulation platform constructs a digital twin field environment based on the real-time environmental parameters, the germplasm phenotypic data, and the historical germplasm observation data. The digital twin field environment includes a soil profile stratification model, a microclimate boundary layer model, and a crop growth kinetics model. The soil profile stratification model stores water transport parameters and nutrient diffusion parameters at different soil depths. The microclimate boundary layer model updates the near-surface radiation balance and turbulent exchange coefficients in a manner driven by daily variation cycles and weather events. The crop growth kinetics model outputs the crop's organ formation rate and biomass allocation ratio with accumulated temperature accumulation and light energy utilization efficiency as core state variables. The Monte Carlo simulation method is implemented as follows: The probability distribution functions of the mutagenic dose and the treatment duration were determined based on historical mutagenic experimental records. Random samples are independently drawn from the probability distribution function of the mutagenic dose and the probability distribution function of the treatment duration to form a set of candidate treatment condition combinations; For each candidate treatment combination, a pre-defined genomic damage simulation subroutine is executed in the digital twin field environment. The genomic damage simulation subroutine generates a set of single nucleotide variant sites and a set of structural variant sites based on the randomness of radiation energy deposition and the competitive pathway of DNA double-strand break repair. The set of single nucleotide variant sites and the set of structural variant sites are input into the phenotypic conversion subroutine in the crop growth kinetic model. The phenotypic conversion subroutine calculates the quantitative perturbation value of the variant on the organ formation rate and the biomass allocation ratio based on the regulatory network position and metabolic pathway affiliation of the gene where the variant site is located. The set of single nucleotide variant sites, the set of structural variant sites, and the corresponding quantitative perturbation values ​​under each candidate treatment condition combination are recorded as the genotype variation and corresponding phenotypic response data.

[0009] Furthermore, in S2, the simulation platform also includes a field heterogeneous noise injection module, which generates a spatially continuous noise field based on the spatial variation semivariance function of the real-time environmental parameters, and superimposes the noise field onto the corresponding phenotypic response data output by the Monte Carlo simulation method. The spatial variation semivariance function is obtained by fitting repeated measurement records from the historical germplasm observation data, and the field heterogeneity noise injection module regenerates an independent noise field in each Monte Carlo simulation iteration.

[0010] Furthermore, in S3, the particle swarm optimization algorithm encodes the parameters to be adjusted in the identification process into a particle position vector, which includes a mutation site significance threshold, a phenotypic response amplitude cutoff coefficient, and a site co-occurrence frequency threshold. After the particle swarm optimization algorithm initializes the particle population, each particle runs the associated site screening strategy based on the genotype variation and corresponding phenotypic response data output by the simulation platform, and calculates the positive detection rate and false positive detection rate of the strategy on the simulation platform. The weighted harmonic mean of the positive detection rate and the false positive detection rate is used as the fitness value. The particle swarm optimization algorithm iteratively updates the optimal position of individual particles and the global optimal position of the swarm based on the fitness value, and adjusts the flight direction and step size of each particle according to the velocity update equation until the rate of change of the fitness value is lower than the preset convergence criterion or the number of iterations reaches the preset upper limit.

[0011] Furthermore, in S3, the associated site screening strategy is composed of a hierarchical filtering architecture and a fine-grained localization architecture connected in series; The hierarchical filtering architecture takes the whole genome variation sites in the genotype variation and corresponding phenotypic response data as the processing objects. According to the chromosome partitioning strategy, the whole genome is divided into several physical intervals of equal length. Then, the variation sites in each physical interval are sorted according to the phenotypic association strength and the candidate sites with the highest ranking are retained. The phenotypic association strength is determined by the ratio of the local phenotypic effect size to the global phenotypic effect size at the location of the variation site. The fine-targeting architecture receives the candidate sites output by the hierarchical filtering architecture, performs linkage disequilibrium decay analysis on all variant sites within a preset flanking length range around each candidate site, and identifies tag sites with independent genetic effects based on the law that the correlation coefficient between paired sites decays with physical distance. During the iteration process, the particle swarm optimization algorithm synchronously adjusts the granularity of the physical interval division and the sort retention ratio in the hierarchical filtering architecture, as well as the preset flanking length and the significance cutoff value of the chain imbalance decay analysis in the fine positioning architecture. The significance cutoff value is dynamically updated by the real-time changes in the fitness value.

[0012] Furthermore, in S4, the evaluation process of the detection effect is implemented in a stratified manner according to the mutagenic dose interval. The genotype variation and corresponding phenotypic response data generated by the simulation platform are divided into low-dose subset, medium-dose subset and high-dose subset according to the mutagenic dose value. Then, the detection recall and detection precision of the associated site screening strategy are calculated in each subset. The detection recall and detection precision are weighted and summarized according to the sample size of each subset to obtain the global detection validity index. The ratio between the global detection validity index and the number of known associated sites pre-labeled in the simulated truth set is used to determine whether the current strategy has reached the preset validity threshold. If the validity threshold is not reached, the process returns to S3 to re-execute the iterative optimization of the particle swarm optimization algorithm.

[0013] Furthermore, in S5, the screening process of the optimal combination of mutagenesis conditions and mutation sites is performed in sequence as a preliminary screening stage of mutagenesis conditions and a fine screening stage of mutation sites. In the preliminary screening stage of mutagenesis conditions, the detection effect corresponding to each candidate mutagenesis condition obtained in S4 is used as the sorting basis, and multiple mutagenesis conditions with the highest detection effect are selected to form a candidate condition pool. In the mutation site screening stage, within the mutation site set corresponding to the candidate condition pool, the stability of each site in the multi-round Monte Carlo simulation on the simulation platform and the consistency of the response between each site and the real-time environmental parameters are used as two-dimensional ranking criteria. The mutation site with the highest comprehensive score under the two-dimensional ranking criteria is selected as the target site in the optimal combination. At the same time, the mutagenic dose and treatment duration corresponding to the target site are recorded as the target conditions in the optimal combination.

[0014] Furthermore, after S5, a post-verification step is performed. When the real-time environmental parameters and germplasm phenotypic data continuously collected by the sensors and detectors deployed in the crop field for the next time window arrive, the post-verification step inputs the target site in the optimal combination of mutagenesis conditions and mutation sites into the genotype-phenotype association physical model, driving the physical model to forward deduce the expected phenotypic change trajectory within the next time window. The expected phenotypic change trajectory is compared with the measured germplasm phenotypic data in the next time window. If the deviation between the two exceeds the preset allowable range, the simulation platform is triggered to re-execute the Monte Carlo simulation method with the real-time environmental parameters of the next time window as new initial conditions. The particle swarm optimization algorithm is driven by the updated data obtained after re-execution to modify the optimal combination of mutation conditions and mutation sites.

[0015] This invention discloses a method for identifying genotype-phenotype association variant sites in space-induced mutagenesis germplasm resources, belonging to the field of bioinformatics. It dynamically maps environmental parameters and germplasm phenotypic data collected in real time by field sensors and detectors to a genotype-phenotype association physical model, and constructs a simulation platform using historical observation data. Monte Carlo simulations generate genotype variations and corresponding phenotypic responses under different mutagenesis treatments. A particle swarm optimization algorithm is employed to formulate an association site screening strategy, adaptively adjusting parameter configurations in the identification process to maximize the validity of variant site identification. The detection effect of the strategy is evaluated hierarchically to screen the optimal combination of mutagenesis conditions and variant sites. This invention solves the problem of balancing positive detection and false positive control in fixed-threshold screening strategies, significantly enhancing the repeatability and stability of identified variant sites in real field environments. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the process of a method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources, as claimed in this embodiment of the invention. Figure 2 This is a second flowchart of a method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources, as claimed in an embodiment of the present invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0018] The technical solutions disclosed in the various embodiments of this application are described in detail below with reference to the accompanying drawings.

[0019] According to the first embodiment of the present invention, referring to Figure 1 This invention claims protection for a method for identifying genotype-phenotype associated variation sites in space-mutated germplasm resources, comprising: S1: Acquire real-time environmental parameters and germplasm phenotypic data from sensors and detectors deployed in the crop field, and dynamically map them to a preset genotype-phenotype association physical model, using a set of differential equations to describe the dynamic coupling relationship between environmental factors and germplasm phenotype. S2, Analyze the real-time environmental parameters, the germplasm phenotypic data and the acquired historical germplasm observation data, construct a simulation platform with input interface for mutagenesis condition parameters and output interface for mutation sites, and use Monte Carlo simulation method to generate genotypic variation and corresponding phenotypic response data under different combinations of mutagenesis doses and treatment conditions; S3. A particle swarm optimization algorithm is used to formulate an associated site screening strategy. The genotype variation and corresponding phenotypic response data output by the simulation platform are used as input to adjust the parameter configuration in the identification process in order to maximize the validity of the variant site identification. S4. After executing the associated site screening strategy, evaluate the detection effect of the strategy on candidate variant sites. The detection effect is measured by comparing the degree of agreement between the preset simulated truth set in the simulation platform and the actual output set of the strategy. S5. Based on the evaluation results of the detection effect, screen and determine the optimal combination of mutagenesis conditions and variant sites to achieve the highest genotype-phenotype association resolution accuracy of crop germplasm under the current environment and genetic background.

[0020] In this embodiment, environmental parameter sensor groups and phenotypic detector groups are installed in the target crop field according to a pre-designed spatial layout. The sensor groups include air temperature sensors, air humidity sensors, photosynthetically active radiation sensors, wind speed and direction sensors, and soil multi-depth conductivity measurement elements. Each element independently outputs readings according to a set acquisition frequency. The detector group includes a multispectral imaging probe and a thermal infrared imaging probe, used to acquire the spatial distribution of the crop canopy's reflectance spectrum in different bands and the spatial distribution of leaf surface temperature. The detectors output image data according to another set acquisition frequency. All sensors and detectors are equipped with wireless transmission modules, which transmit their respective readings, along with device identification and the timestamp of the acquisition time, to the field data aggregation base station in real time.

[0021] After receiving the raw data packets from each channel, the data aggregation base station first performs a timestamp alignment operation. The base station maintains a sliding buffer with a preset time window as the unit, and pairs the sensor readings and detector readings according to the principle of the closest timestamp, so that the environmental parameters and phenotypic data at the same observation time form a complete record vector. After completing the timestamp pairing, the base station performs a spatial coordinate registration operation. The base station has pre-stored the precise spatial coordinates of each sensor and detector, and based on these coordinates, each reading is mapped to the spatial grid position in the field, so that the environmental parameters and phenotypic data correspond to each other in space.

[0022] After completing timestamp alignment and spatial coordinate registration, a dynamic mapping operation is performed to the preset genotype-phenotype association physical model. This physical model consists of an environmental response module and a genetic expression module. The environmental response module receives real-time environmental parameters at the current moment, including environmental temperature, environmental humidity, light intensity, wind speed and direction, and soil electrical conductivity, as input variables. Internally, the module calculates the environmental stress vector under the current environmental conditions based on the crop variety's response to environmental factors at different growth stages. This environmental stress vector includes water stress components, temperature stress components, and radiation stress components, and the values ​​of each component reflect the degree of limitation of the current environmental factors on crop growth.

[0023] The genetic expression module receives the aforementioned environmental stress vector as boundary conditions and, in conjunction with the allele effect values ​​stored in the genetic parameter library of the crop variety, calculates the expected expression range of each phenotypic trait under the current environmental stress. The output of the genetic expression module is an expected phenotypic value vector, covering the expected values ​​or ranges of multiple traits such as plant height, number of grains per ear, and thousand-grain weight. The germplasm phenotypic data measured by the detector is compared item by item with the expected phenotypic values ​​output by the genetic expression module, and the residuals of each trait are calculated. These residuals are recorded in real time and used to drive the weight fine-tuning of the genetic parameter library in the physical model. The fine-tuning method adjusts the weight coefficients of the corresponding gene effect values ​​according to the sign and magnitude of the residuals, so that the subsequent output of the physical model gradually approaches the measured data.

[0024] While continuously receiving real-time data, the local computing server reads historical germplasm observation data from the storage medium. The historical germplasm observation data includes environmental parameter sequences and phenotypic trait records of the same crop variety in the same plot over multiple growing seasons. A simulation platform is built based on the real-time data and historical data.

[0025] The simulation platform includes a digital twin field environment module, which consists of a soil profile stratification model, a microclimate boundary layer model, and a crop growth kinetics model. The soil profile stratification model divides the field soil into several layers according to depth, with each layer storing water transport parameters and nutrient diffusion parameters. The microclimate boundary layer model updates the near-surface radiation balance value on a daily cycle and drives the transient adjustment of the turbulent exchange coefficient based on weather events. The crop growth kinetics model uses accumulated temperature and light energy utilization efficiency as core state variables, outputting the crop's organ formation rate and biomass allocation ratio.

[0026] The Monte Carlo simulation method is implemented in the following order: First, based on historical mutagenesis experiment records, the probability distribution functions of the mutagenesis dose and treatment duration are determined. The dose distribution function is obtained by summarizing and fitting the irradiation dose records used in previous similar space-induced mutagenesis experiments, and the duration distribution function is obtained by summarizing and fitting the corresponding irradiation duration records. Second, random samples are independently drawn from the above two probability distribution functions, and the sampling is repeated a preset number of times to obtain a corresponding number of candidate treatment condition combinations. Each combination contains one dose sample and one duration sample. Third, for each group of candidate treatment condition combinations, a genome damage simulation subroutine is invoked in the digital twin field environment. This subroutine generates an initial damage site distribution based on the randomness of radiation energy deposition in the cell nucleus, and then produces the final set of variant sites, including single nucleotide variants and structural variants, based on the competitive selection of different repair pathways for DNA double-strand breaks. Fourth, the set of variant sites is input into the phenotypic conversion subroutine in the crop growth kinetics model. This subroutine calculates the quantitative perturbation value of the variant on the organ formation rate and biomass allocation ratio based on the regulatory network position and metabolic pathway affiliation of the gene containing each variant site. The fifth step is to record the set of variant sites and quantitative perturbation values ​​corresponding to each combination of candidate treatment conditions, which will serve as the genotypic variation and corresponding phenotypic response data for that combination.

[0027] Using all genotype variations and corresponding phenotypic response data output by the simulation platform as input, the particle swarm optimization algorithm is launched to formulate an associated site screening strategy. The particle swarm optimization algorithm first encodes the parameters that need to be adjusted in the identification process into particle position vectors. These parameters include the significance threshold of the variant site, the phenotypic response amplitude cutoff coefficient, and the site co-occurrence frequency threshold. The particle population is initialized, and the position vector of each particle is randomly initialized within the range of the values ​​of each parameter, and the velocity vector of each particle is randomly initialized within a preset range.

[0028] Each particle runs an associated site screening strategy based on the parameter values ​​in its position vector, which involves performing significance testing on variant sites, filtering phenotypic response amplitudes, and screening co-occurrence frequencies on all simulated data. After screening, the positive detection rate and false positive detection rate of the strategy on the simulation platform are calculated. The positive detection rate and false positive detection rate are then calculated into a weighted harmonic mean according to preset weights, which is used as the fitness value of the particle.

[0029] The particle swarm optimization algorithm iteratively updates the individual optimal position and the global optimal position of the swarm based on the fitness value of each particle, and adjusts the flight direction and step size of each particle according to the velocity update equation. The velocity update equation includes an inertial weight term, an individual cognition term, and a social sharing term. The above iterations continue until the rate of change of the fitness value is lower than the preset convergence criterion in a series of iterations, or the number of iterations reaches the preset upper limit.

[0030] After the particle swarm optimization algorithm converges, the parameter values ​​corresponding to the globally optimal particle and the candidate mutation site output set obtained by running the associated site screening strategy under this parameter configuration are recorded. The candidate mutation site output set is then compared site by site with the preset simulation ground truth set in the simulation platform. The simulation ground truth set is a list of real associated sites pre-marked during the Monte Carlo simulation. The comparison method is as follows: for each site in the candidate mutation site output set, check whether it appears in the simulation ground truth set; for each site in the simulation ground truth set, check whether it is covered by the candidate mutation site output set. The detection recall and detection precision are calculated based on the comparison results, and these two indicators are used as the evaluation results of the detection effect of this strategy on candidate mutation sites.

[0031] Based on the detection effectiveness evaluation results obtained from S4, all candidate treatment condition combinations were sorted from high to low according to their corresponding detection recall, and the top-ranked treatment condition combinations were selected to form the initial screening set. Within the initial screening set, the combinations were further sorted from high to low according to detection precision, and the top-ranked treatment condition combinations were selected. For these treatment condition combinations, their corresponding variant site sets were extracted, and the frequency of simultaneous detection of each site under these conditions was counted. The variant site with the highest simultaneous detection frequency was selected as the optimal target site, and the corresponding mutagenic dose and treatment duration were used as the optimal target conditions. Thus, the mutagenic conditions and variant site combinations that achieve the highest resolution under the current environmental and genetic background were obtained.

[0032] Furthermore, in S1, the real-time environmental parameters are collected wirelessly and synchronously by multiple sets of sensor nodes deployed in different ecological locations in the field. The sensor nodes output environmental temperature, environmental humidity, light intensity, wind speed and direction and soil conductivity parameters at preset time intervals. The detector outputs crop canopy spectral reflectance and leaf surface temperature distribution parameters at preset spatial intervals. The real-time environmental parameters and the germplasm phenotypic data are time-stamp aligned and spatial coordinate registered before being entered into the physical model.

[0033] In this embodiment, multiple sets of sensor nodes are deployed at preset intervals within the field test area of ​​the target crop. Each set of sensor nodes integrates a temperature sensing element, an air humidity sensing element, a photosynthetically active radiation sensing element, a wind speed and direction sensing element, and a soil multi-depth conductivity measurement element. Each element starts measurement simultaneously at a unified acquisition frequency, and the readings, along with the node identifier and the timestamp of the acquisition time, are transmitted to the field data collection base station via a wireless communication protocol to form an ad hoc network.

[0034] Within the same field test area, multiple sets of detector nodes were deployed at predetermined intervals. Each set of detector nodes included a multispectral imaging probe and a thermal infrared imaging probe. The multispectral imaging probe could simultaneously acquire crop canopy reflectance data across multiple spectral bands, forming a multi-band image. The thermal infrared imaging probe acquired leaf surface temperature distribution data, forming a temperature image. Each detector node performed imaging acquisition at a different frequency than the sensor nodes, but aligned with the sensor nodes' acquisition times at the top of the hour. The detector nodes transmitted image data to the field data aggregation base station via wireless communication.

[0035] The field data collection base station is equipped with a high-precision real-time clock chip. This clock chip is periodically calibrated by an external time signal to ensure that the base station time is synchronized with the standard time. After receiving the environmental parameter data packets from the sensor nodes and the image data packets from the detector nodes, the base station extracts the timestamps from them. The base station runs a timestamp alignment process. This process uses a preset time window as the matching range and pairs the sensor readings and detector readings according to their timestamps. Sensor readings and detector readings whose timestamps differ by no more than the time window are included in the recording vector of the same observation time.

[0036] After completing timestamp pairing, the base station performs spatial coordinate registration. The base station stores the precise spatial coordinates of each sensor node and each detector node. For a sensor reading in the recorded vector at a certain moment, the base station determines its position in the field grid based on the sensor's spatial coordinates. For detector image data in the same recorded vector, the base station calculates the actual ground position corresponding to each pixel in the image based on the detector's spatial coordinates and imaging field of view parameters. Then, the spatial position of the sensor reading is superimposed with the position of the corresponding pixel in the detector image to form a complete recorded vector after spatial coordinate registration. The recorded vector after timestamp alignment and spatial coordinate registration is sent to the physical model described in S1 for subsequent processing.

[0037] Further, in S1, the genotype-phenotype association physical model includes an environmental response module and a genetic expression module. The environmental response module receives the real-time environmental parameters and outputs an environmental stress vector. The genetic expression module receives the environmental stress vector and outputs a phenotypic plasticity interval. The dynamic mapping substitutes the real-time environmental parameters at each sampling time into the environmental response module to solve for the transient environmental stress value. Then, the transient environmental stress value is used as a boundary condition to input into the genetic expression module to update the expected phenotypic value at the current time. The residual between the expected phenotypic value and the germplasm phenotypic data measured by the detector is recorded in real time and used to drive the parameter fine-tuning of the physical model.

[0038] In this embodiment, the genotype-phenotype association physical model consists of two sub-modules: an environmental response module and a genetic expression module. The environmental response module stores the response functions of the crop variety to environmental factors at different growth stages. Specifically, the environmental response module receives real-time environmental parameters, including ambient temperature, humidity, light intensity, wind speed and direction, and soil electrical conductivity, as input variables. Internally, the module first calculates the water stress index: comparing the measured soil electrical conductivity with the soil electrical conductivity threshold for the crop variety. If the measured value is lower than the threshold, it is considered a water stress state. The water stress index is calculated linearly based on the ratio of the measured value to the threshold, ranging from no stress to severe stress. Next, the module calculates the temperature stress index: comparing the measured ambient temperature with the optimal growth temperature range for the crop variety. If the measured temperature is higher than the upper limit or lower than the lower limit of the range, the degree of heat or cold stress is calculated based on the extent of the exceedance. The temperature stress index also ranges from no stress to severe stress. Finally, the module calculates the radiation stress index: comparing the measured light intensity with the light saturation point for the crop variety. If the measured value is higher than the light saturation point, the degree of light inhibition is calculated based on the extent of the exceedance. The radiation stress index ranges from no stress to severe stress. The environmental response module combines the three stress indices mentioned above into an environmental stress vector for output.

[0039] The genetic expression module receives the environmental stress vector output by the environmental response module. Internally, the module stores a genetic parameter library for the crop variety, recording the allele type and corresponding phenotypic effect value for each known functional gene. The module's operating logic is as follows: First, it determines the crop's growth stage based on the current number of days in growth. Second, based on the sensitivity of each phenotypic trait to environmental stress at this growth stage, it retrieves the corresponding allele effect value from the genetic parameter library. Third, it uses each component of the environmental stress vector as a weighting coefficient, weighting it onto the corresponding allele effect value to calculate the expected expression value of each phenotypic trait under the current environmental stress. The genetic expression module calculates the expected values ​​of multiple phenotypic traits, constructing the expected phenotypic vector output.

[0040] The specific operation of dynamic mapping is as follows: At each sampling moment, real-time environmental parameters and germplasm phenotypic data obtained from field sensors and detectors are input into the physical model. The model first runs the environmental response module to calculate the environmental stress vector based on the real-time environmental parameters at the current moment. Then, the environmental stress vector is used as the boundary condition and input into the genetic expression module, which updates the expected phenotypic vector at the current moment accordingly. The expected phenotypic vector is compared item by item with the germplasm phenotypic data measured by the detector, and the residual of each phenotypic trait is calculated. These residuals are recorded in real time in the residual log file on the server. The server periodically summarizes the residual logs and adjusts the weight coefficients of the corresponding allele effect values ​​in the genetic parameter library according to the cumulative distribution of the residuals. If the residual of a certain trait is continuously positive for a period of time, the effect value weight of the corresponding gene of that trait is appropriately increased; if the residual is continuously negative, the weight is appropriately decreased. After the above dynamic mapping and parameter fine-tuning, the expected output of the physical model gradually approaches the measured data.

[0041] Furthermore, in S2, the simulation platform constructs a digital twin field environment based on the real-time environmental parameters, the germplasm phenotypic data, and the historical germplasm observation data. The digital twin field environment includes a soil profile stratification model, a microclimate boundary layer model, and a crop growth kinetics model. The soil profile stratification model stores water transport parameters and nutrient diffusion parameters at different soil depths. The microclimate boundary layer model updates the near-surface radiation balance and turbulent exchange coefficients in a manner driven by daily variation cycles and weather events. The crop growth kinetics model outputs the crop's organ formation rate and biomass allocation ratio with accumulated temperature accumulation and light energy utilization efficiency as core state variables. Reference Figure 2 The Monte Carlo simulation method is implemented as follows: The probability distribution functions of the mutagenic dose and the treatment duration were determined based on historical mutagenic experimental records. Random samples are independently drawn from the probability distribution function of the mutagenic dose and the probability distribution function of the treatment duration to form a set of candidate treatment condition combinations; For each candidate treatment combination, a pre-defined genomic damage simulation subroutine is executed in the digital twin field environment. The genomic damage simulation subroutine generates a set of single nucleotide variant sites and a set of structural variant sites based on the randomness of radiation energy deposition and the competitive pathway of DNA double-strand break repair. The set of single nucleotide variant sites and the set of structural variant sites are input into the phenotypic conversion subroutine in the crop growth kinetic model. The phenotypic conversion subroutine calculates the quantitative perturbation value of the variant on the organ formation rate and the biomass allocation ratio based on the regulatory network position and metabolic pathway affiliation of the gene where the variant site is located. The set of single nucleotide variant sites, the set of structural variant sites, and the corresponding quantitative perturbation values ​​under each candidate treatment condition combination are recorded as the genotype variation and corresponding phenotypic response data.

[0042] In this embodiment, the simulation platform constructs a digital twin field environment on a local computing server based on real-time environmental parameters, germplasm phenotypic data, and historical germplasm observation data; the digital twin environment includes three sub-modules.

[0043] The soil profile stratification model divides the field soil profile into several layers according to depth. Each layer stores water transport parameters, including saturated hydraulic conductivity, field capacity, wilting coefficient, and nutrient diffusion parameters, including nitrate nitrogen diffusion coefficient and ammonium nitrogen diffusion coefficient. The soil profile stratification model calculates the vertical transport flux of water between each layer based on the water movement equation and the diffusion flux of nutrients between each layer based on the diffusion equation.

[0044] The microclimate boundary layer model updates the near-surface radiation balance on a daily basis. It incorporates an astronomical algorithm to calculate the diurnal variation of the solar altitude angle based on the field's latitude and longitude coordinates and the current date, thereby calculating the total daily shortwave solar radiation reaching the top of the canopy. The model also uses cloud cover forecast data provided by a weather forecasting interface to attenuate shortwave solar radiation. Simultaneously, the microclimate boundary layer model calculates the longwave radiation balance, determining the net longwave radiation loss based on air temperature and surface temperature. The model also includes a subroutine for calculating the turbulent exchange coefficient, which calculates the turbulent exchange coefficients of momentum, heat, and water vapor based on wind speed gradients and atmospheric stability. When rainfall occurs, the model updates the soil surface moisture content and recalculates the water infiltration rate.

[0045] The crop growth kinetics model uses accumulated temperature and light energy utilization efficiency as core state variables. The accumulated temperature is calculated by subtracting the crop's baseline temperature from the daily average air temperature and summing the difference to the total accumulated temperature since the sowing date. The light energy utilization efficiency is calculated by multiplying the daily total photosynthetically active radiation by the crop's light energy conversion coefficient and then by an environmental stress correction factor. The crop growth kinetics model determines the crop's growth stage based on the accumulated temperature, calculates the daily newly synthesized dry matter based on the light energy utilization efficiency, and then distributes the dry matter to each organ according to the biomass allocation ratio of that growth stage, outputting the daily organ formation rate and cumulative biomass.

[0046] Detailed implementation of Monte Carlo simulation: The first step is to determine the probability distribution function of the mutagenesis dose; to summarize the irradiation dose records of previous similar space-based mutagenesis experiments stored on the local server, convert these records into standardized absorbed dose values, perform statistical analysis on all dose values, and fit the probability distribution function of the dose values; to determine the probability distribution function of the treatment duration; to summarize the irradiation duration records of the above experiments, and fit the probability distribution function of the duration values.

[0047] The second step is to independently extract random samples from the two probability distribution functions mentioned above. Using the server's random number generator, a uniformly distributed random number is first generated. This random number is then used as the input to the cumulative distribution function of the dose distribution function to solve for the corresponding dose value. Another uniformly distributed random number is then generated to solve for the corresponding duration value. The above extraction process is repeated a preset number of times to obtain the corresponding number of candidate processing condition combinations.

[0048] The third step involves executing a genome damage simulation subroutine in a digital twin field environment for each candidate treatment combination. This subroutine first calculates the number of energy deposition events in each cell nucleus at that dose. The expected number of events is obtained by dividing the product of the dose and the cell nucleus volume by the average ionization energy. The actual number of events is then generated based on a random distribution. The location of each energy deposition event on the chromosome is randomly determined according to the chromosome's physical length. Subsequently, the DNA double-strand break repair process is simulated: for each break site, a repair path is randomly selected based on the chromatin state near that site. Different repair paths introduce single nucleotide variants or small fragment insertions / deletions with different error rates, or structural variants due to template selection errors. Through the above simulation, a set of single nucleotide variant sites and a set of structural variant sites are generated for each candidate treatment combination.

[0049] The fourth step involves inputting the set of variant sites into the phenotypic transformation subroutine of the crop growth kinetics model. This subroutine first performs functional annotation on each variant site to determine which functional region of the gene it is located in. If it is located in an exon region, it further determines whether it is a missense mutation or a synonym mutation. Based on the regulatory network position and metabolic pathway affiliation of the gene containing the variant site, it retrieves the effect coefficient of the gene from the genetic parameter database. The effect coefficient is multiplied by the heterozygosity or homozygosity of the variant site to calculate the quantitative perturbation value of the variant on organ formation rate and biomass allocation ratio. This perturbation value is then transferred to the corresponding biomass allocation ratio parameter, adjusting the parameter from the baseline value to the new value.

[0050] The fifth step involves recording the dose value, duration value, set of variant sites, and corresponding quantitative perturbation value of each candidate treatment combination as genotypic variation and corresponding phenotypic response data for that combination. All data is stored in the output database of the simulation platform for subsequent steps.

[0051] Furthermore, in S2, the simulation platform also includes a field heterogeneous noise injection module, which generates a spatially continuous noise field based on the spatial variation semivariance function of the real-time environmental parameters, and superimposes the noise field onto the corresponding phenotypic response data output by the Monte Carlo simulation method. The spatial variation semivariance function is obtained by fitting repeated measurement records from the historical germplasm observation data, and the field heterogeneity noise injection module regenerates an independent noise field in each Monte Carlo simulation iteration.

[0052] In this embodiment, the simulation platform also includes a field heterogeneous noise injection module. The function of this module is to superimpose a spatially continuously distributed noise field onto the phenotypic response data output by the Monte Carlo simulation in order to simulate the phenotypic measurement fluctuations caused by spatial heterogeneity factors such as soil texture, micro-topography, and distribution of pests and diseases in the real field environment.

[0053] This module first fits a spatial variance function based on repeated measurement records from historical germplasm observation data. Historical germplasm observation data includes repeated phenotypic measurement records collected from the same experimental plot over multiple growing seasons using a multi-point sampling method. Each sampling point has multiple repeated measurements at the same time point. For each phenotypic trait, half the square of the difference between the measurement values ​​of any two sampling points is calculated as the semivariogram value for that point pair. The semivariogram values ​​of all point pairs are grouped according to the spatial distance between the point pairs, and the average semivariogram value within each distance group is calculated. Using the average semivariogram value as the ordinate and the spatial distance as the abscissa, a semivariogram function curve is fitted. This curve exhibits typical characteristics such as the nugget effect, sill value, and range variation.

[0054] In each Monte Carlo simulation iteration, the field heterogeneous noise injection module generates a spatially continuous noise field based on the aforementioned semivariogram function. The generation method is as follows: First, at each grid point of the field spatial grid, an initial random value is extracted from a random distribution with a mean of zero and a variance equal to the sill value. Then, based on the spatial correlation structure described in the semivariogram function, the initial random value is smoothed, that is, the value of each grid point is replaced by the weighted average of the values ​​of its surrounding grid points. The weights are determined by the spatial distance between grid points and the semivariogram function, with closer grid points having greater weights. After smoothing, a spatially continuous noise field image with spatial autocorrelation is obtained.

[0055] The generated noise field is superimposed on the corresponding phenotypic response data output by the Monte Carlo simulation method. Specifically, for each simulated plant at a grid point, its phenotypic response data is added to the corresponding value of the noise field at that grid point to obtain the final phenotypic response output. In this way, the phenotypic response data output by the Monte Carlo simulation exhibits spatial heterogeneity consistent with real field observations.

[0056] The field heterogeneous noise injection module regenerates an independent noise field in each round of Monte Carlo simulation iteration. That is, each of the candidate treatment condition combinations is superimposed using its independently generated noise field, thereby ensuring that the noise between different simulation iterations is independent of each other.

[0057] Furthermore, in S3, the particle swarm optimization algorithm encodes the parameters to be adjusted in the identification process into a particle position vector, which includes a mutation site significance threshold, a phenotypic response amplitude cutoff coefficient, and a site co-occurrence frequency threshold. After the particle swarm optimization algorithm initializes the particle population, each particle runs the associated site screening strategy based on the genotype variation and corresponding phenotypic response data output by the simulation platform, and calculates the positive detection rate and false positive detection rate of the strategy on the simulation platform. The weighted harmonic mean of the positive detection rate and the false positive detection rate is used as the fitness value. The particle swarm optimization algorithm iteratively updates the optimal position of individual particles and the global optimal position of the swarm based on the fitness value, and adjusts the flight direction and step size of each particle according to the velocity update equation until the rate of change of the fitness value is lower than the preset convergence criterion or the number of iterations reaches the preset upper limit.

[0058] In this embodiment, the particle swarm optimization algorithm first determines three parameters to be adjusted and their value ranges. The first parameter is the significance threshold of the mutation site, which is used to determine whether the association effect of a certain mutation site reaches a statistically significant level during the screening process. The second parameter is the phenotypic response amplitude cutoff coefficient, which is used to filter out mutation sites with too small phenotypic response amplitudes. The third parameter is the site co-occurrence frequency threshold, which is used to screen stable mutation sites that repeatedly appear in multiple rounds of simulation.

[0059] The particle swarm optimization algorithm encodes the above three parameters into a three-dimensional particle position vector; initializes the particle population with a population size set to a preset value; for each particle, the three components of its position vector are randomly initialized in a uniform distribution within their respective value ranges; the velocity vector of each particle is also a three-dimensional vector, and its three components are randomly initialized in a uniform distribution within a preset range.

[0060] Each particle executes an association site screening strategy based on three parameter values ​​in its position vector. The specific execution method of the screening strategy is as follows: For all genotype variations and corresponding phenotypic response data output by the simulation platform, firstly, the association strength statistic between each variation site and the phenotypic response is calculated, and the significance level of the statistic is compared with the significance threshold parameter of the particle, retaining sites that reach the significance level; then, the phenotypic response amplitude of the retained sites is calculated, and the amplitude is compared with the truncation coefficient parameter of the particle. First, the dispersion index of the phenotypic response amplitude of all sites is calculated, and the truncation coefficient is multiplied by the index as the truncation value, retaining sites with amplitudes greater than the truncation value; finally, the number of times each point appears in all simulation data is counted, the frequency of occurrence is calculated, and the frequency is compared with the co-occurrence frequency threshold parameter of the particle, retaining sites with frequencies higher than the threshold; after the above three-level filtering, the remaining sites are the candidate variation site output set under the particle configuration.

[0061] Calculate the strategy performance under this particle configuration; compare the candidate mutation site output set with the preset simulation truth set in the simulation platform; the simulation truth set contains a known number of real associated sites; count the number of sites in the candidate output set that appear in the simulation truth set, and calculate the positive detection rate; count the number of sites in the candidate output set that do not appear in the simulation truth set, and calculate the false positive detection rate; calculate the weighted harmonic average of the positive detection rate and the false positive detection rate according to preset weights, and use it as the fitness value of the particle.

[0062] The particle swarm optimization algorithm iteratively updates the fitness value of each particle. The algorithm first records the highest fitness value ever achieved by each particle and its corresponding position vector, representing its individual optimal position. Simultaneously, it records the highest fitness value achieved by all particles and its corresponding position vector, representing its globally optimal position. In each iteration, for each particle, its velocity vector is updated according to the velocity update equation: the new velocity vector equals the current velocity vector multiplied by the inertia weight, plus the individual cognition term, the cognition coefficient multiplied by the difference between the individual optimal position and the current position, plus the social sharing term, the social coefficient multiplied by the difference between the globally optimal position and the current position. The position vector is updated as the sum of the new velocity vector and the current position vector. If a component of the updated position vector exceeds its value range, that component is truncated to the nearest boundary value, and the corresponding velocity component is set to zero. This iteration continues until the rate of change of the fitness value falls below a preset criterion in multiple consecutive iterations, or the number of iterations reaches a preset upper limit.

[0063] Furthermore, in S3, the associated site screening strategy is composed of a hierarchical filtering architecture and a fine-grained localization architecture connected in series; The hierarchical filtering architecture takes the whole genome variation sites in the genotype variation and corresponding phenotypic response data as the processing objects. According to the chromosome partitioning strategy, the whole genome is divided into several physical intervals of equal length. Then, the variation sites in each physical interval are sorted according to the phenotypic association strength and the candidate sites with the highest ranking are retained. The phenotypic association strength is determined by the ratio of the local phenotypic effect size to the global phenotypic effect size at the location of the variation site. The fine-targeting architecture receives the candidate sites output by the hierarchical filtering architecture, performs linkage disequilibrium decay analysis on all variant sites within a preset flanking length range around each candidate site, and identifies tag sites with independent genetic effects based on the law that the correlation coefficient between paired sites decays with physical distance. During the iteration process, the particle swarm optimization algorithm synchronously adjusts the granularity of the physical interval division and the sort retention ratio in the hierarchical filtering architecture, as well as the preset flanking length and the significance cutoff value of the chain imbalance decay analysis in the fine positioning architecture. The significance cutoff value is dynamically updated by the real-time changes in the fitness value.

[0064] In this embodiment, the associated site screening strategy consists of a hierarchical filtering architecture and a fine-tuning architecture connected in series, with both adjusting their internal parameters synchronously during the iteration of the particle swarm optimization algorithm.

[0065] The hierarchical filtering architecture takes the whole genome variation sites output by the simulation platform as the processing object; firstly, it executes the chromosome partitioning strategy, dividing the crop's reference genome into several equal-length intervals according to the preset length of each physical interval. For the terminal intervals that are shorter than the preset length, they are merged with the previous interval; after the partitioning is completed, the whole genome is divided into several physical intervals.

[0066] For all variant sites within each physical interval, the phenotypic association strength for each site is calculated. The phenotypic association strength is calculated as follows: for a variant site within the interval, the phenotypic response values ​​of all simulated plants carrying the variant at that site are extracted, and their average value is calculated; the phenotypic response values ​​of all simulated plants not carrying the variant at that site are extracted, and their average value is calculated; the difference between the two average values ​​is taken as the local phenotypic effect size for that site; simultaneously, the median value of the phenotypic effect size for all variant sites across the entire genome is calculated as the global phenotypic effect size; the phenotypic association strength for that site is equal to the ratio of the local phenotypic effect size to the global phenotypic effect size.

[0067] Within each physical interval, all variant sites are sorted from high to low according to the phenotypic association strength mentioned above; the sites with the highest ranking are retained, and the retention ratio is controlled by the ranking retention ratio parameter in the particle position vector; after the hierarchical filtering architecture performs the above operation on all physical intervals, all retained candidate sites are summarized and output to the fine-grained localization architecture; the interval division granularity parameter is adjusted synchronously during the particle swarm optimization process.

[0068] The fine-tuning architecture receives candidate sites output by the hierarchical filtering architecture; for each candidate site, it extracts all variant sites within a preset flanking length range around the site; the preset flanking length is also controlled by the corresponding parameter in the particle position vector; it performs linkage disequilibrium decay analysis on all extracted sites; the linkage disequilibrium decay analysis is calculated as follows: for any two sites within this range, the allele combination frequencies of these two sites in all simulated samples are statistically analyzed, and the correlation coefficient between them is calculated; the correlation coefficients of all site pairs are grouped according to the physical distance between the two sites, and the average correlation coefficient within each distance group is calculated; the decay law of the average correlation coefficient with physical distance is observed.

[0069] The fine-mapping architecture identifies tag sites with independent genetic effects based on the results of linkage disequilibrium decay analysis. The identification rule is as follows: within a preset flanking range, if the correlation coefficients between a group of sites are all higher than a preset threshold, these sites are classified into the same linkage block. From each linkage block, the site with the highest phenotypic association strength is selected as the tag site of that linkage block; the remaining non-tag sites are discarded. After processing by the fine-mapping architecture, the final output sites are the candidate variant site output set under this particle configuration. The significance cutoff value parameter is used in linkage disequilibrium decay analysis to determine whether the correlation coefficient is significantly non-zero. This parameter is dynamically updated by the particle swarm optimization algorithm during the iteration process based on the real-time changes in the fitness value.

[0070] Furthermore, in S4, the evaluation process of the detection effect is implemented in a stratified manner according to the mutagenic dose interval. The genotype variation and corresponding phenotypic response data generated by the simulation platform are divided into low-dose subset, medium-dose subset and high-dose subset according to the mutagenic dose value. Then, the detection recall and detection precision of the associated site screening strategy are calculated in each subset. The detection recall and detection precision are weighted and summarized according to the sample size of each subset to obtain the global detection validity index. The ratio between the global detection validity index and the number of known associated sites pre-labeled in the simulated truth set is used to determine whether the current strategy has reached the preset validity threshold. If the validity threshold is not reached, the process returns to S3 to re-execute the iterative optimization of the particle swarm optimization algorithm.

[0071] In this embodiment, the evaluation process of the detection effect is implemented in layers according to the mutagenic dose interval. First, all genotype variations and corresponding phenotypic response data generated by the simulation platform are divided into three subsets according to the mutagenic dose values: low-dose subset, medium-dose subset, and high-dose subset. The low-dose subset is defined as data with a dose less than a first preset dose value; the medium-dose subset is defined as data with a dose greater than or equal to the first preset dose value and less than or equal to the second preset dose value; and the high-dose subset is defined as data with a dose greater than the second preset dose value.

[0072] Calculate the recall and precision of the associated site screening strategy within each subset. Taking the low-dose subset as an example: count the number of real associated sites contained in the simulated truth set of this subset; run the associated site screening strategy under the globally optimal particle configuration on this subset and output the candidate site set; count the number of sites in the candidate site set that appear in the simulated truth set, and calculate the recall of the low-dose subset as this number divided by the total number of real associated sites in the subset; count the total number of digits in the candidate site set, and calculate the precision of the low-dose subset as the number of sites that appear in the truth set divided by the total number of digits in the candidate site set.

[0073] Perform the same calculations on the medium-dose subset: count the total number of true associated sites and the number of candidate site matches in the medium-dose subset, and calculate the detection recall and detection precision of the medium-dose subset; perform the same calculations on the high-dose subset: count the total number of true associated sites and the number of candidate site matches in the high-dose subset, and calculate the detection recall and detection precision of the high-dose subset.

[0074] The global detection validity index is obtained by weighting and summing the detection recall and detection precision of the three subsets according to the sample size of each subset. Specifically, the global detection recall is equal to the sum of the products of the detection recall of each subset and its sample size, divided by the total sample size; the global detection precision is equal to the sum of the products of the detection precision of each subset and its sample size, divided by the total sample size.

[0075] The ratio between the global detection validity index and the number of known associated loci pre-labeled in the simulated ground truth set is used as the criterion; this ratio is the global detection recall. This ratio is compared to a preset validity threshold. If the ratio reaches or exceeds the threshold, the current strategy is considered to have reached the validity threshold, and there is no need to return to S3 for re-optimization. If the threshold is not reached, the current iteration count is updated, and the algorithm is returned to S3 to re-execute the particle swarm optimization algorithm until the global detection recall reaches or exceeds the threshold.

[0076] Furthermore, in S5, the screening process of the optimal combination of mutagenesis conditions and mutation sites is performed in sequence as a preliminary screening stage of mutagenesis conditions and a fine screening stage of mutation sites. In the preliminary screening stage of mutagenesis conditions, the detection effect corresponding to each candidate mutagenesis condition obtained in S4 is used as the sorting basis, and multiple mutagenesis conditions with the highest detection effect are selected to form a candidate condition pool. In the mutation site screening stage, within the mutation site set corresponding to the candidate condition pool, the stability of each site in the multi-round Monte Carlo simulation on the simulation platform and the consistency of the response between each site and the real-time environmental parameters are used as two-dimensional ranking criteria. The mutation site with the highest comprehensive score under the two-dimensional ranking criteria is selected as the target site in the optimal combination. At the same time, the mutagenic dose and treatment duration corresponding to the target site are recorded as the target conditions in the optimal combination.

[0077] In this embodiment, the screening process sequentially performs the initial screening stage of mutagenic conditions and the fine screening stage of variant sites.

[0078] In the initial screening stage of mutagenic conditions, the detection effect of each group of candidate treatment conditions obtained in S4 is used as the ranking criterion. The detection effect of each group of candidate treatment conditions is represented by the subset detection recall rate corresponding to that group of conditions in the stratified evaluation. That is, the detection recall rate of the low-dose subset is used for each group of conditions in the low-dose subset, the detection recall rate of the medium-dose subset is used for each group of conditions in the medium-dose subset, and the detection recall rate of the high-dose subset is used for each group of conditions in the high-dose subset. All candidate treatment conditions are ranked from high to low according to the above detection effect scores. Multiple mutagenic conditions with the top ranking are selected to form a candidate condition pool, and the number of selected conditions is a preset value.

[0079] The mutation site fine screening stage involves screening the mutation sites within the candidate condition pool based on two ranking criteria.

[0080] The first dimension is the stability of each point in the multiple rounds of Monte Carlo simulation on the simulation platform; for a candidate point, the total number of times it appears in all simulation data is counted, the frequency of occurrence is calculated, and the frequency of occurrence is used as the stability score of the point.

[0081] The second dimension is the response consistency between each locus and real-time environmental parameters. For a candidate locus, the phenotypic response values ​​of all simulated plants carrying the variation of that locus are extracted under different environmental conditions, and the correlation between these response values ​​and the key environmental parameters in the corresponding environmental conditions is calculated. If the correlation is high, it indicates that the phenotypic effect of that locus is significantly regulated by environmental conditions and the response consistency is good. The absolute value of the correlation is used as the response consistency score of that locus.

[0082] A comprehensive score is calculated for all sites in the candidate condition pool according to the two dimensions mentioned above. The comprehensive score is calculated as a weighted sum of the stability score and the response consistency score, with the weight set to a preset ratio. The mutation site with the highest comprehensive score is selected as the target site in the optimal combination. At the same time, the mutagenic dose and treatment duration corresponding to the target site are recorded as the target conditions in the optimal combination. Thus, an optimal combination of mutagenic conditions and mutation sites is obtained.

[0083] Furthermore, after S5, a post-verification step is performed. When the real-time environmental parameters and germplasm phenotypic data continuously collected by the sensors and detectors deployed in the crop field for the next time window arrive, the post-verification step inputs the target site in the optimal combination of mutagenesis conditions and mutation sites into the genotype-phenotype association physical model, driving the physical model to forward deduce the expected phenotypic change trajectory within the next time window. The expected phenotypic change trajectory is compared with the measured germplasm phenotypic data in the next time window. If the deviation between the two exceeds the preset allowable range, the simulation platform is triggered to re-execute the Monte Carlo simulation method with the real-time environmental parameters of the next time window as new initial conditions. The particle swarm optimization algorithm is driven by the updated data obtained after re-execution to modify the optimal combination of mutation conditions and mutation sites.

[0084] In this embodiment, the post-verification step is triggered when the real-time environmental parameters and germplasm phenotypic data collected by the sensors and detectors deployed in the crop field for the next time window arrive; the length of the time window is set to a preset duration.

[0085] When the first batch of real-time data for the next time window arrives, the target locus in the obtained optimal combination is input into the genetic expression module of the genotype-phenotype association physical model. Specifically, in the genetic parameter library, the effect value of the gene corresponding to the target locus is forcibly set to the effect value recorded in the optimal combination. Keeping other parameters of the physical model unchanged, the real-time environmental parameters at the current moment are used as the initial boundary conditions to drive the physical model to forward extrapolate the expected phenotypic change trajectory for each day in the next time window. The forward extrapolation method is as follows: starting from the current moment, the environmental parameters of each time step are input into the physical model sequentially at a preset time step, the expected phenotypic value is updated step by step, and the extrapolation continues until the end of the time window, resulting in an expected phenotypic change trajectory containing multiple time points.

[0086] Meanwhile, in the next time window, the sensors and detectors deployed in the field continue to collect measured germplasm phenotypic data at the usual collection frequency; at the end of the next time window, all measured phenotypic data collected in the time window are summarized to form a trajectory of measured phenotypic changes that includes the same time points.

[0087] The expected phenotypic change trajectory is compared with the measured phenotypic change trajectory at each time point, and the deviation between the two is calculated. The deviation is calculated as follows: for each time point, the absolute difference between the expected phenotypic value and the measured phenotypic value is calculated; the average of the absolute differences at all time points is taken as the average deviation.

[0088] The average deviation is compared with a preset tolerance range, which is set as a preset proportion of the phenotypic dispersion index. If the average deviation exceeds the tolerance range, the deviation is determined to exceed the preset tolerance range.

[0089] The simulation platform is triggered to re-execute the Monte Carlo simulation method using the real-time environmental parameters of the next time window as new initial conditions. The re-execution method is as follows: the environmental parameter sequence actually observed in the next time window is used as the new environmental input, replacing the original environmental parameter probability distribution, and the Monte Carlo simulation process is re-run, including re-extracting candidate treatment condition combinations, re-executing the genome damage simulation, and regenerating phenotypic response data. The updated data obtained after the re-execution drives the particle swarm optimization algorithm to correct the optimal combination, that is, the particle swarm optimization process is re-run using the updated simulation data to obtain the corrected optimal mutagenesis conditions and optimal mutation sites. If the average deviation does not exceed the allowable range, there is no need to trigger a re-simulation, and the current optimal combination is confirmed as valid and stored.

[0090] The embodiments of this application have been described above with reference to the accompanying drawings. However, this application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive.

[0091] Those skilled in the art, inspired by this application, may make many other forms without departing from the spirit and scope of the claims, all of which are within the protection of this application.

Claims

1. A method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources, characterized in that, include: S1: Acquire real-time environmental parameters and germplasm phenotypic data from sensors and detectors deployed in the crop field, and dynamically map them to a preset genotype-phenotype association physical model, using a set of differential equations to describe the dynamic coupling relationship between environmental factors and germplasm phenotype. S2, Analyze the real-time environmental parameters, the germplasm phenotypic data and the acquired historical germplasm observation data, construct a simulation platform with input interface for mutagenesis condition parameters and output interface for mutation sites, and use Monte Carlo simulation method to generate genotypic variation and corresponding phenotypic response data under different combinations of mutagenesis doses and treatment conditions; S3. A particle swarm optimization algorithm is used to formulate an associated site screening strategy. The genotype variation and corresponding phenotypic response data output by the simulation platform are used as input to adjust the parameter configuration in the identification process in order to maximize the validity of the variant site identification. S4. After executing the associated site screening strategy, evaluate the detection effect of the strategy on candidate variant sites. The detection effect is measured by comparing the degree of agreement between the preset simulated truth set in the simulation platform and the actual output set of the strategy. S5. Based on the evaluation results of the detection effect, screen and determine the optimal combination of mutagenesis conditions and variant sites to achieve the highest genotype-phenotype association resolution accuracy of crop germplasm under the current environment and genetic background.

2. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In S1, the real-time environmental parameters are collected wirelessly and synchronously by multiple sets of sensor nodes deployed in different ecological locations in the field. The sensor nodes output environmental temperature, environmental humidity, light intensity, wind speed and direction and soil electrical conductivity parameters at preset time intervals. The detector outputs crop canopy spectral reflectance and leaf surface temperature distribution parameters at preset spatial intervals. The real-time environmental parameters and the germplasm phenotypic data are time-stamp aligned and spatial coordinate registered before being entered into the physical model.

3. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In step S1, the genotype-phenotype association physical model includes an environmental response module and a genetic expression module. The environmental response module receives the real-time environmental parameters and outputs an environmental stress vector. The genetic expression module receives the environmental stress vector and outputs a phenotypic plasticity interval. The dynamic mapping substitutes the real-time environmental parameters at each sampling time into the environmental response module to solve for the transient environmental stress value. Then, the transient environmental stress value is used as a boundary condition to input into the genetic expression module to update the expected phenotypic value at the current time. The residual between the expected phenotypic value and the germplasm phenotypic data measured by the detector is recorded in real time and used to drive the parameter fine-tuning of the physical model.

4. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In S2, the simulation platform constructs a digital twin field environment based on the real-time environmental parameters, the germplasm phenotypic data, and the historical germplasm observation data. The digital twin field environment includes a soil profile stratification model, a microclimate boundary layer model, and a crop growth kinetics model. The soil profile stratification model stores water transport parameters and nutrient diffusion parameters at different soil depths. The microclimate boundary layer model updates the near-surface radiation balance and turbulent exchange coefficients in a diurnal variation cycle and weather event-driven manner. The crop growth kinetics model outputs the crop organ formation rate and biomass allocation ratio with accumulated temperature accumulation and light energy utilization efficiency as core state variables. The Monte Carlo simulation method is implemented as follows: The probability distribution functions of the mutagenic dose and the treatment duration were determined based on historical mutagenic experimental records. Random samples are independently drawn from the probability distribution function of the mutagenic dose and the probability distribution function of the treatment duration to form a set of candidate treatment condition combinations; For each candidate treatment combination, a pre-defined genomic damage simulation subroutine is executed in the digital twin field environment. The genomic damage simulation subroutine generates a set of single nucleotide variant sites and a set of structural variant sites based on the randomness of radiation energy deposition and the competitive pathway of DNA double-strand break repair. The set of single nucleotide variant sites and the set of structural variant sites are input into the phenotypic conversion subroutine in the crop growth kinetic model. The phenotypic conversion subroutine calculates the quantitative perturbation value of the variant on the organ formation rate and the biomass allocation ratio based on the regulatory network position and metabolic pathway affiliation of the gene where the variant site is located. The set of single nucleotide variant sites, the set of structural variant sites, and the corresponding quantitative perturbation values ​​under each candidate treatment condition combination are recorded as the genotype variation and corresponding phenotypic response data.

5. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In S2, the simulation platform further includes a field heterogeneous noise injection module, which generates a spatially continuous noise field based on the spatial variation semivariance function of the real-time environmental parameters, and superimposes the noise field onto the corresponding phenotypic response data output by the Monte Carlo simulation method. The spatial variation semivariance function is obtained by fitting repeated measurement records from the historical germplasm observation data, and the field heterogeneity noise injection module regenerates an independent noise field in each Monte Carlo simulation iteration.

6. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In step S3, the particle swarm optimization algorithm encodes the parameters to be adjusted in the identification process into a particle position vector, which includes a mutation site significance threshold, a phenotypic response amplitude cutoff coefficient, and a site co-occurrence frequency threshold. After the particle swarm optimization algorithm initializes the particle population, each particle runs the associated site screening strategy based on the genotype variation and corresponding phenotypic response data output by the simulation platform, and calculates the positive detection rate and false positive detection rate of the strategy on the simulation platform. The weighted harmonic mean of the positive detection rate and the false positive detection rate is used as the fitness value. The particle swarm optimization algorithm iteratively updates the optimal position of individual particles and the global optimal position of the swarm based on the fitness value, and adjusts the flight direction and step size of each particle according to the velocity update equation until the rate of change of the fitness value is lower than the preset convergence criterion or the number of iterations reaches the preset upper limit.

7. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In S3, the associated site screening strategy is composed of a hierarchical filtering architecture and a fine-grained localization architecture connected in series. The hierarchical filtering architecture takes the whole genome variation sites in the genotype variation and corresponding phenotypic response data as the processing objects. According to the chromosome partitioning strategy, the whole genome is divided into several physical intervals of equal length. Then, the variation sites in each physical interval are sorted according to the phenotypic association strength and the candidate sites with the highest ranking are retained. The phenotypic association strength is determined by the ratio of the local phenotypic effect size to the global phenotypic effect size at the location of the variation site. The fine-targeting architecture receives the candidate sites output by the hierarchical filtering architecture, performs linkage disequilibrium decay analysis on all variant sites within a preset flanking length range around each candidate site, and identifies tag sites with independent genetic effects based on the law that the correlation coefficient between paired sites decays with physical distance. During the iteration process, the particle swarm optimization algorithm synchronously adjusts the granularity of the physical interval division and the sort retention ratio in the hierarchical filtering architecture, as well as the preset flanking length and the significance cutoff value of the chain imbalance decay analysis in the fine positioning architecture. The significance cutoff value is dynamically updated by the real-time changes in the fitness value.

8. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In step S4, the evaluation process of the detection effect is implemented in a stratified manner according to the mutagenic dose interval. The genotype variation and corresponding phenotypic response data generated by the simulation platform are divided into low-dose subset, medium-dose subset and high-dose subset according to the mutagenic dose value. Then, the detection recall and detection precision of the associated site screening strategy are calculated in each subset. The detection recall and detection precision are weighted and summarized according to the sample size of each subset to obtain the global detection validity index. The ratio between the global detection validity index and the number of known associated sites pre-labeled in the simulated truth set is used to determine whether the current strategy has reached the preset validity threshold. If the validity threshold is not reached, the process returns to S3 to re-execute the iterative optimization of the particle swarm optimization algorithm.

9. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 1, characterized in that, In S5, the screening process of the optimal combination of mutagenesis conditions and mutation sites is carried out in sequence, including the initial screening stage of mutagenesis conditions and the fine screening stage of mutation sites. In the initial screening stage of mutagenesis conditions, the detection effect corresponding to each candidate mutagenesis condition obtained in S4 is used as the sorting basis, and multiple mutagenesis conditions with the highest detection effect are selected to form a candidate condition pool. In the mutation site screening stage, within the mutation site set corresponding to the candidate condition pool, the stability of each site in the multi-round Monte Carlo simulation on the simulation platform and the consistency of the response between each site and the real-time environmental parameters are used as two-dimensional ranking criteria. The mutation site with the highest comprehensive score under the two-dimensional ranking criteria is selected as the target site in the optimal combination. At the same time, the mutagenic dose and treatment duration corresponding to the target site are recorded as the target conditions in the optimal combination.

10. The method for identifying genotype-phenotype association variation sites in space-induced mutagenesis germplasm resources according to claim 9, characterized in that, After S5, a post-verification step is performed. When the real-time environmental parameters and germplasm phenotypic data of the next time window continuously collected by the sensors and detectors deployed in the crop field arrive, the post-verification step inputs the target site in the optimal combination of mutagenesis conditions and mutation sites into the genotype-phenotype association physical model, driving the physical model to positively deduce the expected phenotypic change trajectory within the next time window. The expected phenotypic change trajectory is compared with the measured germplasm phenotypic data in the next time window. If the deviation between the two exceeds the preset allowable range, the simulation platform is triggered to re-execute the Monte Carlo simulation method with the real-time environmental parameters of the next time window as new initial conditions. The particle swarm optimization algorithm is driven by the updated data obtained after re-execution to modify the optimal combination of mutation conditions and mutation sites.

Images