Cow ketosis regulatory gene prediction system based on multi-omics analysis and machine learning

Abstract

The invention provides a dairy cow ketosis regulatory gene prediction system based on multi-omics analysis and machine learning, which belongs to the technical field of molecular breeding and disease prevention and control, and comprises a data acquisition and preprocessing module, a feature set establishment module, a machine learning model establishment module and a result evaluation module, the method comprises the following steps: integrating dairy cow genome, transcriptome and metabolome data, screening candidate regulatory genes through whole genome association analysis, gene differential expression analysis, cis-eQTL positioning and co-positioning analysis, and constructing a gene expression feature set; core regulation genes are screened through L1 regularization penalty by means of a Lasso model, weights are distributed, model parameters are optimized in combination with grid search and cross validation, and model performance is evaluated through an ROC curve and an AUC value. According to the invention, an integrated technical system from gene screening to risk prediction is constructed, efficient screening of the core regulation gene and accurate prediction of ketosis risk are realized, and the economic loss of breeding is effectively reduced.

Description

Technical Field

[0001] This invention relates to the fields of molecular breeding and disease control technology, and in particular to a prediction system for regulatory genes of ketosis in dairy cows based on multi-omics analysis and machine learning. Background Technology

[0002] Bovine ketosis is a prevalent metabolic disease in peripartum dairy cows, primarily caused by abnormal fat catabolism resulting from a negative energy balance. Affected cows exhibit symptoms such as decreased appetite, reduced milk production, and weight loss. In severe cases, it can lead to abortion, infertility, and even death. It can also induce complications such as mastitis and metritis, causing significant economic losses to the dairy farming industry. Currently, related research focuses primarily on clinical diagnosis and nutritional regulation, with relatively lagging research into the molecular mechanisms of ketosis development. Regulatory genes play a crucial role in the pathogenesis of ketosis. Screening and identifying relevant regulatory genes is a core prerequisite for revealing the pathogenesis, developing early diagnostic biomarkers, and formulating precise prevention and control measures. Furthermore, early prediction of ketosis risk based on regulatory genes is a key means to reduce morbidity and improve farming efficiency.

[0003] In existing technologies, gene screening mainly relies on single methods such as genome-wide association analysis (GWAS) and expression quantitative trait locus (eQTL) analysis, while disease prediction largely depends on clinical phenotypic indicators. However, these technologies have significant limitations: GWAS struggles to effectively distinguish between causal variants and linkage disequilibrium variants, leading to insufficient accuracy in gene screening; traditional disease prediction methods are easily affected by external factors such as breeding environment and management level, resulting in limited prediction accuracy. Although colocalization analysis can integrate multi-omics data and improve the accuracy of causal gene screening, and Lasso models can combine feature screening and predictive modeling through L1 penalty, possessing unique advantages in biological data mining, there is currently no technology to integrate colocalization analysis into database modules, nor is there an integrated system that uses Lasso models to first screen regulatory genes and then predict ketosis risk. This results in existing technologies being unable to efficiently complete the collaborative process of "gene screening - risk prediction," making it difficult to meet the actual needs of precision prevention and control. Summary of the Invention

[0004] The purpose of this invention is to provide a prediction system for regulatory genes of ketosis in dairy cows based on multi-omics analysis and machine learning. By deeply integrating colocalization analysis into the feature set construction module, and by integrating multi-omics data such as genomics and transcriptomics, combined with GWAS, eQTL analysis and colocalization analysis, it solves the pain point of single gene screening methods in distinguishing causal variations, significantly improving the screening accuracy of candidate regulatory genes. At the same time, it adopts the Lasso model, and through L1 regularization penalty, it achieves the synergy of invalid gene elimination, core gene weight allocation and prediction model construction, breaking through the limitations of traditional prediction methods that are susceptible to environmental interference, and realizing an efficient collaborative process from regulatory gene screening to ketosis risk prediction.

[0005] To achieve the above objectives, this invention proposes a bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning, comprising a data acquisition and preprocessing module, a feature set establishment module, a machine learning model establishment module, and a result evaluation module connected in sequence. These modules work together to predict and validate bovine ketosis regulatory genes, wherein: The data acquisition and preprocessing module includes a phenotype collection unit, a genome data processing unit, and a transcriptome data processing unit. It collects phenotype data of bovine ketosis, outputs a ketosis-related dataset, processes multi-omics data, and outputs a basic data sub-library after quality control, comparison, and quantitative analysis. The feature set construction module includes a genome-wide association analysis unit, a gene differential expression analysis unit, a cis-eQTL analysis unit, and a colocalization analysis unit; The machine learning model building module includes a regulatory gene prediction unit, which constructs a Lasso model based on screened candidate regulatory genes, selects core regulatory genes and assigns weights through L1 regularization penalty, and optimizes model parameters by combining grid search and cross-validation. The results evaluation module uses ROC curves and AUC values ​​to assess model overfitting and the effectiveness of ketosis prediction.

[0006] Preferably, the multi-omics data includes bovine genomic data, transcriptomic data, and metabolomic data, wherein: the genomic data contains SNP site information; the transcriptomic data contains gene expression level information; the metabolomic data contains blood BHBA content information; and the phenotypic data includes the bovine disease status, which is labeled with a binary label: diseased = Case, healthy = Control.

[0007] Preferably, the classification criteria for BHBA concentration are as follows: BHBA concentration ≤ 1.4 mmol / L is considered healthy, 1.4 mmol / L < BHBA ≤ 3 mmol / L is considered subclinical ketosis, and BHBA > 3 mmol / L is considered clinical ketosis. Subsequent phenotypic labeling and dataset construction are carried out based on the ketosis classification results.

[0008] Preferably, the quality control criteria for genomic data are: a pass rate of less than 95% for removing individual genotypes, a pass rate of less than 90% for removing individual SNP genotypes, a minor allele frequency (MAF) of the SNP < 0.05, and a deviation from Hardy-Weinberg equilibrium value (P < 1.0 × 10⁻⁶). -6 The mark.

[0009] Preferably, the processing flow of transcriptome data in the data acquisition and preprocessing module is as follows: after quality control of the sequencing data, it is aligned to the reference genome ARS-UCD2.0, quantitative analysis is performed, and the expression matrix of the annotated genes is obtained.

[0010] Preferably, in the feature set establishment module, the analysis results of the genome-wide association analysis unit, gene differential expression analysis unit, and cis-eQTL analysis unit are output as an intermediate data sub-library. The co-localization analysis results performed using the genome-wide association analysis results and cis-eQTL results, combined with the gene differential expression analysis results, are output as a gene expression feature set, and the candidate genes selected are stored in the gene expression feature set.

[0011] Preferably, the common location analysis steps in the feature set establishment module are as follows: Step S1: Perform basic data organization, organize the basic information of the sample, including but not limited to pathological status and gene expression levels; Step S2: Perform genome-wide association analysis, using the ketosis-associated phenotype as the dependent variable and the pre-processed SNP loci as the independent variable, calculate the association strength P-value between each SNP locus and the ketosis phenotype, and screen for significantly associated SNP loci. Step S3: Using the ketosis classification phenotype as the dependent variable and gene expression level as the independent variable, screen for TOP genes with differential expression under different ketosis states; Step S4: Perform cis-eQTL analysis, with the difference in gene expression levels as the dependent variable and the pre-processed SNP sites as the independent variables. Set the distance threshold of cis-expression quantitative trait loci to 1 MB and screen for significant cis-eQTLs associated with gene expression changes. Step S5: Perform colocalization analysis, integrate GWAS results with matched eQTL information for association verification, calculate posterior probability by setting hypotheses, and quantify the possibility of SNP sites simultaneously regulating ketosis phenotype and gene expression. Step S6: Determine if there is co-localization of SNP sites, record the corresponding SNP sites as target SNP sites that simultaneously affect the ketosis phenotype and gene expression in dairy cows, and select the genes associated with the target SNP sites as candidate regulatory genes. Step S7: Establish a gene expression feature set, integrate GWAS, cis-eQTL and colocalization analysis results, and construct the gene expression feature set through gene annotation, metabolic pathway and functional enrichment analysis.

[0012] Preferably, in the machine learning model, the working steps of the gene prediction unit are as follows: Step T1: Divide the feature set into a training set and a validation set using random sampling to ensure that the two sets of data have the same distribution. The training set is used for model training, feature selection and parameter optimization, and the validation set is used to evaluate model performance. Step T2: Construct a Lasso regression model based on the training set. The model formula is as follows: ; ; in, The set of vectors representing the possible weights of each feature gene. For sample size, For sample number, For the first The true labels for each sample are: diseased = Case, healthy = Control. The model predicts the first The probability that a sample is a positive class. For L1 penalty intensity, The number of characteristic genes, For the first The regression coefficient of each gene, The intercept is... For the first Gene expression vectors of each sample, It is the transpose symbol; Step T3: Feature selection of candidate regulatory genes is performed using L1 regularization penalty. The penalty formula is as follows: ; in, This is an L1 penalty item; Step T4: Automatically compress the weights of genes that do not contribute to the regulation of ketosis to 0 and remove them. Optimize the model regularization parameters using a grid search method and five 10-fold cross-validation iterations. The formula for the grid search method is: ; in, For the set of candidate regularization parameters, For the number of grid points, Number the grid points. The L1 penalty intensity is the value corresponding to the minimum mesh size. The L1 penalty intensity is the value corresponding to the maximum grid size. For linearly equidistant cut points, ; The formula for 5-fold 10-fold cross-validation is: ; in, For health status labels, For a given The average cross-validation performance, For the number of repetitions, For folds, For the first The repetition, the... The AUC value of the fold; Step T5: Select the optimal L1 penalty intensity, the selection rule is as follows: ; in, For optimal penalty intensity, For a given The average AUC obtained on cross-validation, To obtain To achieve the maximum intensity of punishment.

[0013] Preferably, the results evaluation module systematically evaluates the output of the machine learning model building module, focusing on two core objectives: model overfitting and the effectiveness of ketosis prediction. ROC curves and AUC values ​​are used as the core evaluation metrics. Specific steps include: Step W1: Evaluate the model's effectiveness in overfitting and predicting ketosis. Plot ROC curves and calculate AUC values ​​based on the training and validation sets respectively. The ROC curve formula is: ; ; in, The current threshold, , For the first The model predicts the probability for each sample. For the first The true labels for each sample are: diseased = Case, healthy = Cnotrol. For indicator functions, when The value is 1 if the condition is met, otherwise it is 0. The formula for calculating AUC is: ; in, For the horizontal axis variable, Given an FPR value, return the corresponding threshold. , Threshold The actual rate, Integrate along FPR from 0 to 1; Step W2: Compare the AUC values ​​of the training set and the validation set to determine whether overfitting has occurred, and evaluate the effectiveness of ketosis prediction based on the AUC values ​​of the validation set.

[0014] Preferably, the evaluation criteria for the results evaluation module are as follows: Overfitting assessment: A difference of ±0.1 between the AUC values ​​of the training set and the validation set indicates that the model has strong generalization ability; Predictive effectiveness assessment: Comparing the ROC curves and AUC values ​​of the training and validation sets, a validation set AUC value ≥ 0.8 indicates that the core regulatory genes can effectively predict ketosis.

[0015] Therefore, this invention proposes a prediction system for regulatory genes of ketosis in dairy cows based on multi-omics analysis and machine learning, with the following beneficial effects: (1) This invention integrates multi-omics data from the genome, transcriptome, and metabolome, and combines GWAS, gene differential expression analysis, cis-eQTL localization and colocalization analysis to effectively distinguish between causal variation and linkage disequilibrium variation. The candidate regulatory genes screened are more targeted, which solves the problem of insufficient accuracy of single gene screening methods.

[0016] (2) The present invention uses the Lasso model to remove invalid genes through L1 regularization penalty, and combines grid search and 5 times 10-fold cross-validation to optimize parameters, which significantly improves the accuracy of ketosis prediction and is not easily affected by environmental factors.

[0017] (3) This invention integrates colocalization analysis into the feature set building module for the first time, and combines the Lasso model to build an integrated system of "core regulatory gene screening - ketosis risk prediction". No step-by-step operation is required, which greatly improves the analysis efficiency and meets the actual needs of precise prevention and control.

[0018] (4) The core regulatory genes screened by this invention provide key targets for revealing the molecular mechanism of ketosis in dairy cows. Risk prediction based on core genes can realize early warning of ketosis, providing technical support for molecular breeding, precise nutrition regulation and disease prevention and control in dairy cows, and significantly reducing economic losses in breeding. Attached Figure Description

[0019] Figure 1 The structure and overall workflow of a prediction system for regulatory genes of ketosis in dairy cows based on multi-omics analysis and machine learning are shown in the diagram. Figure 2 A schematic diagram of the kernel density of AUC values ​​for grid search and 5-fold 10-fold cross-validation; Figure 3 This is a schematic diagram of the ROC curve for the training set. Figure 4 This is a schematic diagram of the ROC curve for the validation set. Detailed Implementation

[0020] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0021] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning as understood by one of ordinary skill in the art to which this invention pertains.

[0022] Example 1 like Figure 1 As shown, this invention provides a prediction system for regulatory genes of bovine ketosis based on multi-omics analysis and machine learning, comprising a data acquisition and preprocessing module, a feature set establishment module, a machine learning model establishment module, and a result evaluation module connected in sequence, wherein: The data acquisition and preprocessing module includes a phenotype collection unit, a genome data processing unit, and a transcriptome data processing unit. It collects phenotype data of bovine ketosis, outputs a ketosis-related dataset, processes multi-omics data, and outputs a basic data sub-library after quality control, comparison, and quantitative analysis.

[0023] The feature set construction module includes a genome-wide association analysis unit, a gene differential expression analysis unit, a cis-eQTL analysis unit, and a colocalization analysis unit. This module is used to manage various types of data and integrates colocalization analysis functions to realize association mining of multi-omics data and preliminary screening of candidate regulatory genes.

[0024] The machine learning model building module includes a regulatory gene prediction unit, which constructs a Lasso model based on screened candidate regulatory genes, selects core regulatory genes and assigns weights through L1 regularization penalty, and optimizes model parameters by combining grid search and cross-validation. The results evaluation module uses ROC curves and AUC values ​​to assess model overfitting and the effectiveness of ketosis prediction.

[0025] The data acquisition and preprocessing modules collect high-quality raw data and standardize it to provide a reliable data foundation for subsequent analysis. Specifically, the phenotypic collection unit collects postpartum blood samples from Holstein dairy cows, detects BHBA concentrations, and labels them with disease tags (Case / Control) to clarify ketosis classification (healthy, subclinical, clinical), providing a basis for disease phenotypic labeling. The genomic data processing unit performs quality control on whole-genome sequencing or SNP microarray data, removing low-pass-rate, low-minor-allelic-frequency, and markers deviating from Hardy-Weinberg equilibrium, while retaining valid SNP loci information. The transcriptome data processing unit, after quality control of sequencing data, aligns it to a reference genome (ARS-UCD2.0) and generates annotated gene expression matrices through quantitative analysis to ensure the accuracy of gene expression data. The final output includes a ketosis-related dataset (containing phenotypic tags and BHBA concentrations) and a basic data sub-library (quality-controlled SNP loci and gene expression matrices).

[0026] The feature set construction module integrates multi-omics data for association mining, screening candidate regulatory genes related to ketosis from a massive gene pool to construct a high-quality feature set. It integrates information such as sample pathological status, gene expression levels, and SNP sites to construct a standardized data matrix, providing a unified data format for subsequent analyses. The genome-wide association analysis unit uses ketosis phenotype as the dependent variable and SNP sites as independent variables to screen for significantly associated SNP sites, identifying genetic variation regions related to the ketosis phenotype. The differential gene expression analysis unit compares gene expression levels under different disease states, screening for differentially expressed genes and focusing on phenotype-related functional genes. The cis-eQTL analysis unit uses gene expression level as the dependent variable and SNP sites as independent variables, setting a 1MB cis-distance threshold to screen for significant cis-eQTLs affecting gene expression, revealing the regulatory relationship between genetic variations and gene expression. The colocalization analysis unit integrates GWAS and eQTL results, calculates posterior probabilities using the Coloc algorithm, and screens causal SNP-related genes that simultaneously regulate ketosis phenotype and gene expression, improving the specificity of candidate genes. The output includes intermediate data sub-libraries (GWAS, differential expression analysis, cis-eQTL localization results) and gene expression feature sets (expression level data of candidate regulatory genes).

[0027] The machine learning model building module aims to achieve integrated analysis of core regulatory gene screening and ketosis risk prediction using the Lasso model, transforming candidate genes into a practical prediction model. Random sampling is used to divide the gene expression feature set into training and validation sets, ensuring consistent data distribution and guaranteeing the objectivity of model validation. A regression model is built based on the training set, automatically eliminating genes with no regulatory effect (weights compressed to 0) through L1 regularization penalty, and assigning non-zero weights to effective genes. A grid search combined with five 10-fold cross-validations is used to select the optimal regularization parameters. This improves the model's generalization ability; it outputs genes with non-zero weights from the model and predicts the risk of ketosis on the validation set based on core genes. The output includes a list of core regulatory genes (including weights), the trained Lasso prediction model, and the prediction results on the validation set.

[0028] The results evaluation module scientifically assesses the reliability and predictive performance of the model, providing a basis for the practicality of the technical solution. It plots ROC curves and calculates AUC values ​​based on the training and validation sets. By comparing the differences in AUC values ​​between the two sets, it determines whether the model is overfitting (a difference within ±0.1 indicates strong generalization ability). Using a validation set AUC ≥ 0.8 as the criterion, it verifies the effectiveness of ketosis prediction based on core regulatory genes. The final output is a model evaluation report, including ROC curves, AUC values, overfitting assessment conclusions, and prediction effectiveness conclusions.

[0029] The modules work together to predict and validate genes regulating ketosis in dairy cows, specifically including the following steps: Step 1: Data collection and preprocessing, including the following steps: Step 1.1: Characterization and Data Collection: Blood samples were collected from the tail vein of 1030 Holstein dairy cows postpartum. BHBA concentration in the blood samples was measured using a handheld FreeStyle Optium Neo H ketometer. A BHBA concentration ≤ 1.4 mmol / L was considered healthy; 1.4 mmol / L < BHBA ≤ 3 mmol / L was considered subclinical ketosis; and BHBA > 3 mmol / L was considered clinical ketosis. Genetic assessment was performed using ssGBLUP. A mixed model was developed using the GLM procedure implemented in R4.3.2, with blood BHBA as the dependent variable, to determine the fixed effects in the model. The model is as follows: ; in, Blood BHBA concentration, For fixed-effect field - year - quarter (Level 104) Indicates parity (1, 2, 3+) Use the number of days of lactation as a covariate. It is an additive genetic effect. For permanent environmental effects, As random residuals, the breeding value was used as the dependent variable for subsequent genome-wide association analysis.

[0030] Step 1.2: Genomic Data Processing: Genomic DNA was extracted from bovine blood using the phenol-chloroform method. DNA quality was assessed using UV spectrophotometry and gel electrophoresis. Genotyping was performed using the Neogen Biotechnology Geneseek Genomic Profiler (GGP) Bovine 100 K SNP chip. Quality control (QC) was performed using Plink software, removing genotypes with a pass rate below 95%, single SNP genotype pass rate below 90%, minor allele frequency (MAF) of SNPs < 0.05, and those deviating from Hardy-Weinberg equilibrium (P < 1.0 × 10⁻⁶). -6 The mark of ).

[0031] Blood samples from 50 dairy cows, including blood ketone levels at 5, 7, and 12 days postpartum, were sent to Shijiazhuang Borui Biotechnology Co., Ltd. for testing. Resequencing libraries were constructed from the qualified DNA using the GenoBaits® DNA Library Prep Kit. These libraries were then sequenced using the BGI DNBSEQ-T7 / Illumina NovaSeq 6000 sequencing platform in PE150 mode, yielding approximately 7X genotypic coverage. MultiQC software was used for quality assessment, removing low-quality bases and adapter sequences. BWA was used to align to the reference genome ARS-UCD2.0, and finally, 12,354,736 SNPs were retained for cis-eQTL analysis.

[0032] Step 1.3: Transcriptome Data Processing: Total RNA was extracted from each sample using Trizol reagent. RNA integrity and DNA contamination were analyzed using 1% agarose gel electrophoresis. RNA purity (OD260 / 280 and OD260 / 230 ratios) and concentration were detected using a NanoPhotometer spectrophotometer and a Qubit 2.0 fluorometer. RNA integrity was assessed using an Agilent Bioanalyzer 2100 system, ensuring that the absorbance (OD260 / 280) of all RNA samples was between 1.80 and 2.06, and that an RNA integrity score greater than 8 was used to verify RNA quality. Illumina was used for further analysis. The USA sample preparation kit was used to construct a cDNA library for each sample using approximately 1 μg of RNA. The mRNA was first synthesized into cDNA in a reverse transcriptase system, then randomly fragmented, and the fragmented cDNA was used as a template. Finally, cDNA with a length of 250-300 bp was selected, amplified, and purified to obtain the final library. The obtained transcriptome sequencing data were quality controlled, aligned to the reference genome (ARS-UCD 2.0) using HISAT2, and quantitative analysis was performed using HTseq to obtain the expression matrix of 32,626 annotated genes.

[0033] Step 2: Establish the feature set, including the following steps: Step 2.1: Perform genome-wide association analysis: Based on all GEBVs in the genetic assessment, calculate the effect of each SNP using the FarmCPU method, implemented in R version 4.3.2 using rMVP

[75] for GWAS. FarmCPU benefits from random effects models (REM) based on The optimal SNP selection formula is as follows: ; The results showed that nine significant SNPs were distributed in BTA 2, 6, 9, 11, 13 and 23. Genes located within 500 kb of significant SNPs were considered as potential candidate genes for the traits under study.

[0034] Step 2.2: Perform differential gene expression analysis: Based on the BHBA concentrations measured at 5, 7, and 12 days postpartum in 50 individuals, they were categorized into normal, subclinical, and clinical groups according to standards. Simultaneously, based on the principle of a single variable, two expression matrices were created for differential gene expression analysis: one for postpartum days and the other for pathological status. Specifically, the analysis was conducted at 7 days postpartum vs. 5 days postpartum, 12 days postpartum vs. 5 days postpartum, 12 days postpartum vs. 7 days postpartum, and between clinical vs. normal, clinical vs. subclinical, and subclinical vs. normal groups. Limma was used to perform differential expression analysis on the six control groups, using the following formula: ; in, For genes The log2 representation vector, For designing the matrix ( (including grouping, batching, etc.) For genes The vector of coefficients to be estimated, For the residual vector, It follows a normal distribution. This represents the residual variance.

[0035] The difference between the other party The compression formula is: ; in, For genes The sample variance For residual degrees of freedom, , The global prior mean. These are the a priori degrees of freedom.

[0036] Screening for gene expression changes ≥2 or ≤0.5 and The genes.

[0037] Step 2.3: Perform Cis-eQTL analysis: Using resequencing genomic data from 50 Holstein cattle, the MatrixEQTL model was set to modelLINEAR, with a distance threshold of 1 MB for cis-expressed quantitative trait loci (cisDist < 1 × 10⁻⁶). 6 Furthermore, the error covariance matrix (errorCovariance<-numeric()) was not used; the formula is as follows: ; in, For the sample Gene Standardized expression level, For the sample SNP The additive encoding (0 / 1 / 2). The effect size of cis-eQTL (to be estimated). For the sample On genes The residual.

[0038] The p-value threshold for significant cis-eQTL is 0.00001 (pvOutputThreshold_cis < 1 × 10⁻⁶). -5Only significant cis-eQTLs with FDR < 0.05 were retained, and the three selected cis-eQTLs were annotated.

[0039] Step 2.4: Perform co-location analysis: Based on the genome-wide association analysis (GWA) results of estimating breeding values ​​using BHBA in bovine blood, regions of 1 MB flanking nine significant SNPs were analyzed. This step used 612 SNPs. Since the located cis-eQTLs had little overlap with loci within this region, cis-eQTLs from bovine blood, mammary gland, and liver tissues were downloaded from FarmGTEx. The posterior probability 4 (PPH4) was applied in the COLOC algorithm to determine whether these two traits shared genetic variation. The relevant formula is as follows: ; ; in, The number of valid SNPs in the region. for In the assumption The Bayes Factor (obtained from GWAS) is below. Assumption The prior probability, Assumption The posterior probability, PPH4>0.8, indicates that GWAS and eQTL share genetic variations.

[0040] Step 2.5: Establishing a gene expression feature set: Integrating SNPs and QTLs located by genome-wide association analysis, differential gene expression analysis, cis-eQTL mapping, and co-location analysis, a total of 77 candidate regulatory genes were screened through gene annotation and (GO) and KEGG enrichment analysis. The different expression levels of these 77 genes at different days in 50 individuals were extracted to construct a gene expression feature set.

[0041] Step 3: Build a machine learning model, including the following steps: Step 3.1: Divide the established gene expression feature set into a training set and a validation set. The training set includes gene expression levels at 5 and 7 days postpartum, while the validation set includes gene expression levels at 12 days postpartum. Use BHBA ≤ 1.4 mmol / L as the Control and BHBA > 1.4 mmol / L as the Case. Use the R package glment to build a machine learning model, as shown in the following formula: ; ; ; in, The set of vectors representing the possible weights of each feature gene. For sample size, For sample number, For the first The true labels for each sample are: diseased = Case, healthy = Control. The model predicts the first The probability that a sample is a positive class. For L1 penalty intensity, The number of characteristic genes, For the first The regression coefficient of each gene, The intercept is... For the first Gene expression vectors of each sample, This is the transpose symbol.

[0042] Step 3.2: Use 5 rounds of 10-fold cross-validation and design a grid search (range: 0.1-0.01, 50 grids) to determine the optimal hyperparameters. The relevant methods are as follows: ; ; ; in, For the set of candidate regularization parameters, For the number of grid points, , Number the grid points. , The L1 penalty strength corresponds to the minimum / maximum of the mesh. , For linearly equidistant cut points, , For health status labels, For a given The average cross-validation performance, For the number of repetitions, , For folds, , For the first The repetition, the... The AUC value of the fold For optimal penalty intensity, For a given The average AUC obtained on cross-validation, To obtain To achieve the maximum intensity of punishment; The optimal Substituting into the final model (complete training set), using the optimal L1 penalty level, we obtain 16 genes with non-zero weights, and the output AUC=0.914. The AUC of these 16 genes in the validation set is 0.814.

[0043] Step 4: Evaluate the results. Use kernel density estimation (KDE) to visualize the grid search kernel. The AUC distribution of 5-fold 10-fold cross-validation is shown below. Figure 2 The calculation formula is: ; ; ; in, The current threshold, , For the first The model predicts the probability for each sample. For the first The true labels for each sample are: Control = negative, Case = positive. For indicator functions, when If the call is positive, select 1; otherwise, select 0. For the horizontal axis variable, Given an FPR value, return the corresponding threshold. , Threshold The actual rate, Integrate along FPR from 0 to 1; Using the R package ggplot for visualization is optimal. The ROC curves on the final model validation set are as follows: Figure 3 and Figure 4 As shown.

[0044] It is worth noting that all contents not described in detail in this invention are existing technologies and are well known to those skilled in the art.

[0045] Therefore, this invention proposes a prediction system for regulatory genes of ketosis in dairy cows based on multi-omics analysis and machine learning. It constructs an integrated technical system that combines multi-omics data integration, precise candidate gene screening, core gene identification, and synergistic ketosis risk prediction. Co-localization analysis is deeply integrated into the feature set construction process. By integrating genomic, transcriptomic, and metabolomic data, combined with genome-wide association analysis (GWAS), differential gene expression analysis, cis-eQTL localization, and co-localization analysis, it effectively distinguishes between causal variations and linkage disequilibrium variations, addressing the pain point of insufficient accuracy in single-gene screening methods. The Lasso model, through L1 regularization penalty, achieves synergy in invalid gene removal, core gene weight allocation, and prediction model construction. Combined with grid search and five 10-fold cross-validation iterations to optimize parameters, it overcomes the limitations of traditional prediction methods that are susceptible to environmental interference. Through dual evaluation using ROC curves and AUC values, the validation set AUC ≥ 0.8, ensuring prediction accuracy and model generalization ability. This invention significantly improves the efficiency of screening regulatory genes for ketosis in dairy cows and the accuracy of risk prediction, providing key technical support for the analysis of the molecular mechanism of ketosis pathogenesis, the development of early diagnostic biomarkers, molecular breeding, and precision prevention and control.

[0046] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A predictive system for regulatory genes of ketosis in dairy cows based on multi-omics analysis and machine learning, characterized in that, The system includes sequentially connected modules for data acquisition and preprocessing, feature set establishment, machine learning model establishment, and result evaluation. These modules work together to predict and validate regulatory genes for ketosis in dairy cows. The data acquisition and preprocessing module includes a phenotype collection unit, a genome data processing unit, and a transcriptome data processing unit. It collects phenotype data of bovine ketosis, outputs a ketosis-related dataset, processes multi-omics data, and outputs a basic data sub-library after quality control, comparison, and quantitative analysis. The feature set construction module includes a genome-wide association analysis unit, a gene differential expression analysis unit, a cis-eQTL analysis unit, and a colocalization analysis unit; The machine learning model building module includes a regulatory gene prediction unit, which constructs a Lasso model based on screened candidate regulatory genes, selects core regulatory genes and assigns weights through L1 regularization penalty, and optimizes model parameters by combining grid search and cross-validation. The results evaluation module uses ROC curves and AUC values ​​to assess model overfitting and the effectiveness of ketosis prediction.

2. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: The multi-omics data includes bovine genomic data, transcriptomic data, and metabolomic data. Among them, the genomic data contains SNP site information; the transcriptomic data contains gene expression level information; the metabolomic data contains BHBA content information in the blood; and the phenotypic data includes the disease status of the bovines, which is labeled with a binary label: diseased = Case, healthy = Control.

3. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 2, characterized in that: The classification criteria for BHBA concentration are as follows: BHBA concentration ≤ 1.4 mmol / L is considered healthy, 1.4 mmol / L < BHBA ≤ 3 mmol / L is considered subclinical ketosis, and BHBA > 3 mmol / L is considered clinical ketosis. Based on the ketosis classification results, subsequent phenotypic labeling and dataset construction are carried out.

4. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: The quality control criteria for genomic data were as follows: a pass rate of less than 95% for removing individual genotypes, a pass rate of less than 90% for removing individual SNP genotypes, a minor allele frequency (MAF) of <0.05 for SNPs, and a deviation from Hardy-Weinberg equilibrium (P < 1.0 × 10⁻⁶). -6 The mark.

5. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: The data acquisition and preprocessing module's transcriptome data processing workflow is as follows: after quality control of the sequencing data, it is aligned to the reference genome ARS-UCD2.0, quantitative analysis is performed, and the expression matrix of the annotated genes is obtained.

6. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: In the feature set construction module, the analysis results of the genome-wide association analysis unit, gene differential expression analysis unit, and cis-eQTL analysis unit are output as an intermediate data sub-library. The co-localization analysis results of the genome-wide association analysis and cis-eQTL results, combined with the gene differential expression analysis results, are output as a gene expression feature set, which stores the candidate genes selected.

7. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: The steps of localization analysis in the feature set establishment module are as follows: Step S1: Perform basic data organization, organize the basic information of the sample, including but not limited to pathological status and gene expression levels; Step S2: Perform genome-wide association analysis, using the ketosis-associated phenotype as the dependent variable and the pre-processed SNP loci as the independent variable, calculate the association strength P-value between each SNP locus and the ketosis phenotype, and screen for significantly associated SNP loci. Step S3: Using the ketosis classification phenotype as the dependent variable and gene expression level as the independent variable, screen for TOP genes with differential expression under different ketosis states; Step S4: Perform cis-eQTL analysis, with the difference in gene expression levels as the dependent variable and the pre-processed SNP sites as the independent variables. Set the distance threshold of cis-expression quantitative trait loci to 1 MB and screen for significant cis-eQTLs associated with gene expression changes. Step S5: Perform colocalization analysis, integrate GWAS results with matched eQTL information for association verification, calculate posterior probability by setting hypotheses, and quantify the possibility of SNP sites simultaneously regulating ketosis phenotype and gene expression. Step S6: Determine if there is co-localization of SNP sites, record the corresponding SNP sites as target SNP sites that simultaneously affect the ketosis phenotype and gene expression in dairy cows, and select the genes associated with the target SNP sites as candidate regulatory genes. Step S7: Establish a gene expression feature set, integrate GWAS, cis-eQTL and colocalization analysis results, and construct the gene expression feature set through gene annotation, metabolic pathway and functional enrichment analysis.

8. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: In machine learning models, the working steps of regulating gene prediction units are as follows: Step T1: Divide the feature set into a training set and a validation set using random sampling to ensure that the two sets of data have the same distribution. The training set is used for model training, feature selection and parameter optimization, and the validation set is used to evaluate model performance. Step T2: Construct a Lasso regression model based on the training set. The model formula is as follows: ； ； in, The set of vectors representing the possible weights of each feature gene. For sample size, For sample number, For the first The true labels for each sample are: diseased = Case, healthy = Control. The model predicts the first The probability that a sample is a positive class. For L1 penalty intensity, The number of characteristic genes, For the first The regression coefficient of each gene, The intercept is... For the first Gene expression vectors of each sample, It is the transpose symbol; Step T3: Feature selection of candidate regulatory genes is performed using L1 regularization penalty. The penalty formula is as follows: ； in, This is an L1 penalty item; Step T4: Automatically compress the weights of genes that do not contribute to the regulation of ketosis to 0 and remove them. Optimize the model regularization parameters using a grid search method and five 10-fold cross-validation iterations. The formula for the grid search method is: ； in, For the set of candidate regularization parameters, For the number of grid points, Number the grid points. The L1 penalty intensity is the value corresponding to the minimum mesh size. The L1 penalty intensity is the value corresponding to the maximum grid size. For linearly equidistant cut points, ; The formula for 5-fold 10-fold cross-validation is: ； in, For health status labels, For a given The average cross-validation performance, For the number of repetitions, For folds, For the first The repetition, the... The AUC value of the fold; Step T5: Select the optimal L1 penalty intensity, the selection rule is as follows: ; in, For optimal penalty intensity, For a given The average AUC obtained on cross-validation, To obtain To achieve the maximum intensity of punishment.

9. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: The results evaluation module systematically evaluates the output of the machine learning model building module, focusing on two core objectives: model overfitting and the effectiveness of ketosis prediction. ROC curves and AUC values ​​are used as the core evaluation metrics. Specific steps include: Step W1: Evaluate the model's effectiveness in overfitting and predicting ketosis. Plot ROC curves and calculate AUC values ​​based on the training and validation sets respectively. The ROC curve formula is: ； ； in, The current threshold, , For the first The model predicts the probability for each sample. For the first The true labels for each sample are: diseased = Case, healthy = Cnotrol. For indicator functions, when The value is 1 if the condition is met, otherwise it is 0. The formula for calculating AUC is: ； in, For the horizontal axis variable, Given an FPR value, return the corresponding threshold. , Threshold The actual rate, Integrate along FPR from 0 to 1; Step W2: Compare the AUC values ​​of the training set and the validation set to determine whether overfitting has occurred, and evaluate the effectiveness of ketosis prediction based on the AUC values ​​of the validation set.

10. The bovine ketosis regulatory gene prediction system based on multi-omics analysis and machine learning according to claim 1, characterized in that: The evaluation criteria for the results evaluation module are as follows: Overfitting assessment: A difference of ±0.1 between the AUC values ​​of the training set and the validation set indicates that the model has strong generalization ability; Predictive effectiveness assessment: Comparing the ROC curves and AUC values ​​of the training and validation sets, a validation set AUC value ≥ 0.8 indicates that the core regulatory genes can effectively predict ketosis.

Images