Single-cell pathway analysis method developed based on Overlapping Group Lasso

Abstract

This invention proposes a single-cell pathway analysis method based on Overlapping Group Lasso. The model's logistic form is used to regress gene expression data and cell state for each cell. A Stability Selection algorithm is employed, performing multiple sampling regressions on the data and retaining the pathways with the most stable screening frequency as the final result. The specific steps are as follows: Step 1: Mathematical expression and preprocessing of the input cell gene expression matrix, pathway data, and cell labels; Step 2: Establishing a regression model; Step 3: Solving the regression model; Step 4: Stability selection. The model selects independent variable groups at the pathway level and, based on the characteristics of Overlapping Group Lasso, decomposes the effects of overlapping genes in different genomes, effectively avoiding the influence of overlapping gene effects.

Description

Technical Field

[0001] This invention proposes a single-cell pathway analysis method based on Overlapping Group Lasso, which can be directly applied to the field of pathway analysis related to single-cell sequencing technology. Background Technology

[0002] Single-cell RNA-sequencing, based on the central dogma and existing computer algorithms, can translate a complex mRNA sequence into a cell-gene expression matrix. Each row represents a cell, each column represents a gene, and the elements in the matrix represent the expression level of a gene in a particular cell. This technology was initially applied to situations where biological materials are scarce (such as the embryonic development stage).

[0003] In recent years, this technology has developed rapidly, enabling more precise characterization of gene expression in thousands of cells at a given moment. With the increase in data volume, many new research questions have emerged, such as integrating data from multiple single cells, identifying the type of each cell, and eliminating data noise.

[0004] One of the research questions is pathway analysis—based on gene expression data from two or more different types of cells, identifying genomes (pathways) with significant differences in gene expression between two cell populations and biological significance. The most commonly used method for pathway analysis is GSEA[1]. This method first ranks the expression level of each gene relative to cell type, then checks the gene ranking of each pathway for statistical testing and assigns a score to the pathway. Finally, after regularizing all pathway scores, it obtains the pathways with statistically significant differences in gene expression between the two datasets. Since pathway analysis can transform complex high-dimensional gene expression data into interpretable results, it has become an important part of the single-cell data analysis process.

[0005] The background technology relies on the Overlapping Group Lasso algorithm [2] and the Stability Selection algorithm [3] in sparse learning methods. There is a core idea in sparse learning: although the amount of data is getting bigger and bigger, the underlying generation mechanism may actually be controlled by only a few simple parameters. Based on this idea, when faced with a huge dataset, we often need to find the most important variables.

[0006] Before introducing the Overlapping Group Lasso algorithm, it is necessary to introduce the Lasso algorithm and the Group Lasso algorithm [4]. The Lasso algorithm mainly explores which independent variables are most important to the dependent variable among multiple independent variables. It adds an L1 regularization term containing the parameters to be found after the loss function of linear regression and uses the geometric properties of the L1 regularization term to achieve the purpose of screening a single variable. When it is known that a single variable is in different variable groups, the Group Lasso algorithm adds a regularization term to the parameters of each group of independent variables and performs regression calculation on the variable group as a unit. It is often used to screen the independent variable group related to the dependent variable. The Overlapping Group Lasso algorithm improves on this basis so that it can also screen the independent variable group that is important to the dependent variable when there are the same variables between the variable groups. In this invention, the pathway (genome) is a variable group. There are some common independent variables between each pathway. One of the important steps of this invention is to use the Overlapping Group Lasso algorithm to screen the genomes with significant differences in gene expression in two cell populations.

[0007] After selecting variables or variable groups, algorithm errors or inaccurate parameter selection often result in meaningless or inaccurate selections. Therefore, this invention employs a Stability Selection algorithm to ensure the accuracy of the final selected variables or variable groups with a high probability. Specifically, this involves multiple samplings and probability estimations of all sampled variables or variable groups, followed by setting appropriate thresholds to identify more stable variables or variable groups.

[0008] The following problems exist in the existing technology:

[0009] 1. Since the Overlapping Group Lasso is only a linear model, it only captures the linear relationship between gene expression levels and cell state, but the non-linear relationship between the two is not taken into account in this model.

[0010] 2. When dealing with large datasets, the Overlapping Group Lasso computation is slow due to the slow convergence speed of the proximal gradient algorithm and the excessive number of experiments set in the Stability Selection algorithm.

[0011] 3. Since it is unknown which genes and pathways actually affect the dependent variable in real single-cell datasets, the pathway screening results obtained are only of a certain reference value.

[0012] Most existing models are improvements on previous models, but the amount of gene expression data and pathway data in those earlier models was insufficient, leading them to consider this effect unimportant or negligible. However, with the increase in data volume, this effect can no longer be ignored.

[0013] References:

[0014] [1]Meinshausen N,Bühlmann P.Stability selection[J].Journal of the Royal Statistical Society Series B: Statistical Methodology, 2010,72(4):417-473.

[0015] [2]Jacob L,Obozinski G,Vert J P.Group lasso with overlap and graphlasso[C] / / Proceedings of the 26th annual international conference onmachinelearning.2009:433-440.

[0016] [3]Meinshausen N,Bühlmann P.Stability selection[J].Journal of the Royal Statistical Society Series B: Statistical Methodology, 2010,72(4):417-473.

[0017] [4] Yuan M, Lin Y. Model selection and estimation in regression with grouped variables [J]. Journal of the Royal Statistical Society Series B: StatisticalMethodology, 2006, 68 (1): 49-67. Summary of the Invention

[0018] Existing pathway analysis methods (such as GSEA) mostly examine pathways on a single genome basis, ignoring the impact of shared genes across different genomes. Other methods following GSEA also have similar issues. Currently, few pathway analysis tools address the impact of overlapping genes.

[0019] This invention takes into account that the same gene has different expression mechanisms in different pathways, and aims to decompose the effects of these genes to re-screen for important pathways and avoid screening pathways that are incorrect.

[0020] This invention uses the logistic form of the Overlapping Group Lasso model and the StabilitySelection algorithm to screen gene expression data from two cell types and predict cell states accordingly. Figure 1 The inputs and outputs of the invention are simply illustrated in the text.

[0021] First, compared to previous methods that analyzed individual genes (such as Lasso, GSEA, etc.), we used the logistic regression of the Overlapping Group Lasso model on a genome-wide basis to regress gene expression data and cell state for each cell. Because previous methods ignored the overlap effect caused by identical genes across genomes, the Overlapping Group Lasso model addresses this issue by introducing latent variables for overlapping genes.

[0022] Furthermore, since Lasso and subsequent Overlapping Group Lasso models both face the problem of sensitivity to penalty term parameters, this invention employs the Stability Selection algorithm, which performs multiple sampling regressions on the data and retains the path with the most stable selection frequency as the final result.

[0023] This invention discloses a single-cell pathway analysis method based on Overlapping Group Lasso, comprising the following steps:

[0024] Step 1: Perform mathematical expression and preprocessing on the input cell gene expression matrix, pathway data and cell labels.

[0025] set up Where m is the number of cells, n is the number of genes, and the row vector x i Let x′ be the gene expression vector of the i-th cell. i This represents its transpose, and the elements in the matrix are... Let represent the expression level of the j-th gene in the i-th cell. The matrix data is log-normalized before being input into the model. Let there be a set of K pathways. The k-th path is then represented as a set. All pathways are uniquely labeled as {0, ..., K-1}, |G k | represents the number of genes in the k-th pathway; this model only considers |G.k | ≥15 (i.e., pathways with more than 15 genes) cases. Gene pathway vector Where n is the number of genes, p j ∈{0, 1, ..., K-1} represents the pathway number to which the j-th gene belongs. This vector represents the pathway numbers to which all genes belong. Cell state vector m is the number of cells, y i ∈{0,1} is a marker indicating whether the i-th cell is diseased (if y i =0, then the cell is a normal mouse cell; if y = 0, then the cell is a normal mouse cell; i If the value is 1, then the cell is a diabetic mouse cell.

[0026] Step 2: Build a regression model. Embed pathway information into the preprocessed gene expression data, and then model these data using the logistic form of the Overlapping Group Lasso model;

[0027] Based on the assumptions made in the previous step, the input matrix is The output is Let the regression model be as follows:

[0028]

[0029] in, This represents the expression matrix of all genes in the k-th pathway (i.e., the data of all genes in the k-th pathway in the corresponding columns of X). This represents the regression coefficient of each gene in the k-th pathway, where c is a constant term. Let be Gaussian noise, representing variables other than gene expression that affect cell state. Since the pathway to which each gene belongs is known, let To replicate the sample matrix of overlapping genes in each pathway, the optimization problem can be expressed as:

[0030]

[0031] in, For loss function, As a penalty term, this term aims to screen out pathways that affect cell state, ensuring that the regression coefficients of genes in the affected pathways are not zero; while the regression coefficients of genes in the unaffected pathways are all zero. This approach decomposes the effect of each gene in different pathways, without fixing the effect of the gene in all pathways.

[0032] Let v be the latent vector representing the regression coefficient of the gene in the k-th pathway. If the j-th gene is in this pathway, then during optimization, this gene will be in v. k The value is 0 if it is in the middle, and not 0 otherwise.

[0033]

[0034] This model decomposes the coefficient β of all genes into the sum of K latent vectors (as shown in Equation 3). The core idea is to decompose the influence of overlapping genes on cell state into each pathway containing that gene to avoid the following situation: assuming pathway A and pathway B have overlapping genes g:

[0035] 1. In one pathway (e.g., A), all coefficients for genes are 0, while in another pathway (e.g., B), all coefficients except for gene g are non-zero. This indicates that the effect of pathway A in the overlapping genes covers the effect of pathway B.

[0036] 2. The pathway gene coefficients for both A and B are not zero, but the coefficient for gene g is the same. This indicates that the overlapping gene effects of pathways A and B are entangled and have not been resolved.

[0037] Because y i Since ∈{0,1}, we need to perform Logistic regression to transform the linear output into a probability model.

[0038]

[0039] According to Equation 2, this problem is a minimization problem, so the loss function should be the negative log-maximum likelihood estimate:

[0040]

[0041] From equations 2 and 5, let The optimization problem can be expressed as:

[0042]

[0043] Only need to solve for V, and then through The total regression coefficient can be obtained, and more importantly, each column vector in V can represent the effect of each pathway on the cellular diabetes state.

[0044] Step 3: Solve the regression model. Based on the model in Step 2, select an appropriate parameter range and use the proximal gradient descent algorithm to solve the regression model to obtain the gene expression coefficients for each pathway;

[0045] Since this optimization problem contains non-smooth terms, this invention employs the proximal gradient descent algorithm to handle the aforementioned optimization problem (an algorithm for handling optimization problems with non-smooth terms). In the t-th iteration, the smooth terms are... Gradient descent on β:

[0046]

[0047] in, The prediction results are used before each iteration. Then, the minimum value problem of the non-smooth term is solved by using the proximal operator. The proximal operator is calculated once for each of the K genomes to constrain them.

[0048]

[0049] Based on the above calculations, the proximal gradient algorithm applied to this optimization problem is shown in Table 1:

[0050] Table 1: Proximal gradient descent algorithm combined with Overlapping Group Lasso

[0051]

[0052]

[0053] Step 4: Stability Selection. A stability selection algorithm is added to the algorithm in Step 3. Data is repeatedly sampled and optimal parameters are tested to evaluate the importance of each pathway. Finally, the gene expression coefficients of important pathways are obtained, and irrelevant pathways are screened out.

[0054] Stability selection is an algorithm used to screen variables in sparse models, given p genomes and l distinct parameter values ​​for an Overlapping Group Lasso model within a given range. For each parameter value λ i By repeatedly sampling data M times to fit the model, we can obtain the number of times, N, the j-th independent variable is selected in the M fitted models under a certain parameter. Based on N and M, we can calculate the selection frequency of the j-th independent variable. Usually according to f j (λ i A line graph can be plotted to visualize the screening process, showing the relationship between each parameter value and the screening frequency for each genome. Based on this graph, a frequency screening threshold th∈[0,1] can be set to screen variables. Then the j-th genome can be considered as a stable independent variable that has been selected.

[0055] The beneficial effects of this invention compared to the prior art are as follows:

[0056] 1. The model screens out independent variable groups based on pathways, and based on the characteristics of Overlapping Group Lasso, it decomposes the effects of overlapping genes in different genomes, effectively avoiding the impact of overlapping gene effects.

[0057] 2. The Overlapping Group Lasso model combined with the stability selection algorithm can make the screening pathways accurate and prevent the subjective adjustment of the model hyperparameter λ from affecting the final screening results. This combination makes the screening pathway process interpretable.

[0058] 3. In addition to the pathway filtering function of the Overlapping Group Lasso model, this model also has good predictive performance, achieving a prediction accuracy of over 90% in both simulated and real datasets. Attached Figure Description

[0059] Figure 1 This is a simplified flowchart of the present invention. The input consists of two different types of single-cell data, which, after passing through the Overlapping Group Lasso and Stability Selection algorithms, yields effective pathways.

[0060] Figure 2 Comparison of regression coefficients between two datasets.

[0061] Figure 3 For T1D pathway q under different λ values k value.

[0062] Figure 4 For T2D pathway q under different λ values k value.

[0063] Figure 5 shows the results of variable selection for the two datasets using the stability regression algorithm. Detailed Implementation

[0064] This invention compares the results after applying it to two mouse single-cell datasets. The data used were α-single-cell data from mice aged 16-18 weeks, and were divided into two parts:

[0065] 1. Single-cell alpha data of T1D (type 1 diabetes) / normal mice: a total of 4756 cells, of which 3383 were from T1D mice and 1373 were from normal mice.

[0066] 2. Single-cell alpha data of T2D (type 2 diabetes) / normal mice: a total of 6587 cells, including 5214 cells from T2D mice and 1373 cells from normal mice.

[0067] Then, pathways with fewer than 15 genes and those lacking any of the aforementioned 1278 genes were removed from the MsigDB mouse M2 pathway database, resulting in a final count of 152 pathways. Next, the expression data of overlapping genes appearing in different pathways in all cells (one column in the data) were copied by b.j Next, b j Let be the number of times the j-th gene appears in all pathways. Then, each column of the gene expression matrix is ​​recombined according to the pathway order, resulting in a total of 5233 columns. The final model training data matrix is ​​denoted as .

[0068] get Then, it is logarithmically standardized. Let... The mean vector of each column is The standard deviation vector of each column is Then for each cell (i.e.) Gene expression data (each row) Perform the following conversion:

[0069]

[0070] After preprocessing the data, the learning rate α was set to 0.1, and the parameter λ ranged between [0.01, 0.02, ... 0.1]. With different fixed parameters λ, for T1D or T2D data, gene expression data of 1373 α cells (the same number as all normal mouse α cells) were randomly extracted from the α cells of diseased mice. These were then combined with the gene expression data of normal mouse α cells to obtain two training datasets: T1D and T2D, each containing 2746 α cells. The ratio of normal mouse α cells to diseased mouse α cells in each dataset was 1:1. Model fitting was then performed on both training datasets, and the estimated parameters were... The parameter estimation for each pathway is based on v k Arranging them sequentially yields the coefficients for each gene, such as... Figure 2 As shown.

[0071] Figure 2 The horizontal axis represents all genes in The corresponding positions (including overlapping genes after duplication) have a number equal to The column dimension is 5233; the vertical axis represents the parameter estimate for each gene. The top and bottom figures represent the regression results of T1D and normal mouse α cells, respectively. The red vertical lines separate different pathways, and the different colors of the β values ​​represent the regression results of fixing two parts of the data under different λ values. The larger the λ value, the fewer gene pathways are screened.

[0072] The important pathways identified after regression analysis of the two datasets differed under different λ values. For T1D, the regression results were relatively stable under different λ values; however, for T2D, the regression results fluctuated under different λ values. The pathways identified under different λ values ​​and the coefficients of each gene within these pathways were extracted from the two datasets. The l1 norm of the coefficients of all genes within each pathway was calculated and averaged, denoted as:

[0073]

[0074] q k The magnitude of λ can reflect the influence of each pathway on the diabetic state. For regression parameters of T1D and normal mouse data, the non-zero q values ​​at each λ value are used. k By comparing the values, we can obtain the following: Figure 3 As shown.

[0075] Figure 3 The model was shown to screen pathways corresponding to q in T1D and normal mouse α-cell datasets. k Value. The horizontal axis in the graph represents the name of each pathway, and the vertical axis represents q. k The values ​​are represented by different colors, where each color represents the result at a different λ value. Among the five pathways, the q values ​​of REACTOME_IMMUNE_SYSTEM and WP_MECHANISMS_ASSOCIATED_PLURIPOTENCY are much higher than those of the other pathways. These two pathways represent the immune system pathway and the cell differentiation potential pathway of α cells, respectively. These results are consistent with the existing medical literature's understanding of mouse T1D.

[0076] For the regression parameters of T2D and normal mouse data, the non-zero q values ​​at each λ value are... k By comparing the values, we can obtain the following: Figure 4 As shown.

[0077] Figure 4 This shows the q corresponding to the pathway selected by the model. k The values ​​are represented by different colors, indicating the results at different λ values. These pathways have higher q values ​​compared to other pathways.

[0078] (1)REACTOME_DEGRADATION_OF_THE_EXTRACELLULAR_

[0079] MATRIX (Extracellular Degeneration-Related Pathway)

[0080] (2) REACTOME_INNATE_IMMUNE_SYSTEM (Internal Immune System Related Pathways)

[0081] (3) REACTOME_IMMUNE_SYSTEM (Immune System Related Pathways)

[0082] (4) REACTOME_SIGNALING_BY_GPCR (G receptor protein signaling pathway)

[0083] (5) REACTOME METABOLISM OF PROTEINS (Protein metabolism-related pathways)

[0084] (6) REACTOME_CELLULAR_RESPONSES TO STIMULI (Cellular response to external stimuli related pathway)

[0085] (7) WP_ADIPOGENESIS_GENES (Adipogenesis-related pathways)

[0086] The results above show that the active pathways of α cells in T1D and T2D diabetic mice are different from those in normal mice, which helps to improve the pathological understanding of the two types of diabetes.

[0087] After obtaining the above results, a stability screening algorithm was used. The parameter of the Overlapping Group Lasso penalty term was λ, and its parameter screening range was set to [0.01, 0.02, ... 0.1]. For each value of λ, all α data in the T1D or T2D portion were sampled 20 times, with the sample size in each sampling being the same as the number of normal mouse cells (i.e., 1373 cells). Then, regression was performed on the gene expression data of 2746 cells, and the result after each sampling fitting was... If ||v k If ||1=0, then the path is not selected; if ||v k If ||1≠0, then the path is considered selected.

[0088] For the k-th path, a frequency function with respect to parameter λ can be obtained: p k (λ). This function value represents the ||v|| of the pathway in the 20-sample fitting model at different λ values. k The probability that ||1≠0, the value of this function is related to Images such as Figure 5a and Figure 5b As shown:

[0089] Figure 5a and Figure 5b This is called the Frequency Path for each variable. Figure 5a The results obtained by using a stability selection algorithm on T1D and normal mouse cells are shown. Figure 5aThe results are shown for T2D and normal mouse cells obtained using the same algorithm. The x-axis represents... The vertical axis is p k (λ), where black dashed lines represent paths not selected, colored solid lines represent selected paths, and red dashed lines represent frequency filtering thresholds. After plotting the Frequency Path for each variable, a filtering threshold can be determined to filter the paths: for example, in Figure 5, the filtering threshold is set to 0.7. If the maximum filtering frequency of the k-th path within the range of the λ parameter is max(p k If (λ) > 0.7, the pathway is considered a prominent pathway, meaning that the pathway has a significant difference in activity in α cells of normal mice and diabetic mice (T1D / T2D).

[0090] The stability selection algorithm offers improvements over the previous method of simply adjusting the parameter λ. For T1D, the selected pathways are consistent with those obtained through conventional Overlapping Group Lasso, demonstrating stable and significant differences in activity of these pathways in T1D mouse α cells, thus enhancing the model's interpretability. For T2D, the selected pathways differ somewhat from those obtained through conventional Overlapping Group Lasso. Using the stability selection algorithm, a new pathway, REACTOME_GENE_EXPRESSION_TRANSCRIPTION, was identified. This pathway primarily guides the production of mRNA and proteins, suggesting a significant difference in transcriptional activity between T2D and normal mouse α cells.

Claims

1. A single cell pathway analysis method developed based on Overlapping Group Lasso, characterized in that: The gene expression data of each cell is regressed with the cell state using the Logistics form of Overlapping Group Lasso model; the Stability Selection algorithm is used to regress the data by multiple sampling and the pathways with the most stable screening frequency are reserved as the final result; Specifically, the steps are as follows: Step one: mathematically express and pretreat the input cell gene expression matrix, pathway data and cell label; Step two: embed the pathway information in the pretreated gene expression data, and then use the Logistics form of Overlapping Group Lasso model to model the data; Step three: according to the model in step two, select a suitable parameter range and use the proximal gradient descent algorithm to solve the regression model to obtain the gene expression coefficient of each pathway; Step four: add the stability selection algorithm on the basis of the algorithm in step three, repeatedly sample the data and test the optimal parameters to evaluate the importance of each pathway, and finally obtain the gene expression coefficient of the important pathway and screen out the irrelevant pathway; where, since the pathway to which each gene belongs is known, let Let S be the sample matrix of the overlapping genes in each pathway. Then the optimization problem is expressed as: wherein, is a loss function, is a penalty term, which makes the regression coefficients of genes in the influential pathways not all zeros, and the regression coefficients of genes in the non-influential pathways all zeros, so as to screen out the pathways having influence on the cell state; and the method decomposes the effect of each gene in different pathways and does not fix the effect of the gene in all pathways. wherein, represents the latent vector of the kth pathway gene regression coefficient, which is 0 if the jth gene is not in the pathway, and otherwise is not 0: k wherein, The coefficient β of all genes is decomposed into the sum of K latent vectors.

2. The single-cell pathway analysis method based on Overlapping Group Lasso according to claim 1, characterized in that: In step one, set where m is the number of cells, n is the number of genes, row vector x i is the gene expression vector of the ith cell, x' i denotes its transpose, and the elements of the matrix are denotes the expression level of the ith cell at the jth gene, and the entire matrix is log-regularized before inputting the model.

3. The single-cell pathway analysis method based on Overlapping Group Lasso according to claim 2, characterized in that: Let G = {G0, G1,..., GK-1} be the set of K pathways The kth pathway is represented as a set All pathways are uniquely labeled as {0,..., K-1}, |G k | represents the number of genes in the kth pathway.

4. The single-cell pathway analysis method based on Overlapping Group Lasso according to claim 3, characterized in that: Only consider the case |G k |≥15; gene pathway vector where n is the number of genes, p j ∈{0, 1,..., K-1} is the pathway number that the jth gene belongs to, and this vector represents the numbers of all the pathways that the genes belong to; cell state vector m is the number of cells, y i ∈{0, 1} is the label of whether the ith cell is diseased.

5. The single-cell pathway analysis method based on Overlapping Group Lasso according to claim 1, characterized in that: In step two, the input matrix is The output is Let the regression model be as follows: wherein, represents the expression matrix of all genes in the kth pathway, i.e., the data of all genes in the kth pathway corresponding to the columns in X, represents the regression coefficient of each gene in the kth pathway, c is the constant term, is the Gaussian noise, representing variables other than gene expression that have an impact on the cell state.

6. The single-cell pathway analysis method based on Overlapping Group Lasso according to claim 1, characterized in that: Let path A and path b have overlapping genes g: If the coefficient of one of the pathway genes is all 0, and the coefficient of the other pathway gene is not 0 except that the coefficient of the gene is 0; this case indicates that the path A effect of the overlapping gene covers the path B effect; If the coefficients of the pathway genes A and B are not all 0, but the coefficients of the genes g are the same; this case indicates that the overlapping gene effect of path A and path b of the overlapping gene is entangled and not decomposed; Since y i ∈ {0,1}, a logistic regression is needed to convert the linear output into a probability: The loss function is the negative log maximum likelihood estimate: Let The optimization problem is then expressed as: Only V needs to be solved, and then The total regression coefficients are obtained, and each column vector in V represents the effect of each pathway on the diabetic state of the cell.

7. The single-cell pathway analysis method based on Overlapping Group Lasso according to claim 1, characterized in that: In step three, a proximal gradient descent algorithm is used to handle the optimization problem; in the t-th iteration, let the smoothing term Gradient descent on β: where, is the prediction result before each iteration, and then the minimum value problem of the non-smooth term is solved by the proximal operator, and the proximal operator is calculated once for each of the K genomes to constrain them:

8. The single-cell pathway analysis method based on Overlapping Group Lasso developed according to claim 1, characterized in that: In step four, stability selection is an algorithm used to screen variables for sparsity models, given p genomes and l different parameter values ​​for the Overlapping Group Lasso model. For each parameter value λ i The data is repeatedly sampled M times to fit the model, and the number of times N is selected in the M model fittings for the j-th independent variable under a certain parameter is obtained. The selection frequency of the j-th independent variable is then calculated based on N and M. According to f j (λ i Plot a line graph of each parameter value versus the screening frequency for each genome to visualize the screening process. Set a frequency screening threshold th∈[0,1] to screen variables. Then the j-th genome is considered to be a stable independent variable that has been selected.

Images