CN115995262B - Method for analyzing corn genetic mechanism based on random forest and LASSO regression - Google Patents

Abstract

The invention belongs to the field of system biologyThe field of genetics, in particular to a method for analyzing a corn genetic mechanism based on random forests and LASSO regression. The invention carries out linkage disequilibrium analysis on genotype data after corn pretreatment and deletes r ² After SNP loci with a value greater than 0.2, the remaining loci form a new feature subset; determining important influencing factors, wherein the influencing factors form a new feature subset; establishing a LASSO analysis model through the feature subset and the lipidomic data, and predicting each lipid metabolite phenotype; establishing a LASSO analysis model using predicted lipid metabolite data values in combination with phenotypic data, predicting the phenotypic data using a determinant coefficient R ² The model is evaluated. The invention analyzes the genetic mechanism of the corn lipid group through screening representative feature subsets and carrying out multi-group data combined analysis, and provides a theoretical basis for the genetic improvement of the corn high-yield inbred line.

Description

A method for analyzing the genetic mechanism of maize based on random forest and LASSO regression

Technical Field

The present invention belongs to the field of systems biology\genetics, and specifically relates to a method for analyzing the genetic mechanism of corn based on random forest and LASSO regression.

Background Art

Genome-wide association study (GWAS) is an effective method to analyze the genetic variation of complex quantitative traits. Thousands of related loci have been successfully identified by establishing direct associations between genotypes and human diseases and plant quantitative traits. In nature, plants are affected by developmental stages, tissues and environmental stimuli to jointly produce a large number of metabolites, ranging from about 100,000 to 1 million, so the metabolome is often used as a "bridge" between genotypes and plant phenotypes. Metabolic phenotypes are quantitative traits, and phenotypic data vary greatly. Metabolome data are used as intermediate phenotypes. Association analysis with genotype data can locate more genes and may also locate rare SNP sites. However, the variants identified from GWAS can only explain part of the total genetic variation, and the predictive ability of phenotypes is limited.

Machine learning and deep learning methods have shown great potential in processing biological data. In biology, machine learning algorithms are mainly used in two aspects: one is to make accurate predictions where experimental data is lacking, and use these predictions to guide future research work; the other is to use machine learning algorithms to explain the phenotype of potential genotypes and gain new insights from them, thereby deepening our understanding of biology. Single nucleotide polymorphism (SNP) refers to the polymorphism of DNA sequence caused by changes such as insertion or deletion of a single nucleotide (A, T, C, G) at the whole genome level, and is the most widespread genetic variation in the population. Each SNP can explain a small part of the phenotypic variation, but a large number of SNPs together may account for a larger proportion. When using high-dimensional genomic data, the risk of overfitting is high. Feature selection techniques are used to identify SNPs associated with complex traits, allowing researchers to focus their research on the most influential SNPs.

Summary of the invention

In view of the problem that traditional association analysis methods can only explain part of the heritability and cannot identify micro-effect SNPs, the present invention provides a multi-omics analysis method based on random forest and LASSO regression to analyze the genetic mechanism of corn lipidome. The method evaluates the importance of features through random forest filtering method, screens influential feature subsets, and further establishes a two-layer LASSO regression model using the feature subsets to predict the phenotype. The feature subsets obtained by the random forest algorithm and the feature subsets and phenotype group data obtained by the first LASSO regression analysis are analyzed to analyze the genetic mechanism of corn lipidome, thereby providing ideas and methods for breeding high-yield and tolerant corn.

In order to achieve the above-mentioned object of the invention, the specific technical solution adopted by the present invention is:

A random forest and LASSO regression model is based on the following steps:

(1) Preprocessing of corn multi-omics (genomics, lipidome, phenotype) data, such as removing outliers, filling missing values, and data standardization;

(2) Linkage disequilibrium analysis was performed on the preprocessed genotype data. The degree of linkage disequilibrium between SNPs was calculated using the plink technology and represented by the standardized index r2 ^. After deleting SNP sites with ^r2 (linkage disequilibrium parameter) values greater than 0.2, the remaining sites constituted a new feature subset.

(3) A random forest model is established through the new feature subset to rank the importance of variables and determine the important influencing factors. These influencing factors constitute a new feature subset.

(4) Establishing the LASSO analysis model based on the above feature subsets to predict the phenotype of each lipid metabolite;

(5) The predicted lipid metabolite data values were used to establish the LASSO analysis model, predict the phenotypic data, and evaluate the model using the coefficient of determination ^R2 .

Furthermore, the method for obtaining the genotype data is as follows: the genotype file is a file in hapmap format, each row represents a material, each column represents a SNP site, and the original value is the genotype of the SNP. Before subsequent analysis, the data needs to be re-encoded to convert the genotype into a numerical type using a 0-1-2 encoding method.

Specifically, the present invention provides a multi-omics analysis method based on random forest and LASSO analysis, and the specific operation steps are:

(1) The genotype file is a hapmap format file. Each row represents a material, each column represents a SNP site, and the original value is the SNP genotype. Before subsequent analysis, the data needs to be recoded to convert the genotype into a numerical type using a 0-1-2 encoding method.

(2) Check whether there are missing values in the genotype data, lipid group data, and phenotype group data. If so, use the KNN algorithm to interpolate the genotype data and use the mean to fill in the lipid group data and phenotype group data;

(3) Considering the linkage disequilibrium (LD) problem between SNPs, the features with linkage disequilibrium parameter ^r2 greater than 0.2 were deleted using plink technology, and the remaining features were used to construct feature subsets;

(4) Input the feature subset from the previous step into the random forest model, use the 10-fold cross-validation method to evaluate the accuracy, output the results, and save them as a new subset;

(5) The feature subset obtained in the previous step was used to construct a LASSO regression model, and the parameters were estimated using the 5-fold cross-validation method to predict each lipid metabolite;

(6) Use the predicted values of lipid metabolites in the previous step to predict the phenotype and use the coefficient of determination ^r2 to evaluate the model; integrate the results of steps 4 and 6 to analyze the genetic mechanism of the corn lipidome.

Furthermore, the steps of evaluating variables in the random forest model in step (4) of the above specific operation are:

Step 1: Use the bootstrap sampling method to select n samples from the sample as the training set;

Step 2: Generate a decision tree and each node using the sample set; repeat the sampling of samples and variables with replacement, repeat k times, and generate k decision trees;

Step 3: Use the obtained random forest to predict the test sample, and use the voting method to determine the prediction result. The category where the majority of k decision trees belongs is the prediction result.

Step 4: Use random forest and 10-fold cross validation to screen important features.

At present, traditional association analysis methods have discovered thousands of SNP sites associated with phenotypes. For most complex quantitative traits, these single SNP sites have little ability to cause phenotypic variation, and only a small number of people can use these sites to explain phenotypic variation. Researchers believe that quantitative phenotypic variation is caused by the interaction between genes and genes, and genes and the environment, while association analysis methods can only detect the correlation between a single SNP site and the phenotype, and lack the ability to detect the interaction of multiple genes.

The beneficial effects of the present invention compared with the prior art are: the present invention can handle the complex interactions of multiple genes and environmental factors, and does not need to assume that the data conforms to a specific distribution or model; it can also discover new associations and laws, improve prediction accuracy and interpretability; and can automatically perform operations such as feature selection and dimensionality reduction, reducing tedious manual operations. In addition, compared with traditional association analysis methods, the SNP sites identified by the present invention can explain more phenotypic variations.

BRIEF DESCRIPTION OF THE DRAWINGS

Figure 1 is a flow chart for analyzing the genetic mechanism of corn based on random forest and LASSO regression models.

DETAILED DESCRIPTION

In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

Example 1 The genetic materials used in the data of this article include a maize association population containing 513 inbred lines (YangX, Gao S, Xu S, et al. Characterization of a global germplasm collection and its potential utilization for analysis of complex quantitative traits inmaize [J]. Molecular Breeding, 2011, 28: 511-526.) and a RIL (B73/By804) population containing 197 inbred lines, which are derived from the normal line B73 and the high-oil line By804 (Xiao Y, Tong H, Yang X, et al. Genome-wide dissection of the maize ear genetic architecture using multiple populations [J]. New Phytologist, 2016, 210 (3): 1095-1106.). The materials were planted in the field test area of Huazhong Agricultural University, Wuhan.

In this embodiment, a method for predicting plant phenotype based on genome and lipidome is provided. FIG1 is a flow chart of analyzing the genetic mechanism of corn based on random forest and LASSO regression model, comprising the following steps:

(1) Data preprocessing

The genotype file is a 541*1145174 non-numerical matrix (as shown in Table 1), containing missing values, with rows being SNPs and columns being samples. Each value in the matrix represents the allele of the SNP in the material, which has four forms: A, T, C, and G. The lipidome file is a 479*154 numerical matrix (as shown in Table 2), containing missing values, with rows being materials and columns being lipid metabolites. The phenotypic data is a 541*6 numerical matrix (as shown in Table 3), containing missing values, with rows being materials and columns being traits.

Table 1 Genotype raw data

。

.

Table 2 Lipidome raw data

。

.

Table 3 Raw phenotypic data

。

.

First, the genotype file was recoded and converted into a numerical file for subsequent calculations. The hmp format file was converted into a vcf format file using the tassel software on the Linux system. The converted vcf file was converted into the plink (.ped, .map) format using the vcftools software. The plink tool was used to calculate the frequency of alleles at each site in the material. The allele with a high frequency was regarded as the major allele and coded as 1. The other allele was regarded as the minor allele and coded as 2. Missing values were represented by NA (as shown in Table 4).

Table 4 Recoded genotype data

。

.

Non-normally distributed data often contain many outliers, which may have a significant impact on the modeling results. Therefore, before subsequent calculations, the lipidome data and phenotypic data were subjected to Boxcox processing (as shown in Tables 5 and 6) to make the data conform to or be close to normal distribution.

Table 5 Data obtained after normalization of lipid group raw data

。

.

Table 6 Data obtained after standardization of original phenotypic data

。

.

The genotype data contains missing values, which are filled using the KNN algorithm. The similarity between each sample and other samples is calculated using the Euclidean distance. For each missing value, the five nearest neighbor samples are found, and the average of the known genotypes of the five samples is used as the estimated value of the missing value. This step is repeated until all missing values are filled (as shown in Table 7). The calculation formula for the Euclidean distance is:

。

.

d(x,y) represents the Euclidean distance between a point x( _x1 , _x2 , ..., _xn ) and a point y( _y1 , _y2 , ..., _yn ) in n-dimensional space, where _xi represents the i-th coordinate of the first point; _yi represents the i-th coordinate of the second point; and n is the vector dimension.

The mean of the standardized lipidome data and phenotypic data was used to fill in missing values (as shown in Tables 8 and 9), and the mean of each lipid or trait in the population was calculated and used as the estimated value of the missing value.

Table 7 Genotype data after KNN filling missing values

。

.

Table 8 Lipidome data after mean imputation

。

.

Table 9 Phenotypic data after using mean to fill missing values

。

.

(2) Filter redundant features based on linkage disequilibrium between SNPs

Linkage disequilibrium (LD) refers to the non-random association between two or more loci, which are located on the same chromosome and are often inherited to offspring at the same time. This association may lead to a high correlation between characteristics, thus affecting the accuracy and explanatory power of the model. According to the genotype data matrix, the linkage disequilibrium coefficient (r ² ) between each pair of SNPs is calculated; for each SNP, one or more adjacent SNPs with the highest LD degree are selected and let them represent the SNP; the LD value between each representative SNP and other representative SNPs is calculated. If the LD value of a representative SNP with any other representative SNPs exceeds the preset threshold, the representative SNP is deleted. The threshold is selected as 0.2, and it is considered that there is linkage disequilibrium between SNPs with r ² greater than 0.2. This method can effectively reduce redundant information and noise in the data set and improve the efficiency and accuracy of genetic association studies. A total of 88,766 SNP loci were obtained (as shown in Table 10). The following calculation takes the 100_mM_NaCL_Na ⁺ phenotype as an example.

Table 10 Variables obtained after removing LD

。

.

(3) Build a random forest model and evaluate feature importance

Use the feature file obtained after the previous step to build a random forest model, evaluate the feature importance, and further extract the features. The genotype file (a matrix containing 88,766 SNP sites) obtained from step (2) is used as the input data set for this step; a random forest model is built on the data set, which contains multiple decision trees, and different training sets and test sets are generated using bagging and random sampling techniques; for each decision tree, the best split variable is selected at the node, and the mean decrease accuracy of the variable at the node is calculated as its importance score. These scores can be weighted averaged by various methods to obtain the importance ranking of each variable in the entire model; the variables are sorted according to these rankings, and the t variables with the highest scores are selected as the final selection results to form the feature subset. t is determined by the cross-validation method, and the cross-validation adopts the 10-fold cross-validation method. After the above steps are calculated, 13,456 variables with high scores are finally obtained (as shown in Table 11).

Table 11. Important variables obtained using the random forest method

。

.

(4) Establishing a regression model for lipid metabolite data

Using the feature subset obtained in step (3), the LASSO regression model is constructed to calculate the predicted value of each lipid metabolite. The feature subset obtained in step (3) is used as the input data set of this model; then, the cv.glmnet function in the glmnet package is used to construct the LASSO regression model. This function uses the cross-validation method to select the optimal regularization parameter λ and obtain the corresponding coefficient estimates; then, the variables are sorted according to the coefficient estimates and the important feature subsets are screened out; the response variable is predicted using the feature subset, and the prediction results are output and saved for the next calculation; the above three steps of model construction, subset screening, and response variable prediction are repeated until each lipid metabolite in the lipidome data is used as a response variable to build a model and obtain the prediction results (as shown in Table 12).

Table 12 Predicted values of lipid metabolites

。

.

(5) Establishing regression model based on phenotypic data

The prediction results of each lipid metabolite obtained in step (4) were used to form a feature subset, and a LASSO regression model was constructed to calculate the predicted value of the phenotype.

The feature subset obtained in the fourth step was used as the input data set of this model; then, the cv.glmnet function in the glmnet package was used to construct the LASSO regression model. This function uses a 5-fold cross-validation method to select the optimal regularization parameter λ and obtain the corresponding coefficient estimates; then, the variables are sorted according to the coefficient estimates and the important feature subsets are screened out (as shown in Table 13); the model is evaluated using the test data, and relevant indicators (such as root mean square error RMSE, determination coefficient ^R2, etc.) are calculated to measure the model performance. In this section, we calculated that the RMSE value is 0.8001393 and the determination coefficient ^R2 value is 0.2808084, indicating that the combination of genotype data and lipidome data can explain 26% of the phenotypic variation.

Table 13 Screened important lipid metabolites

。

.

(6) In step (5), we obtained lipid metabolites that were significantly correlated with the phenotype. Based on these lipid metabolites, we can find the corresponding significant SNP sites in step (4). Based on these SNP sites, we can use annovar and other methods to annotate and mine candidate genes (as shown in Table 14). We can further analyze the genetic mechanism of the maize lipidome. For example, the gene Zm00001d027395 (NAC78) of the maize NAC transcription factor family may respond to salt stress by affecting the content of DGDG36_4 metabolites.

Table 14 Candidate genes annotated according to significant sites

。

.

Claims (5)

1. The method for analyzing the corn genetic mechanism based on random forest and LASSO regression is characterized by comprising the following steps of:

(1) Preprocessing a plurality of groups of chemical data of corn;

(2) Carrying out linkage disequilibrium analysis on the pretreated genotype data, deleting SNP loci with linkage disequilibrium parameter values greater than 0.2, and forming a new feature subset by the residual loci;

(3) Establishing a random forest model through the new feature subset and the phenotype data, sequencing the importance of the variables, and determining important influence factors, wherein the influence factors form the new feature subset;

(4) Establishing a LASSO analysis model through the feature subset and the lipidomic data, and predicting each lipid metabolite phenotype;

(5) Establishing a LASSO analysis model using predicted lipid metabolite data values in combination with phenotypic data, predicting the phenotypic data using a determinant coefficient R ² And (5) evaluating the model and analyzing the genetic mechanism of the corn.

2. The method of claim 1, wherein the plurality of sets of chemical data of step (1) are genomic data, lipidomic data, and phenotypic data, respectively.

3. The method of claim 1, wherein the genotype data obtained in step (2) is obtained by: the genotype file is a hapmap format file, each row represents a material, each column represents a SNP locus, the original value is the genotype of the SNP, the data needs to be recoded before the subsequent analysis is carried out, the genotype is converted into a numerical type, and a coding mode of 0-1-2 is adopted.

4. The method according to claim 1, characterized by the specific operating steps of:

(1) The genotype file is a hapmap format file, each row represents a material, each column represents a SNP locus, the original value is the genotype of the SNP, the data is required to be recoded before the subsequent analysis is carried out, the genotype is converted into a numerical type, and a coding mode of 0-1-2 is adopted;

(2) Checking whether the genotype data, the lipidosome data and the phenotype group data have missing values, if so, interpolating the genotype data by using a KNN algorithm, and filling the lipidosome data and the phenotype group data by using a mean value;

(3) Taking into account the problem of linkage disequilibrium between SNPs, the linkage disequilibrium parameter r was deleted using the plink technique ² Features greater than 0.2, constructing a feature subset from the remaining features;

(4) Inputting the feature subset of the previous step into a random forest model, adopting a 10-fold cross validation method to evaluate accuracy, outputting a result, and storing the result into a new subset;

(5) The feature subset obtained in the last step is used for constructing a LASSO regression model, a 5-time cross validation method is adopted for estimating parameters, and each lipid metabolite is predicted;

(6) Predicting the phenotype using the predicted value of the lipid metabolite of the previous step, and using the determination coefficient r ² Evaluating the model; and (3) integrating the results of the step (4) and the step (5), and analyzing the genetic mechanism of the corn lipidome.

5. The method of claim 4, wherein the step of evaluating the variables of the random forest model in step (4) is:

the first step: selecting n samples from the samples by using a sampling method with a put-back function as a training set;

and a second step of: generating a decision tree and each node by using the sample set; repeatedly sampling the sample and the variable with the replaced sample and the variable, and repeating the sampling for k times to generate k decision trees;

and a third step of: predicting a test sample by using the obtained random forest, determining a prediction result by using a voting method, and if the category of the mode in the k decision trees is the prediction result;

fourth step: important features were screened using random forests and 10 fold cross validation.

Images