Molecular navigation and orientation design method for agricultural biological gene regulation and control region and application of molecular navigation and orientation design method

Abstract

The invention discloses a molecular navigation and orientation design method for an agricultural biological gene regulation region and application. The method comprises the following steps: firstly, aiming at a target gene, screening out an open chromatin region in a regulation window of the target gene; secondly, on the basis of local genetic variation, the accessibility of each candidate OCR is decomposed into a cis component (CAcis) and a trans component (CAtrans); thirdly, performing two-channel correlation analysis on the expression quantity of the target gene and CAcis and CAtrans respectively; and finally, when and only when statistics of the same OCR in the two channels is significant and regression coefficient symbols are consistent, determining that the OCR is a causal cis-form regulation region, and marking the OCR as a positive or negative regulation element according to the regulation region. The method is strong in causality and clear in direction, can directly output a directional design scheme linked with genome editing, and provides an efficient and reliable technical approach for precise genetic improvement of agricultural biological complex characters.

Description

Technical Field

[0001] This invention relates to the fields of plant molecular biology and bioinformatics, specifically to a method and application of molecular navigation and directional design of gene regulatory regions in agricultural organisms. Background Technology

[0002] The spatiotemporal specificity of gene expression is the molecular basis for determining plant growth, development, morphogenesis, and environmental adaptation. This precise regulation of expression is mainly accomplished by two types of elements working together: one is trans-acting factors encoding proteins (such as transcription factors), and the other is cis-regulatory elements (CREs) located in the DNA sequence. CREs include promoters, enhancers, and repressors, which recruit transcription machinery by binding to trans-acting factors, thereby activating or inhibiting the transcription of downstream target genes. Numerous studies have shown that genetic variation in cis-regulatory elements located in non-coding regions is an important source of phenotypic differences between individuals and a key selection target in plant domestication and modern breeding improvement.

[0003] For example, in rice, FZP , OsLG1 , GW5 and qSH1 Variations in the cis-regulatory regions of non-coding genes have been shown to influence key agronomic traits such as panicle type, grain type, and shattering tendency by regulating gene expression levels. Unlike coding region variations, which lead to protein function loss or harmful pleiotropic effects, modifications to non-coding cis-regulatory regions can achieve "fine-tuning" of gene expression levels. This allows for the optimization of agronomic traits without disrupting basic gene function, providing a broader and more flexible scope for genetic improvement. (High-yield rice genes) Gn1a Ideal plant type gene IPA1 Examples such as [examples of examples] demonstrate how variations in their cis-regulatory regions can moderately alter expression levels to increase breeding yield. In recent years, the development of genome editing technologies, exemplified by the CRISPR / Cas system, has made precise targeted modification of these non-coding regulatory regions possible, showcasing enormous application potential.

[0004] However, current genetic modification targeting non-coding regions still faces significant challenges. First, the precise location, regulatory scope, and functional strength of the cis-regulatory landscape of most genes, particularly distal regulatory elements (such as enhancers) far from the transcription start site (TSS), remain unclear. Existing research largely focuses on promoter regions near genes, often employing multi-site "trial" editing designs, lacking systematic guiding principles, resulting in low efficiency and unpredictable improvement effects. Second, whether a regulatory element promotes or inhibits gene expression—that is, its regulatory "direction"—is often unknown. This makes the selection of editing targets somewhat arbitrary, failing to achieve predictable trait improvements.

[0005] At the methodological level, while three-dimensional genomics techniques (such as Hi-C) help reveal spatial interactions in chromatin and delineate potential "promoter-enhancer" pairs, their resolution is limited, and spatial proximity is not entirely equivalent to direct functional regulation, making it difficult to precisely locate specific element sequences, let alone determine their regulatory direction. On the other hand, association analysis methods based on population data, such as expression quantitative trait loci (eQTL) mapping, can link gene expression differences to genetic variation sites on the genome. However, due to the widespread linkage disequilibrium (LD) effect, they can usually only locate a large genomic region, making it difficult to identify the true causal variations within the region and the regulatory elements they contain.

[0006] In summary, there is an urgent need in this field for a methodology that can systematically and accurately locate the causal cis-regulatory regions of genes across the entire genome, clarify their positive / negative regulatory directions, and directly guide genome editing practices, in order to break through the current bottleneck in non-coding region genetic improvement and accelerate the process of precision design breeding of complex crop traits. Summary of the Invention

[0007] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a method and application for molecular navigation and targeted design of gene regulatory regions in agricultural organisms. This invention integrates gene expression profiles, chromatin accessibility maps, and genomic variation data from multiple populations. Through an innovative causal inference framework, it achieves high-precision localization of the causal cis-regulatory regions of target genes and determination of the regulatory direction (positive / negative). This information is then transformed into a systematic method, software system, and application for operable genome editing design schemes. This invention can be used not only for the genetic improvement of major crops such as rice, maize, and wheat, but also for basic biological research in other plants, such as gene regulatory network analysis and exploration of domestication mechanisms.

[0008] To achieve the above objectives, the technical solution designed by the present invention is as follows: This invention provides a method for molecular navigation and targeted design of gene regulatory regions in agricultural organisms, comprising the following steps: S1. Obtain whole-genome gene expression profiles, chromatin accessibility maps, and whole-genome genotype data of an associated population of a certain agricultural organism; an associated population refers to a biological population with extensive genetic variation and phenotypic diversity, which can be a natural population or an artificially created population for genetic research. S2. Based on whole-genome gene expression profile data, the expression level of each gene in each sample of the population is quantified, and a gene expression matrix is ​​constructed. S3. Based on chromatin accessibility map data, identify open chromatin regions across the entire genome, quantify the accessibility signal intensity of each open chromatin region in each sample of the population, and construct an OCR accessibility matrix. S4. Using the transcription start sites of all genes in the gene expression matrix as anchors, define the upstream and downstream linear windows of each gene as its candidate cis-regulatory range. For each target gene, extract all OCRs whose coordinates lie within its candidate cis-regulatory range from the OCR accessibility matrix, and use them as candidate regulatory regions for that gene. S5. Perform whole-genome filtering on the whole-genome genotype data and construct a whole-genome kinship matrix K_G; S6. For each candidate regulatory region of a gene, construct a local kinship matrix K_local using the genetic variations in the genomic regions adjacent to the candidate regulatory region in the whole genome genotype data; S7. Chromatin accessibility analysis was performed on the candidate regulatory regions of each gene based on the whole genome phylogenetic matrix and the local phylogenetic matrix; and cis-heritability testing was performed on the candidate regulatory regions of each gene to screen out the OCRs of each gene that "have significant cis-heritability". S8. For each OCR that “has significant cis-heritability” selected in step S7, apply the statistical model used in the cis-heritability test to decompose the chromatin accessibility of the OCR into the cis component CA_cis determined by local genetic variation and the trans component CA_trans determined by other factors. S9. For each gene and its corresponding OCR that "has significant cis-heritability", form a "gene-OCR" pair, and perform a dual-channel association analysis of cis and trans channels for each "gene-OCR" pair: Cis-channel: The model is in the form of y = μ + β_cis CA_cis + u_G + ε; Inverse channel: The model is in the form of y = μ + β_trans CA_trans + u_G + ε; Where y is the expression level of the gene component in the gene-OCR pair; CA_cis and CA_trans are the cis and trans components of the OCR component in the gene-OCR pair, respectively; β_cis and β_trans are the fixed-effect regression coefficients of CA_cis and CA_trans on the gene expression level y, respectively; μ is the intercept term; u_G is the random effect used to correct for population structure; and ε is the model residual. S10. Based on the above dual-channel correlation analysis, the following judgment is made: When the P-value of the OCR portion in the same "gene-OCR" pair is less than the statistical significance threshold in both of the above two models, and its regression coefficients β_cis and β_trans have the same sign, the OCR is determined to be the causal cis-regulatory region of the gene. If both regression coefficients β_cis and β_trans are positive, it is a positive control region; if both regression coefficients β_cis and β_trans are negative, it is a negative control region.

[0009] Furthermore, it also includes the following steps: (1) Calculate the z-score statistics z_cis and z_trans corresponding to the two models of the cis-channel and trans-channel, and calculate the comprehensive score. Sort the multiple causal cis-regulatory regions of a gene according to the comprehensive score. (2) Output a list of targets including the coordinates of the regulatory region, the direction of regulation, the statistical significance and the priority score, and generate a specific genome editing design scheme based on it to achieve the expected regulation of the expression level of the target gene; The design scheme includes: (a) To achieve upregulation of gene expression: design strategies designed to disrupt or weaken the function of an identified negative regulatory region, or design strategies designed to enhance the function of an identified positive regulatory region.

[0010] (b) To achieve downregulation of gene expression: design strategies designed to disrupt or weaken the function of an identified positive regulatory region, or design strategies designed to enhance the function of an identified negative regulatory region.

[0011] Furthermore, in step S1, the agricultural organism is any one of rice, corn, wheat, barley, sorghum, soybean, rapeseed, cotton, tomato, and potato; The number of members in the associated group is not less than 200; Whole-genome gene expression profiles were obtained through RNA-seq sequencing; chromatin accessibility profiles were obtained through ATAC-seq sequencing; and whole-genome genotype data were obtained through whole-genome resequencing.

[0012] Furthermore, the agricultural organism is rice; in step S2, the specific steps for quantifying the expression level of each gene in each sample and constructing the gene expression matrix are as follows: S21. Using rapid quantitative tools, RNA-seq sequencing reads are quantified based on the annotated transcript sequences of the target species' reference genome to count the number of reads for each gene. S22. Calculate the TPM of each gene and combine the TPM values ​​of different transcripts of the same gene as the total TPM of the gene. Perform log2(TPM+1) transformation on the total TPM value of the gene and screen genes expressed in at least 5% of the samples for downstream analysis to obtain the original gene expression matrix. S23. The original gene expression matrix is ​​corrected using a latent factor model, and the corrected expression values ​​of each gene in all samples are subjected to a rank-based inverse normal transformation to obtain the gene expression matrix.

[0013] Furthermore, in step S3, the specific steps for identifying open chromatin regions across the entire genome and quantifying the accessibility signal intensity of each open chromatin region in each sample to construct an OCR accessibility matrix are as follows: S31. After aligning the ATAC-seq sequencing reads to the reference genome of the target species, the alignment results are filtered and refined, and high-quality aligned sequencing fragments are selected. Then, Tn5 transposon cleavage site offset correction is performed. S32. Using a 250 bp sliding window with a step size of 100 bp, the adjusted transposon cut sites are counted across the entire genome to quantify chromatin accessibility and obtain the original "window × sample" counting matrix. S33. Perform quantile normalization on the original "window × sample" count matrix; S34. Based on the normalized signal value, calculate an index representing the overall accessibility level of each window, sort all windows based on this index, and select the 20% of windows with the strongest signal as open chromatin regions. S35. For the selected open chromatin regions, their normalized count values ​​are logarithmically transformed after adding a pseudo-count of 1. S36. The latent factor model is used to remove the influence of latent confounding factors, and the accessibility value of each corrected open chromatin region is subjected to a rank-based inverse normal transformation to obtain the OCR accessibility matrix.

[0014] Furthermore, in step S4, the physical distance of the linear window is 5 kb to 1000 kb; In step S5, the condition for whole-genome filtering is that the minor allele frequency (MAF) > 0.05. In step S6, the range of adjacent genomic regions is 5 kb to 1000 kb.

[0015] Furthermore, in step S7, the specific steps for chromatin accessibility analysis are as follows: S71. For each candidate regulatory region of a gene, fit the following variance component model to analyze its accessibility: x = μ + u_cis + u_G + ε; Where x is chromatin accessibility, μ is the intercept term, u_cis follows a distribution of N(0, σ²_cis K_local), representing a cis-random effect determined by local genetic variation, u_G follows a distribution of N(0, σ²_G K_G), representing a random effect determined by the whole genome genetic background and used to control population structure and background polygenic effects, and ε follows a distribution of N(0, σ²_e I), representing the residual; The specific steps for screening each gene for "significant cis-heritability" using OCR are as follows: S72. Obtain the test P value of the cis variance component σ²_cis for each candidate regulatory region of each gene, perform multiple test correction on all P values, and define the candidate regulatory region with FDR < 0.05 as "having significant cis heritability" OCR. In step S8, the specific steps for obtaining the cis component CA_cis and the trans component CA_trans are as follows: S81. For the OCR that "has significant cis-heritability" in step S72, extract the best linear unbiased predictor value of its cis-random effect from the variance component model in step S71. This value is defined as the cis-component CA_cis of the OCR accessibility. Calculate the residual x. CA_cis, the residual is defined as the trans component CA_trans; The other factors refer to all factors other than local genetic variations.

[0016] Furthermore, in step S10, the statistical significance threshold is 1×10⁻⁶. -4 .

[0017] Furthermore, in step (1), the comprehensive score S is calculated according to the following formula: S = (z_cis² z_trans²) / (z_cis² + z_trans²).

[0018] The present invention also provides an application of the method in plant genetic improvement.

[0019] The beneficial effects of this invention are: The method of this invention is based on a core biological assumption: a true causal cis-regulatory element, whose activity changes, whether stemming from cis-genetic variation in its own sequence or from perturbations of nonlocal factors such as the abundance of trans-acting factors in the cellular environment, should produce a directionally consistent regulatory effect on the transcription of the target gene. This invention, through its unique framework of "cis / trans decomposition + dual-channel consistency test," achieves the following significant beneficial effects: 1. Strong causal inference capability and high accuracy: This invention cleverly utilizes perturbations from two different sources, cis and trans, as a "quasi-experiment." A true causal regulatory element should exhibit consistent effects on gene expression under both perturbations. This stringent filtering condition effectively eliminates a large number of false positive signals caused by tightly linked but non-functional sites or upstream co-regulatory factors, thereby significantly improving the causality and accuracy of the identification results.

[0020] 2. Clear regulatory directionality and strong guidance: This invention can not only "find" the regulatory region, but also "understand" its function—whether it promotes or inhibits. This directional information is key to achieving predictable trait improvement. It transforms genome editing from a "trial and error" mode to a "directed design" mode, greatly improving research and development efficiency and success rate.

[0021] 3. Broad Applicability and Good Extensibility: The core logic of this invention does not depend on specific species or specific data generation technologies. This method can be applied as long as the genome, transcriptome, and chromatin accessibility data of a population are available. Its algorithm modules (such as the cis / trans decomposition model) are replaceable and can naturally integrate more complex genetic variation information (such as structural variations), ensuring the universality and forward-looking nature of the method.

[0022] 4. Close integration of industry, academia, research, and application, with high application value: This invention directly addresses the practical needs of molecular breeding, and its final output is a rich list of editing target designs that can be directly used to guide experiments. This invention establishes a technological link from massive omics data to precise breeding decisions, and can serve as a core engine to drive a new breeding model of "biological big data + genome editing," possessing significant application value and commercial prospects. Attached Figure Description

[0023] Figure 1 This is a comparative analysis of the method of this invention with conventional association methods, and a graph showing the results of whole-genome statistical analysis. In the figure, A shows the comparison of the number of genes with significant associations identified by different association strategies; B is a graph comparing the distribution of the number of significant OCRs associated with each gene; C is a distance distribution map of significant OCRs identified by different strategies relative to their target gene transcription start sites (TSS); D is a genomic location annotation diagram of the "most significant causal OCR" determined by the method of this invention; Figure 2 The method of this invention is for the ideal plant type gene in rice. IPA1 Application case analysis diagrams on the above; In the diagram, A is... IPA1 Causal association signal scan diagram of gene loci; B and C are locked. IPA1 dOCR dual-channel correlation verification result image; Figure 3 To Figure 2 Identified in IPA1 Experimental verification results of CRISPR / Cas9 editing in the distal causal regulatory region; In the diagram, A represents the locked state. IPA1 A schematic diagram of CRISPR / Cas9 editing in the remote causal regulation region (dOCR); B represents wild type and IPA1 CR-dOCR Comparison of plant morphology at maturity of homozygous mutants; C represents wild type and IPA1 CR-dOCR Comparison of main ear morphology of homozygous mutants at maturity; D represents wild type and IPA1 CR-dOCR Quantitative statistical chart of effective panicle number in homozygous mutants; E represents wild type and IPA1 CR-dOCR Quantitative statistical chart of the number of spikelets per panicle in homozygous mutants; The gray bars represent the mean, the error bars represent the standard error, and each scatter point represents the data of an independent plant; n is the sample size, and the p-value is calculated by a two-tailed Student's t-test. Detailed Implementation

[0024] The present invention will now be described in further detail with reference to specific embodiments, so that those skilled in the art can understand it.

[0025] Example 1: A method for molecular navigation and directional design of gene regulatory regions in agricultural organisms This embodiment selects the key agronomic trait gene IPA1 (gene ID: LOC_Os08g39890, ideal plant type gene) in rice as a representative case to perform high-precision localization of the causal regulatory region and determination of the regulatory direction. It is known that high expression of IPA1 can lead to a decrease in the number of tillers and an increase in the number of grains per panicle.

[0026] I. Data Acquisition and Processing This embodiment uses a natural population comprising 275 different rice varieties (https: / / ngdc.cncb.ac.cn / bioproject / browse / PRJCA012684). For each variety, tissue samples were collected when the young panicles were 1-2 mm in size.

[0027] 1. Data Acquisition (1) RNA-seq sequencing was performed on the above 275 samples to obtain whole-genome gene expression profile data; (2) ATAC-seq sequencing was performed on 219 of the samples to obtain chromatin accessibility map data; (3) More than 10 million high-quality genetic variations (mainly SNPs and InDels) were obtained from 275 samples through whole-genome resequencing as whole-genome genotype data.

[0028] 2. Gene expression treatment (1) Quantitative analysis: Using rapid quantitative tools (such as Salmon), based on the annotated transcript sequences of the Nipponbare reference genome (available from https: / / rice.uga.edu / ), the raw reads after RNA-seq sequencing of 275 samples were quantified to count the number of reads for each gene.

[0029] (2) Standardization and Filtering: The gene read counts obtained in the previous step were standardized. First, the TPM (Transcripts Per Million) value was calculated to correct for differences in library size between different samples, and the TPM values ​​of different transcripts of the same gene were merged as the total TPM value for that gene. Then, the total TPM value was transformed by log2(TPM+1). Finally, genes expressed in at least 5% of the varieties (TPM>0.1) were screened, resulting in 30,814 genes, thus forming a gene expression matrix that had undergone preliminary standardization.

[0030] (3) Correction and Transformation: To eliminate unknown confounding factors in the population and ensure the data meets the requirements of subsequent statistical analysis, the gene expression matrix formed in the previous step was further processed. First, the matrix was corrected using the latent factor model (PEER) to remove the systematic influence of the first three latent factors. Then, the corrected expression values ​​of each gene in all varieties were subjected to a rank-based inverse normal transformation to make its distribution conform to a normal distribution, finally obtaining a fully corrected gene expression matrix for downstream analysis.

[0031] 3. Chromatin accessibility data processing (1) Data preprocessing and alignment: Raw reads from ATAC-seq sequencing were aligned to the Nipponbare reference genome. Subsequently, the alignment results were filtered and refined, including merging reads from technically duplicated samples, removing PCR duplicates and reads derived from mitochondria and chloroplasts, and selecting high-quality alignments (MAPQ>30). To accurately determine the cleavage center of the Tn5 transposase, the transposon cleavage sites at both ends of the selected high-quality alignment fragments were specifically adjusted according to the strand in which they were located: 4 bp forward for the positive strand and 5 bp backward for the negative strand.

[0032] (2) Window counting and initial matrix generation: A 250 bp sliding window with a step size of 100 bp is used to count the transposon cut sites adjusted in the previous step across the whole genome to quantify chromatin accessibility, thereby generating an initial “window × sample” counting matrix.

[0033] (3) Standardization and OCR Screening: The initial count matrix generated in the previous step is processed. First, quantile normalization is used to correct for global technical differences between different samples. Then, based on the normalized signal values, an index representing the overall accessibility level of each window (i.e., the 95th quantile of the signal value of that window among all samples) is calculated. All windows are sorted according to this index, and the 20% of windows with the strongest signals are selected and defined as Open Chromatin Regions (OCRs) for subsequent analysis.

[0034] (4) Final Correction and Transformation: The selected OCRs and their corresponding normalized counts are processed. First, a pseudo-count of 1 is added to the counts, followed by a log2 transformation, i.e., log2(a+1), where a is the normalized count. Next, the data is corrected using a latent factor model (PEER) to remove the influence of unknown confounding factors from the first five latent factors. Finally, a rank-based inverse normal transformation is performed on the accessibility values ​​of each corrected OCR to make its distribution conform to a normal distribution, ultimately obtaining the OCR accessibility matrix for downstream analysis.

[0035] II. Candidate Regulatory Region Definition, Relationship Matrix Construction, and Cis / Trans Component Decomposition 1. Candidate Range and Relationship Matrix Construction (1) Using the transcription start site (TSS) of gene IPA1 in the gene expression matrix of Example 1 as the anchor point, a linear window with a distance of 100 kb upstream and downstream of gene IPA1 is defined as its candidate cis-regulatory range. For gene IPA1, all OCRs whose coordinates are located within its candidate cis-regulatory range are extracted from the OCR accessibility matrix and used as the candidate regulatory region (candidate OCR) of the gene. A linear window refers to a region of fixed length extending upstream (opposite to the direction of gene transcription) and downstream (in the direction of gene transcription) along the DNA "line" molecule, with the transcription start site of a gene as the reference point. This region is the "linear window".

[0036] (2) The whole genome genotype data (including SNPs and InDel) were filtered by the whole genome, i.e., the minor allele frequency MAF>0.05, and the whole genome kinship matrix K_G was constructed.

[0037] (3) For each candidate OCR, the corresponding local kinship matrix K_local is constructed by using the genetic variation data within its nearest ±100 kb range.

[0038] 2. Decomposition of cis / trans components of chromatin accessibility (LMM-BLUP example) (1) For each candidate OCR, fit the following variance component model to analyze its accessibility (denoted as x): x = μ + u_cis + u_G + ε; Where μ is the intercept term, u_cis follows a distribution of N(0, σ²_cis K_local), representing a cis-random effect determined by local genetic variation, u_G follows a distribution of N(0, σ²_G K_G), representing a random effect determined by the whole genome genetic background and used to control population structure and background polygenic effects, and ε follows a distribution of N(0, σ²_e I), representing the residual.

[0039] (2) The p-value of the cis-variance component σ²_cis of each candidate OCR was obtained through the likelihood ratio test (LRT). Subsequently, multiple test corrections were performed on all p-values. Discovery Rate (FDR) was used to define candidate OCRs with FDR < 0.05 as “having significant cis-heritability” (i.e., their σ²_cis is significantly greater than zero) and they were retained for downstream analysis.

[0040] (3) For the aforementioned OCRs that "have significant cis-heritability", the best linear unbiased predictor value (BLUP of u_cis) of the cis-random effect is extracted from the aforementioned variance component model. This value is defined as the cis-component CA_cis of the OCR accessibility. The cis-component represents the portion of accessibility variation that is entirely determined by the local genetic variation of the OCR. Subsequently, the residual x is calculated. CA_cis, the residual is defined as the trans component CA_trans, which comprehensively reflects the influence of all other factors (such as distant genetic loci, environmental signals, cell state, and fluctuations in the abundance of shared transcription factors) on the accessibility of OCR, except for local genetic variations.

[0041] III. Dual-channel correlation, causal determination, and priority ranking 1. For gene IPA1 and its corresponding candidate OCRs, “gene IPA1-candidate OCR” pairs are formed. The “gene IPA1-OCR” pairs in which the OCR portion passes the cis-heritability test in step 2 (i.e., FDR < 0.05) are selected. These selected OCRs are defined as “having significant cis-heritability” OCRs.

[0042] Subsequently, a two-channel association analysis was performed on each "gene IPA1-OCR" pair containing an "significantly cis-heritable" OCR. In this embodiment, this analysis was achieved by fitting the following two independent mixed linear models: Channel 1 (Civic Channel): y = μ + β_cis CA_cis + u_G + ε Channel 2 (Inverse Channel): y = μ + β_trans CA_trans + u_G + ε Where y represents the expression level of the gene portion in the "IPA1-OCR" pair; CA_cis and CA_trans represent the cis and trans components of the OCR portion in the "IPA1-OCR" pair, respectively; β_cis and β_trans are the fixed-effect regression coefficients of CA_cis and CA_trans on the gene expression level y, respectively; μ is the intercept term; u_G is the random effect used to correct for population structure, and its covariance matrix is ​​given by the whole-genome kinship matrix K_G; ε is the model residual.

[0043] Causal determination: For the same "gene IPA1-OCR" pair, when the P-value of its OCR portion in both of the above models is less than a preset statistical significance threshold (for example, set to 1×10 in this embodiment), -4 If the regression coefficients β_cis and β_trans have the same sign, then the OCR is determined to be the causal cis-regulatory region of the gene.

[0044] If both regression coefficients β_cis and β_trans are positive, it is a positive control region; if both are negative, it is a negative control region.

[0045] 2. Sorting: Calculate the z-score statistics z_cis and z_trans for the two models, where z_cis is the statistic of the cis component in the cis-channel model and z_trans is the statistic of the trans component in the trans-channel model, and then apply the formula S = (z_cis²) The comprehensive score S-score is calculated using z_trans²) / (z_cis² + z_trans²) and used to rank multiple causal cis-regulatory regions of the IPA1 gene.

[0046] 3. Output a target list including regulatory region coordinates, regulatory direction, statistical significance, and priority score, and generate a specific genome editing design scheme based on this list to achieve the desired regulation of target gene expression levels. The design scheme includes: (1) To achieve upregulation of gene expression: strategies can be designed to disrupt or weaken the function of an identified negative regulatory region, or strategies can be designed to enhance the function of an identified positive regulatory region.

[0047] (2) To achieve downregulation of gene expression: strategies can be designed to disrupt or weaken the function of an identified positive regulatory region, or strategies can be designed to enhance the function of an identified negative regulatory region.

[0048] It should be noted that the above dual-channel correlation analysis was performed on all candidate windows one by one at a high-resolution level using a "sliding window" approach. After completing the correlation analysis on all windows, to facilitate global statistics (such as... Figure 1 As shown in the diagram, and in addition to the explanation of biological functions, this invention also includes a result merging step. Specifically, multiple sliding windows that are physically adjacent or overlapping on the genome and all show significant causal association signals for the same gene are merged and defined as a single functional unit, namely a causal cis-regulatory region or a signal peak. In subsequent statistical or case validation (such as...),... Figure 2 Typically, the window with the strongest correlation signal (i.e., the highest S-score) within the causal cis-regulatory region or peak is selected as the "representative window" of the functional unit.

[0049] IV. Experimental Results 1. The results are as follows Figure 2 As shown, a high-confidence distal causal OCR (dOCR) was identified 7.2 kb upstream of the TSS of the IPA1 gene. Figure 2 In Figure A, in the upper orbital plot, each orange dot represents a sliding window, and its Y-axis position indicates the causal association strength (S-score) calculated using the method of this invention. Multiple significant windows cluster together to form a signal peak, collectively identified as a distal causal OCR (dOCR). The lower orbital plot shows the genome background, illustrating the original ATAC-seq signal, the physical location of the OCR, and... IPA1 Gene structure.

[0050] exist Figure 2 B and Figure 2 In C, the X-axis of the scatter plot represents the chromatin accessibility of a representative window of the dOCR (i.e., the red dots in Figure A), where... Figure 2 B is the cis component (CAcis). Figure 2 C represents the trans component (CAtrans); the Y-axis is... IPA1 Gene expression levels. The results showed that both exhibited a significant and consistent positive correlation (linear mixture model LMM). β >0, P <1×10 -4 This satisfies the causal determination criterion of the present invention, thus proving that the dOCR is... IPA1 It is a positive regulatory element of IPA1. Based on this, this embodiment determines that the dOCR is a positive regulatory element (enhancer) of IPA1, and proposes the hypothesis that: disrupting this element will lead to a decrease in IPA1 expression, thereby causing a reduction in the number of grains per ear, and may increase the number of tillers.

[0051] 2. Functional verification of computational predictions was performed through gene editing experiments.

[0052] To verify the above predictions, this embodiment, following the method in step three, utilizes CRISPR / Cas9-mediated genome editing technology to design target sites in the core region of the dOCR. These target sites are located 6890 bp and 7367 bp upstream of the TSS of the IPA1 gene. Figure 3 In diagram A, the upper part is the reference genome, indicating the target sites of the two gRNAs (gRNA-1: AGAGCCTCCATATCTCAGTT, gRNA-2: TGTGGCGGCACGGCAAAGCT). The lower part shows the obtained homozygous edited alleles ( IPA1 CR-dOCR This resulted in deletions of 3 bp and 36 bp in the target region, respectively. Successfully obtained the name... IPA1 CR-dOCR The homozygous mutant strain carries two deletions, 3 bp and 36 bp, respectively, in the target region.

[0053] The results are as follows Figure 3 B~ Figure 3 As shown in E, this homozygous mutant IPA1 CR-dOCR Phenotypic analysis showed that, compared with wild-type (WT) plants, IPA1 CR-dOCR Homozygous plants showed a significant increase in the number of effective panicles (24.9%, P = 0.0016), while the number of grains in the main panicle decreased significantly (27.4%, P = 8.87 × 10⁻⁶). -10 This phenotypic change (increased tillering and reduced grain number per ear) is consistent with the known phenotype of weakened IPA1 function. The experimental results strongly demonstrate that disrupting this predicted positive regulatory element does indeed lead to the corresponding phenotype of weakened gene function, thus verifying the causality of the predicted regulatory region and the accuracy of its positive regulatory direction as determined by the method of this invention.

[0054] in conclusion: The experimental results strongly confirm that disrupting the distal positive regulatory element predicted by the method of this invention does indeed lead to the expected phenotype consistent with weakened gene function. This not only verifies the causal and positive regulatory function of this dOCR on IPA1, but more importantly, it serves as a representative case that powerfully demonstrates the accuracy and strong application value of the overall method of this invention in accurately identifying functional regulatory elements at the single-gene level.

[0055] Example 2 Based on the method in Example 1, the remaining 30,813 expressed genes obtained in step one were subjected to high-precision localization of causal cis-regulatory regions and determination of regulatory direction.

[0056] Comparative Example Direct association analysis method This direct association analysis method does not perform cis / trans decomposition. Its specific implementation steps are as follows: Step 1: Construct the analysis object.

[0057] For each expressed gene obtained in step one of Example 1, based on the candidate cis-regulatory range defined in Example 1 (i.e., 100 kb upstream and downstream of the transcription start site TSS), all corresponding candidate OCRs are extracted from the OCR accessibility matrix finally generated in Example 1. Each gene and each candidate OCR are paired to form a "gene-OCR" pair, which serves as the basic unit for subsequent association analysis.

[0058] Step 2: Construct and fit the correlation model.

[0059] For each gene-OCR pair formed, the following mixed linear model is constructed and fitted using the undecomposed chromatin accessibility value of the OCR (denoted as x1) and the expression level of the corresponding gene (denoted as y1): y1 = μ1 + β x1 + u_G + ε; Where y is the expression level of the gene portion in the "gene-OCR" pair, and its value comes from the gene expression matrix after the final processing in Example 1; x1 is the undecomposed accessibility value of the OCR portion in the "gene-OCR" pair, and its value comes from the OCR accessibility matrix after the final processing in Example 1; μ is the intercept term; β is the fixed-effect regression coefficient of x1 on y1, representing the direct association effect of OCR accessibility on gene expression; u_G is the random effect used to correct for population structure, and its covariance matrix is ​​given by the whole genome kinship matrix K_G constructed in step two of Example 1; ε is the model residual.

[0060] Step 3: Identify significant associations.

[0061] For each "gene-OCR" pair in step two, extract the p-value corresponding to the fixed effect coefficient β. Compare this p-value with a preset significance threshold (1×10⁻⁶). -4 The values ​​are compared. When the P-value is less than the threshold, it is determined that the OCR has a significant direct association with the gene.

[0062] Example 3: Result Comparison and Analysis To highlight the advantages of the method of the present invention, the results of Example 2 and the results of the comparative example are analyzed.

[0063] The results obtained through the comparative example will be compared and analyzed with the results obtained in Example 2. The results are as follows: Figure 1 As shown, the comparison clearly highlights the significant advantages of this invention in improving the specificity, reliability, and remote identification capability of regulatory relationship identification.

[0064] 1. Specificity comparison First, regarding the specificity of the association, the direct association analysis method identified a high proportion of genes with potential cis-associated OCRs, reaching 73.2%. Figure 1 A, corresponding to column "CA" in the figure), and 45.8% of the genes were associated with more than 5 potential OCRs ( Figure 1 B). This broad and non-specific result contains a large number of non-causal false positives caused by linkage disequilibrium or co-accessibility, making it difficult to focus subsequent functional validation. In contrast, the method of this invention uses a rigorous causal inference framework to precisely filter the proportion of associated genes to 20.5% (B). Figure 1 A, corresponding to column "CAcis & trans" in the figure), of which only 6.7% of genes are associated with more than 5 causal OCRs ( Figure 1 B). This clearly demonstrates that the present invention can effectively filter out the vast majority of false positive signals, greatly improving the specificity of regulatory relationship identification.

[0065] 2. Reliability and Functional Dependency Verification Secondly, to assess the biological reliability of the causal OCR identified in this invention, its enrichment around the transcription start site (TSS) was analyzed. For example... Figure 1 As shown in Figure C, the causal OCR (red curve) determined by the method of this invention forms a much sharper and stronger enrichment peak near the TSS than the comparative figure (cyan curve, "CA" curve). This is highly consistent with the general understanding in the field of biology (that key cis-regulatory elements are likely distributed near the TSS of genes). Therefore, this strong TSS enrichment feature strongly confirms from the perspective of functional relevance that the OCR identified by this invention is a real, functional regulatory element, rather than random statistical noise, thus demonstrating the high reliability of the results of the method of this invention.

[0066] 3. Demonstration of remote control element identification capability Based on the proven high specificity and high reliability, this invention demonstrates a unique ability to identify remote control elements. For example... Figure 1 As shown in Figure D, spatial distribution analysis of the high-confidence causal OCRs identified in this invention revealed that a significant proportion (28.1%) of the elements were located in distant regions far from the gene. This is because... Figure 1C confirms the overall reliability of the results of this invention. Therefore, this reliability basis allows the invention to determine that these identified remote signals are genuine long-range causal control elements, rather than illusions that cannot be distinguished from background noise in traditional methods. This demonstrates that the present invention effectively solves the major technical problem of the difficulty in accurately identifying remote control elements due to signal attenuation and noise interference in traditional methods.

[0067] In summary, the method of this invention represents a leap from "broad and fuzzy" association analysis to "precise and reliable" causal inference. It not only liberates researchers from a massive amount of potential correlation signals, accurately pointing to a few highly reliable candidate regulatory elements, but also empowers researchers to identify and trust distant regulatory elements, effectively avoiding the predicament of "trial and error" amidst a large number of false positives, thus solving a core pain point in the field.

[0068] All other parts not described in detail are existing technologies. Although the above embodiments have provided a detailed description of the present invention, they are only some embodiments of the present invention, not all embodiments. People can obtain other embodiments based on these embodiments without creative effort, and these embodiments all fall within the protection scope of the present invention.

Claims

1. A method for molecular navigation and directional design of gene regulatory regions in agricultural organisms, characterized in that: Includes the following steps: S1. Obtain genome-wide gene expression profiles, chromatin accessibility maps, and genome-wide genotype data of a related population of a certain agricultural organism; S2. Based on whole-genome gene expression profile data, the expression level of each gene in each sample of the population is quantified, and a gene expression matrix is ​​constructed. S3. Based on chromatin accessibility map data, identify open chromatin regions across the entire genome, quantify the accessibility signal intensity of each open chromatin region in each sample of the population, and construct an OCR accessibility matrix. S4. Using the transcription start sites of all genes in the gene expression matrix as anchors, define the upstream and downstream linear windows of each gene as its candidate cis-regulatory range. For each target gene, extract all OCRs whose coordinates lie within its candidate cis-regulatory range from the OCR accessibility matrix, and use them as candidate regulatory regions for that gene. S5. Perform whole-genome filtering on the whole-genome genotype data and construct a whole-genome kinship matrix K_G; S6. For each candidate regulatory region of a gene, construct a local kinship matrix K_local using the genetic variations in the genomic regions adjacent to the candidate regulatory region in the whole genome genotype data; S7. Chromatin accessibility analysis was performed on the candidate regulatory regions of each gene based on the whole genome kinship matrix and the local kinship matrix; and cis-heritability testing was performed on the candidate regulatory regions of each gene to screen out the OCRs of each gene that "have significant cis-heritability". S8. For each OCR that "has significant cis-heritability" selected in step S7, apply the statistical model used in the cis-heritability test to decompose the chromatin accessibility of the OCR into the cis component CA_cis determined by local genetic variation and the trans component CA_trans determined by other factors. S9. For each gene and its corresponding OCR that "has significant cis-heritability", form a "gene-OCR" pair, and perform a dual-channel association analysis of cis and trans channels for each "gene-OCR" pair: Cis-channel: The model is in the form of y = μ + β_cis CA_cis + u_G + ε; Inverse channel: The model is in the form of y = μ + β_trans CA_trans + u_G + ε; Where y is the expression level of the gene component in the "gene-OCR" pair; CA_cis and CA_trans are the cis and trans components of the OCR component in the "gene-OCR" pair, respectively; β_cis and β_trans are the fixed-effect regression coefficients of CA_cis and CA_trans on the gene expression level y, respectively; μ is the intercept term; u_G is the random effect used to correct for population structure; and ε is the model residual. S10. Based on the above dual-channel correlation analysis, the following judgment is made: When the P-value of the OCR portion in the same "gene-OCR" pair is less than the statistical significance threshold in both of the above models, and its regression coefficients β_cis and β_trans have the same sign, the OCR is determined to be the causal cis-regulatory region of the gene. If both regression coefficients β_cis and β_trans are positive, it is a positive control region; if both regression coefficients β_cis and β_trans are negative, it is a negative control region.

2. The method according to claim 1, characterized in that: It also includes the following steps: (1) Calculate the z-score statistics z_cis and z_trans corresponding to the two models of the cis-channel and trans-channel, and calculate the comprehensive score. Sort the multiple causal cis-regulatory regions of a gene according to the comprehensive score. (2) Output a list of targets including the coordinates of the regulatory region, the direction of regulation, the statistical significance and the priority score, and generate a specific genome editing design scheme based on it to achieve the expected regulation of the expression level of the target gene.

3. The method according to claim 2, characterized in that: In step S1, the agricultural organism is any one of rice, corn, wheat, barley, sorghum, soybean, rapeseed, cotton, tomato, and potato. The number of members in the associated group is not less than 200; Whole-genome gene expression profile data were obtained through RNA-seq sequencing; Chromatin accessibility map data were obtained through ATAC-seq sequencing; whole-genome genotype data were obtained through whole-genome resequencing.

4. The method according to claim 3, characterized in that: The agricultural organism in question is rice; In step S2, the expression level of each gene in each sample is quantified, and the specific steps for constructing the gene expression matrix are as follows: S21. Using rapid quantitative tools, RNA-seq sequencing reads are quantified based on the annotated transcript sequences of the target species' reference genome to count the number of reads for each gene. S22. Calculate the TPM of each gene and combine the TPM values ​​of different transcripts of the same gene as the total TPM of the gene. Perform log2(TPM+1) transformation on the total TPM value of the gene and screen genes expressed in at least 5% of the samples for downstream analysis to obtain the original gene expression matrix. S23. The original gene expression matrix is ​​corrected using a latent factor model, and the corrected expression values ​​of each gene in all samples are subjected to a rank-based inverse normal transformation to obtain the gene expression matrix.

5. The method according to claim 3, characterized in that: In step S3, the specific steps for identifying open chromatin regions across the entire genome and quantifying the accessibility signal intensity of each open chromatin region in each sample to construct the OCR accessibility matrix are as follows: S31. After aligning the ATAC-seq sequencing reads to the reference genome of the target species, the alignment results are filtered and refined, and high-quality aligned sequencing fragments are selected. Then, Tn5 transposon cleavage site offset correction is performed. S32. Using a 250 bp sliding window with a step size of 100 bp, the adjusted transposon cut sites are counted across the entire genome to quantify chromatin accessibility and obtain the original "window × sample" counting matrix. S33. Perform quantile normalization on the original "window × sample" count matrix; S34. Based on the normalized signal value, calculate an index representing the overall accessibility level of each window, sort all windows based on this index, and select the 20% of windows with the strongest signal as open chromatin regions. S35. For the selected open chromatin regions, their normalized count values ​​are logarithmically transformed after adding a pseudo-count of 1. S36. The latent factor model is used to remove the influence of latent confounding factors, and the accessibility value of each corrected open chromatin region is subjected to a rank-based inverse normal transformation to obtain the OCR accessibility matrix.

6. The method according to claim 2, characterized in that: In step S4, the physical distance of the linear window is 5kb to 1000kb; In step S5, the condition for whole-genome filtering is that the minor allele frequency (MAF) > 0.

05. In step S6, the range of adjacent genomic regions is 5 kb to 1000 kb.

7. The method according to claim 2, characterized in that: In step S7, the specific steps for chromatin accessibility analysis are as follows: S71. For each candidate regulatory region of a gene, fit the following variance component model to analyze its accessibility: x = μ + u_cis + u_G + ε; Where x is chromatin accessibility, μ is the intercept term, u_cis follows a distribution of N(0, σ²_cis K_local), representing a cis-random effect determined by local genetic variation, u_G follows a distribution of N(0, σ²_G K_G), representing a random effect determined by the whole genome genetic background and used to control population structure and background polygenic effects, and ε follows a distribution of N(0, σ²_e I), representing the residual; The specific steps for screening each gene for "significant cis-heritability" using OCR are as follows: S72. Obtain the test P value of the cis variance component σ²_cis for each candidate regulatory region of each gene, perform multiple test correction on all P values, and define the candidate regulatory region with FDR < 0.05 as "having significant cis heritability" OCR; In step S8, the specific steps for obtaining the cis component CA_cis and the trans component CA_trans are as follows: S81. For the OCR that "has significant cis-heritability" in step S72, extract the best linear unbiased predictor value of its cis-random effect from the variance component model in step S71. This value is defined as the cis-component CA_cis of the OCR accessibility. Calculate the residual x. CA_cis, the residual is defined as the trans component CA_trans; The other factors refer to all factors other than local genetic variations.

8. The method according to claim 7, characterized in that: In step S10, the statistical significance threshold is 1×10⁻⁶. -4 .

9. The method according to claim 2, characterized in that: In step (1), the comprehensive score S is calculated according to the following formula: S = (z_cis² z_trans²) / (z_cis² + z_trans²)。 10. The application of the method of claim 1 in agricultural biological genetic improvement.

Images