A method for detecting SNP combinations in GWAS data based on the crested porcupine algorithm based on genetic frequency control pool

Abstract

The method of detecting SNP combinations in GWAS data based on the crowned porcupine algorithm of the genetic frequency control pool belongs to the field of intersection of computer data processing and genetics. Four different individual generation strategies are designed, namely the global guidance strategy, the collaborative perturbation strategy, the local exploration strategy and the convergence acceleration strategy; a gene sampling module with a frequency suppression mechanism, namely the genetic frequency control pool, is introduced; a method of automatically inferring the maximum detectable interaction order based on the sample size is introduced; a significance detection mechanism combining the K2 score and the G test is adopted; the method of the present invention significantly improves the computational efficiency and detection accuracy, and effectively identifies high-order SNP combinations related to complex diseases.

Description

Technical Field

[0001] The present invention belongs to the technical field of intersection of computer data processing and genetics, and particularly relates to a method for detecting SNP combinations in GWAS data using a crested porcupine algorithm based on a genetic frequency control pool. Background Art

[0002] Complex diseases, such as hypertension, coronary heart disease, and diabetes, are typically not caused by a single gene mutation but rather by the long-term interaction of multiple genetic and environmental factors. These diseases exhibit significant multifactorial and multi-locus characteristics, making their prevention, diagnosis, and treatment highly complex. To understand their genetic mechanisms, genome-wide association studies (GWAS) are widely used to identify single nucleotide polymorphisms (SNPs) that are significantly associated with disease phenotypes. GWAS, by testing associations across millions of SNPs in large sample sizes, effectively reveal the potential relationships between numerous single nucleotide polymorphisms and disease.

[0003] However, studies have found that main effect analysis based solely on a single SNP is difficult to fully explain phenotypic variation. Interactions between SNPs (i.e., epistasis) are widespread, and such interactions cannot always be captured by univariate analysis. For example, some SNP sites are not statistically significant when evaluated alone, but when analyzed in combination with other specific sites, they may show significant phenotypic correlations. This "nonlinear" genetic effect is one of the keys to identifying complex disease mechanisms in GWAS studies. Therefore, how to efficiently identify high-order SNP combinations in large-scale data has become an important issue that needs to be addressed in the current fields of genetic statistics and computational biology.

[0004] Currently, the methods used for SNP epistasis detection can be divided into three main categories:

[0005] The first category is exhaustive methods, which identify potential interaction sites by traversing and scoring all possible SNP combinations. For example, BOOST uses Boolean operations to accelerate the analysis of binary SNP combinations and is suitable for detecting second-order interactions in large samples. TEAM constructs a minimum spanning tree model to compress the search space and calculate the statistical association between combinations and phenotypes. Multifactor Dimensionality Reduction (MDR) converts multidimensional genotype information into a one-dimensional high- and low-risk classification, thereby evaluating the relationship between SNP combinations and disease. These methods have strong detection capabilities on small datasets, but their time complexity increases exponentially when searching across the entire genome, greatly limiting their practical application.

[0006] The second category involves random sampling, which samples SNP combinations by building probabilistic models. Typical examples include epiMODE, which uses a Gibbs sampling framework combined with a modular strategy to detect apparent interactions; BEAM, which uses Bayesian inference to classify SNPs into three categories: main effect, interaction, and no effect, and uses MCMC sampling to find the optimal grouping; and epiForest, which uses random forests to construct multiple decision tree models and statistically analyze the frequency of SNP interactions to assess their importance. These methods reduce search overhead to a certain extent, but their detection capabilities are limited by the initial sampling and parameter settings, making the stability of the results difficult to guarantee.

[0007] The third category is heuristic and swarm intelligence search methods, which typically introduce biological information, prior rules, or biomimetic optimization frameworks to achieve a balance between search efficiency and accuracy. SNPRuler screens representative SNP combinations by mining prediction rules; AntEpiSeeker uses an ant colony optimization algorithm to iteratively update the solution space under the guidance of pheromones and heuristic factors, initially identifying potential interaction pairs before performing local exhaustive searches; MACOED constructs a multi-objective ant colony framework, combining logistic regression and Bayesian networks to enhance the ability to detect interaction combinations through complementary objective functions. Although these methods have made some progress in search efficiency and detection performance, they still suffer from problems such as slow convergence and sensitivity to parameters, which limit their practicality on very large-scale datasets.

[0008] Current methods for detecting SNP combinations associated with complex diseases in GWAS data have the following main shortcomings:

[0009] Most detection algorithms based on swarm intelligence optimization only include one individual generation method, while complex diseases usually involve multi-level, nonlinear genetic interactions. It is difficult to fully explore potential pathogenic SNP combinations relying on a single strategy, which limits the diversity of the understanding space and the robustness of the algorithm.

[0010] Most existing methods require users to manually specify the order of epistatic interactions to be analyzed. However, in practical applications, the interaction structure of complex diseases is often unknown, making it difficult for users to accurately judge, resulting in inconvenience in use and the possibility of missing out on true pathogenic interactions.

[0011] Current mainstream algorithms mostly focus on detecting second-order epistatic combinations, with limited support for higher-order (i.e., third-order and above) interaction structures. Experimental verification is also mainly based on second-order simulation data, resulting in poor detection of high-order SNP combinations.

[0012] While swarm intelligence algorithms (such as particle swarm optimization and ant colony algorithms) possess strong global search capabilities and robustness, they also suffer from the drawback of being easily trapped in local optima in practical applications. Specifically, during the search process, information exchange between individuals can cause the swarm to quickly converge on a specific area, reducing the exploration of diverse solutions. If algorithm parameters are improperly set or there is a lack of an effective global escape mechanism, the swarm can easily linger near the local optimum, making it difficult to escape and find the global optimal solution, ultimately limiting search results to the local optimum. Summary of the Invention

[0013] The technical problem to be solved by the present invention is to provide a method for detecting SNP combinations in GWAS data based on the crested porcupine algorithm based on genetic frequency control pool, referred to as the GFP_CPO algorithm. This algorithm significantly improves the computational efficiency and detection accuracy by introducing frequency sampling management, individual multi-strategy generation and adaptive perturbation mechanism, and effectively identifies high-order SNP combinations associated with complex diseases.

[0014] The method for detecting SNP combinations based on the genetic frequency-controlled pool crested porcupine algorithm in GWAS data is named the GFP-CPO algorithm and includes the following steps, which are performed in sequence:

[0015] Step 1: Load the GWAS data into memory; save the information in binary format and use two three-dimensional arrays to manage the genotype data of the case group and the control group respectively. The first dimension of the array represents different SNP sites, the second dimension is used to distinguish various genotype variants, and the third dimension lists the individuals with the SNP genotype.

[0016] Step 2: Set the maximum upper order number mo and the maximum length parameter ml of the contingency table to determine the longest SNP combination length processed by the GFP-CPO algorithm during the search process;

[0017] Step 3: Initialize the population individuals and build a genetic frequency control pool module to record the frequency of each SNP in the population;

[0018] Step 4: Construct an improved convergence acceleration strategy, including four individual generation strategies, namely global guidance strategy, collaborative perturbation strategy, local exploration strategy and convergence acceleration strategy; design a random strategy selection mechanism to select one of the four strategies with a medium probability before generating a new individual each time; generate individuals from the genetic frequency control pool according to frequency guidance , the direction between the individual and the current optimal solution is used as the disturbance direction;

[0019] Step 5: Calculate the value of the objective function of the newly generated individual, use the K2 function as the goal of population evolution, count the occurrence frequencies of various genotype combinations under different disease states, calculate the conditional probability, and analyze the correlation between SNP combinations and diseases.

[0020] The calculation formula for the maximum upper order number mo in step 2 is:

[0021]

[0022] Where: mo represents the maximum upper order number, log represents the logarithmic function with the natural constant as the base, m0 is the number of samples in the control group, m1 is the number of samples in the case group, and min is the function that takes the minimum value;

[0023] The maximum length parameter ml of the contingency table is used to limit the number of non-zero cells in the contingency table;

[0024]

[0025] Where: ml is the maximum length of the contingency table.

[0026] The population individual initialization method described in step 3 is: set the population size parameter PopSize, determine the number of SNPs contained in a single individual based on the automatically calculated maximum upper order number mo, and calculate the required total number of SNPs TotalSNPs, using the formula:

[0027]

[0028] If the number of SNPs N is greater than or equal to TotalSNPs, the required number is randomly allocated from the SNPs to the population; if N is less than TotalSNPs, it means that the existing SNP resources cannot meet the independent needs of each individual, and repeated extraction is performed to ensure that each individual obtains the required number of SNPs;

[0029] The genetic frequency control pool module performs frequency statistics on the SNP information of all individuals in the initial population, and the number of times each SNP appears in the entire population is recorded as f i ; and normalized to frequency distribution, the formula is:

[0030]

[0031] Where p i is the normalized frequency of the ith SNP, is the sampling probability in the current population, and f i is the number of times SNP i appears in the entire initial population, and j is the sum index used to traverse the numbers of all SNP sites.

[0032] The global guidance strategy described in step 4 uses the difference vector between the current individual and the global optimal individual in the GFP_CPO algorithm to guide the individual's forward direction. The formula is:

[0033]

[0034] Where, is the scaling factor, set to 0.1, used to control the individual movement amplitude; is a random variable in the interval [0,1], introducing perturbation diversity; Represents the individual solution with the best fitness in the current iteration; represents the position vector of the i-th individual in the t-th iteration;

[0035] The collaborative perturbation strategy generates two individuals from the genetic frequency control pool according to the frequency guidance, constructs a direction vector and drives the forward direction of the current individual, and constructs a new perturbation direction based on the frequency information. The formula is:

[0036]

[0037] Where: It represents the vector representation of the two bootstrap individuals generated by the genetic frequency control pool according to the frequency in the tth iteration;

[0038] The local exploration strategy uses the relative position between the current individual and the optimal solution to construct the perturbation direction, and controls the perturbation amplitude by introducing a diffusion factor to perform directed and bounded local exploration. The formula is:

[0039]

[0040] Where: F is the diffusion factor, which is a constant that controls the amplitude of local disturbance and is a random variable in the interval [0.05, 0.5];

[0041] The convergence acceleration strategy generates individuals from the genetic frequency control pool according to the frequency guidance. , taking the direction between the individual and the current optimal solution as the perturbation direction; by introducing the convergence perturbation coefficient Control the offset amplitude and perform a fine search strategy. The formula is:

[0042]

[0043] Where: To converge the disturbance coefficient, control the degree of proximity to the frequency control direction, and take the value as a random variable in the interval [0.05, 0.1].

[0044] The calculation formula of the K2 function in step 5 is:

[0045]

[0046] Where: Indicates SNP combination and disease states A measure of the correlation between The set of genotypes representing SNP combinations and the set of disease states are {0, 1}, where 0 represents normal and 1 represents disease; Indicates a specific genotype The number of samples, Indicates that both genotype In a sick state The number of samples;

[0047] When the k2 value of the newly generated individual is less than the average k2 value, the algorithm performs subsequent analysis and optimization steps;

[0048] When the k2 value of a newly generated individual is greater than or equal to the average k2 value, the newly generated individual does not enter the subsequent analysis and optimization steps, the algorithm updates the clustering, and continues to generate new individuals for testing;

[0049] The average k2 value uses the mean function. If the k2 value is calculated for the first time, the k2 value of the current individual is set to the average k2 value. If it is not the first calculation, the average k2 value is updated using the incremental average formula, which is:

[0050]

[0051] Where oldMeanK2 is the average k2 value before the update, nofc is the number of times the k2 value has been calculated before, and k2 is the k2 value of the current individual. After each call to the mean function, nofc increases, indicating that the total number of individuals participating in the average calculation increases.

[0052] The G-test significance test is used to analyze the association between SNP combinations and diseases described in step 5. The p-value of the G-test is less than 0.05. The p-value is the probability of significance. When the p-value is less than 0.05, it indicates that the difference is significant enough, and the two variables are considered to be not independent and there is a significant association. There is a statistically significant correlation between the newly generated individual and the disease state. The newly generated individual is added to the population; the calculation formula is:

[0053]

[0054] Where: g(X,Y) is the p-value of the G-test significance test of SNP combination X and sample state Y, F(X,Y) is the degree of freedom of the independence test, mx is the number of samples with SNP combination genotype x, m y is the number of samples with sample status y, m is the total number of samples, E(x,y) is the expected number of samples with SNP combination genotype x and sample status y, pvalueOfG It is a function that calculates the p-value of a statistic under the chi-square distribution.

[0055] Through the above design scheme, the present invention can bring the following beneficial effects: the present invention designs four different individual generation strategies, namely global guidance strategy, collaborative perturbation strategy, local exploration strategy and convergence acceleration strategy. Each strategy has different information perturbation models and search characteristics, so that multiple mechanisms can be dynamically called for collaborative optimization during the algorithm operation, thereby improving adaptability to complex genetic structures and the breadth of search.

[0056] To avoid relying on user pre-defined criteria, this paper introduces a method that automatically estimates the maximum detectable interaction order based on sample size. The algorithm automatically determines the upper limit of the interaction order within resource constraints and covers all SNP combinations with an order less than or equal to that limit during the detection process, improving ease of use and comprehensiveness of results.

[0057] This study uses a significance detection mechanism combining the K² score and the G test, enabling the algorithm to systematically analyze high-order SNP combinations while retaining sensitivity to low-order interactions. This method accurately identifies high-order (>=3) pathogenic structures, enhancing the ability to model and identify complex disease mechanisms.

[0058] This paper introduces a gene sampling module with a frequency suppression mechanism, namely the genetic frequency control pool. This module records the frequency of each SNP in the generation history and dynamically adjusts its selection probability during sampling, giving SNPs with lower occurrence frequencies higher sampling priority. This expands the understanding space, improves individual diversity, and effectively alleviates the problem of the algorithm being trapped in local optimality. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] The present invention will be further described below with reference to the accompanying drawings and specific embodiments:

[0060] Figure 1 Schematic diagram of the method for detecting SNP combinations in GWAS data based on the crested porcupine algorithm of the present invention based on the genetic frequency control pool.

[0061] Figure 2 This is a comparison chart of the results of different algorithms on the Epi2 100 dataset according to the specific implementation of the present invention.

[0062] Figure 3 This is a comparison chart of the results of different algorithms on the Epi3 100 dataset according to the specific implementation of the present invention. DETAILED DESCRIPTION

[0063] The method of detecting SNP combinations in GWAS data based on the crested porcupine algorithm with genetic frequency control pool is implemented in Java. Figure 1 As shown, the following steps are included:

[0064] Step 1: Read GWAS data

[0065] The GFP_CPO algorithm was constructed, employing an efficient data-reading method to quickly load genome-wide association study (GWAS) data into memory. This method draws on the processing principles of the BOOST method, storing information in a binary encoding format and using two three-dimensional arrays to manage genotype data for case and control groups, respectively. The first dimension of the array represents different SNP loci, the second dimension distinguishes between various genotype variants, and the third dimension lists the individuals with the genotype for each SNP. This encoding structure not only saves significant memory space but also accelerates the calculation of the K² value.

[0066] Step 2: Calculate the maximum epistasis order and the maximum length of the contingency table based on the number of samples in the data

[0067] In the algorithm of the present invention, the maximum order of epistasis (mo) is defined. This parameter limits the maximum interaction dimension of detectable SNPs and reflects the length of the longest SNP combination that the algorithm can handle during the search process. The setting of this parameter depends on the number of samples in the data. When the SNP combination is too long, the number of samples corresponding to certain combined genotypes may be insufficient, thereby affecting the effectiveness of the evaluation function. To this end, the present invention uses formula (1) to calculate mo, ensuring that during the evaluation function calculation process, the number of samples corresponding to each genotype combination is not less than the minimum value based on the natural constant e, thereby preventing the function from failing due to sparse samples.

[0068] (1)

[0069] Where mo represents the maximum epistasis order, log represents the logarithmic function with the natural constant as the base, m0 is the number of samples in the control group, m1 is the number of samples in the case group, and min is the function that takes the minimum value.

[0070] To control the complexity of the contingency table generated when evaluating the association between SNP combinations and diseases, a maximum contingency table length parameter (ml) is further defined to limit the number of nonzero entries in the contingency table. In actual calculations, even if the SNP combinations have the same length, individuals may be missing in certain genotype combinations, resulting in different numbers of nonzero entries in the contingency table. Most evaluation functions are biased by the length of the contingency table, meaning that a higher correlation score is likely to be assigned when there are more nonzero entries. To mitigate this bias, an upper limit on the contingency table length is introduced into the algorithm to ensure fairness in the evaluation of different combinations.

[0071] (2)

[0072] Where ml is the maximum length of the contingency table, and the definitions of min, m0, and m1 are the same as those in formula (1).

[0073] Step 3: Initialize the population and build a genetic frequency control pool

[0074] In the present invention, the initial stage of the algorithm operation requires the completion of population initialization and the construction of the genetic frequency control pool. First, the population size parameter PopSize is set. It can be set by the user through the parameter setting, and its default value is 40. And based on the automatically calculated maximum epistasis order mo, the number of SNPs contained in a single individual is determined, thereby calculating the required total number of SNPs TotalSNPs, using the following formula:

[0075] (3)

[0076] If N (the total number of available SNPs) is greater than or equal to TotalSNPs, the required number of SNPs can be randomly allocated directly from these SNPs to the population, ensuring the uniqueness and diversity of each individual in the population. If N is less than TotalSNPs, it means that the existing SNP resources are insufficient to meet the independent needs of each individual. In this case, repeated extraction is required to ensure that each individual obtains the required number of SNPs, although this may slightly reduce genetic diversity.

[0077] After completing individual initialization, the genetic frequency control pool module is constructed. The genetic frequency control pool is used to record the frequency of each SNP in the population and serve as a reference for the SNP selection probability in the subsequent individual generation process. Its basic principle is as follows: the SNP information of all individuals in the initial population is counted, and the number of times each SNP appears in the entire population is recorded as f i . And normalized to frequency distribution, the formula is as follows:

[0078] (4)

[0079] Among them, p iis the normalized frequency of the ith SNP (i.e., its sampling probability in the current population), f i is the number of times SNP i appears in the entire initial population, and j is the sum index used to traverse the numbers of all SNP sites.

[0080] The introduction of genetic frequency control pool has the following advantages:

[0081] (1) Guide the individual generation process to focus on SNP sites that have been less sampled in history but may be potentially important, thereby increasing the breadth of exploration in the search space;

[0082] (2) While maintaining individual diversity, avoid search bias caused by frequent selection of high-frequency SNPs and prevent the algorithm from falling into local optimality. After each individual is assigned a specific SNP combination, their K2 function is calculated.

[0083] After initialization is completed, the system will calculate the K2 score of each individual as the initial fitness indicator, providing a basic evaluation basis for subsequent clustering, generation, selection and other modules.

[0084] Step 4: Select a strategy to generate new individuals

[0085] To enhance the search space coverage and improve the algorithm's ability to escape local optima, this paper proposes four individual generation strategies with different search characteristics, corresponding to the perception and response behaviors of animals in nature: a global guidance strategy, a coordinated perturbation strategy, a local exploration strategy, and a convergence acceleration strategy. Each strategy features different perturbation directions and combinations, and can be used interchangeably during the evolutionary process to guide the search process in a diversified manner.

[0086] In order to prevent the algorithm from falling into a single search path, the present invention designs a random strategy selection mechanism, which selects one of the four strategies with equal probability before generating a new individual each time.

[0087] In the global guidance strategy, the GFP_CPO algorithm uses the difference vector between the current individual and the global optimal individual to guide the individual's direction. The specific update formula is as follows:

[0088] (5)

[0089] in, is the scaling factor, set to 0.1, used to control the individual movement amplitude; is a random variable in the interval [0,1], introducing perturbation diversity; Represents the individual solution with the best fitness in the current iteration; represents the position vector of the i-th individual in the t-th iteration.

[0090] In the collaborative perturbation strategy, in order to enhance the guidance and diversity of population evolution, the algorithm directly generates two individuals from the genetic frequency control pool according to frequency guidance, constructs a direction vector and drives the direction of the current individual. This strategy combines frequency information to construct a new perturbation direction, which helps to improve the genetic guidance ability and the ability to escape the local optimum. The formula is as follows

[0091] (6)

[0092] It represents the vector representation of the two bootstrap individuals generated by the genetic frequency-controlled pool according to the frequency in the t-th iteration.

[0093] In the local exploration strategy, to improve the current individual's search ability near the optimal solution, the algorithm uses the relative position between the current individual and the optimal solution to construct the perturbation direction and introduces a diffusion factor to control the perturbation amplitude, thereby achieving directed but bounded local exploration behavior. Unlike other strategies, this strategy does not rely on other individuals or random samples. Instead, it expands the local search range by deviating from the optimal solution in the opposite direction, which helps to enhance the algorithm's local development ability and stable convergence performance. The formula is as follows:

[0094] (7)

[0095] F is the diffusion factor, which is a constant that controls the amplitude of local disturbance and is a random variable in the interval [0.05, 0.5].

[0096] In the improved convergence acceleration strategy proposed in this invention, in order to improve the search accuracy of the algorithm near the optimal solution and make full use of historical genetic information, the perturbation direction is no longer based on the unguided global random vector construction, but is generated from the genetic frequency control pool according to the frequency guidance. , the direction between the individual and the current optimal solution is used as the perturbation direction. By introducing the convergence perturbation coefficient By controlling the offset amplitude, a refined search strategy based on "frequency control history + optimal guidance" is implemented, effectively improving the stability and optimization depth of the algorithm in the later convergence stage. The update formula is as follows:

[0097] (8)

[0098] in, To converge the disturbance coefficient, control the degree of proximity to the frequency control direction, and take the value as a random variable in the interval [0.05, 0.1].

[0099] Step 5: Calculate the value of the objective function of the newly generated individual

[0100] The GFP_CPO algorithm proposed in this paper uses the K2 function as the target for population evolution. The K2 function is an effective tool for analyzing the correlation between SNP combinations and diseases, and is particularly suitable for processing genome-wide association study (GWAS) data. Based on Bayesian networks and information theory principles, it can accurately assess the dependencies between SNPs and between SNPs and diseases. Through a greedy search strategy, the K2 function can learn data-driven models from a large number of potential network structures while maintaining computational efficiency. Furthermore, the constructed Bayesian network model not only quantifies the relationships between variables but also presents them intuitively in graphical form, greatly enhancing the model's interpretability. These properties make the K2 function particularly important in genetic data analysis, contributing to a deeper understanding of how genetic variation influences disease. The lower its value, the stronger the relationship between the corresponding SNP combination and the disease.

[0101] The specific calculation method is as follows:

[0102] (9)

[0103] In the above formula, Indicates SNP combination and disease states A measure of the correlation between A set of genotypes representing all possible SNP combinations. For example, for a two-site SNP combination, where each SNP can have a genotype of 0, 1, or 2, will include the following combinations: {(0,0), (0,1), (0,2), (1,0), (1,1), (1,2), (2,0), (2,1), (2,2)}. A set representing disease states, usually {0,1}, where 0 represents normal and 1 represents disease.

[0104] Indicates a specific genotype The number of samples, Indicates that both genotype In a sick state The number of samples.

[0105] This function uses the principles of Bayesian networks to calculate the conditional probability by statistically analyzing the frequency of various genotype combinations under different disease states, thereby assessing the association between SNP combinations and disease states. This method fully considers the frequency of various genotypes and their relationship to disease states, providing an effective mathematical model for evaluating the impact of SNP combinations on disease, making it particularly suitable for revealing the genetic mechanisms of complex diseases.

[0106] When the k2 value of the newly generated individual is less than the average k2 value calculated previously (i.e. ), the algorithm will proceed to subsequent analysis and optimization steps. This means that only when a new individual shows a stronger association than the average observed previously will it be considered for further identification of possible disease-associated SNP combinations. If the k² value of a newly generated individual fails to fall below the average, that individual will not enter the subsequent optimization process, and the algorithm will directly update the cluster and continue to generate new individuals for testing.

[0107] The average k2 value is updated using the mean function. If this is the first k2 calculation (nofc is 0), the current individual's k2 value is directly set to the average k2 value. For non-initial calculations, the average k2 value is updated using the following incremental averaging formula:

[0108] (10)

[0109] In the formula, oldMeanK2 is the average k² value before the update, nofc is the number of previous k² calculations, and k² is the k² value for the current individual. With each call to the mean function, nofc increments, indicating an increase in the total number of individuals involved in the average calculation. This approach ensures that the average k² value can be dynamically adjusted based on new data, helping the algorithm optimize its exploration of disease-associated SNP combinations based on the latest information.

[0110] Step 6: Do you need to add the newly generated individuals to the population?

[0111] The statistical significance of the newly generated individuals is further evaluated by calculating their G-test values. G-test is a statistical test used to determine the correlation between SNP combinations and disease states. If the p-value of G-test is less than 0.05, this indicates that there is a statistically significant correlation between the newly generated individuals and the disease state. In this case, according to the design of the algorithm, the newly generated individuals are considered to be good enough and should be added to the population. This decision is based on the possibility that the new individuals may provide deeper insights into disease mechanisms or increase the genetic diversity of the population, thereby achieving more efficient search and better solution space exploration in the evolutionary algorithm. The specific formula is as follows:

[0112] (11)

[0113] (12)

[0114] (13)

[0115] In the formula, g(X,Y) is the p-value of the G-test independence test of SNP combination X and sample status Y, F(X,Y) is the degrees of freedom of the independence test, mx is the number of samples with SNP combination genotype x, my is the number of samples with sample status y, m is the total number of samples, E(x,y) is the expected number of samples with SNP combination genotype x and sample status y, pvalueOfG is a function that calculates the pvalue of the statistic under the chi-square distribution, and the definitions of other symbols are the same as above.

[0116] In summary, newly generated individuals are added to the population only after they pass the G-test significance test to replace those with poor performance. Their genes are also added to the gene pool to participate in future individual generation. This strategy aims to continuously optimize and update the population to improve the algorithm's overall performance and problem-solving capabilities.

[0117] In a specific embodiment, the experiments and results of the present invention on simulated data are as follows:

[0118] The simulation datasets used to study the performance of the proposed GFP_CPO algorithm are described in detail. These datasets, totaling 18, were generated using the EpiGen simulation tool under various disease model parameters. These datasets encompass complex phenotypic generation scenarios with varying interaction types, epistasis structures, and noise levels.

[0119] To comprehensively evaluate the performance of the GFP_CPO algorithm in the complex genetic interaction detection task, we conducted systematic experiments on the aforementioned 18 simulated datasets and compared its performance with several currently available SNP detection algorithms. The generation configurations for each dataset and the corresponding pathology model parameters are detailed in Table 1.

[0120] Table 1 Introduction to the simulated dataset

[0121]

[0122] Epi2 100 dataset: This dataset consists of nine sub-datasets, each corresponding to a specific second-order epistasis model setup. Specifically, the model types include three common interaction modes: exponential, joint-dominant, and multiplicative. Each model is set with three different parameter strengths (alpha = 10, 30, and 60), resulting in a total of 3 × 3 = 9 simulation model combinations. 100 independent simulation data files were generated for each model, for a total of 900 files. Each simulation file contains 2,000 individuals (1,000 diseased and 1,000 normal) and 100 SNPs (including 2 pathogenic and 98 unrelated loci). This dataset, generated using the EpiGen tool, is primarily used to evaluate the algorithm's ability to identify pathogenic SNP combinations under a second-order, no-marginal-effect interaction structure.

[0123] Epi3 100 Dataset: This dataset also consists of 12 sub-datasets and uses the same model setup as the Epi2 100 dataset. However, Epi3 100 simulates a third-order pathogenic SNP combination, meaning each data file contains three pathogenic SNPs and 97 unrelated SNPs. Each model also generates 100 simulation files, for a total of 900 files. This dataset is used to validate the algorithm's detection performance under third-order, non-marginal interaction structures. It is particularly challenging and valuable for exploring complex, higher-order disease mechanisms.

[0124] To evaluate and compare the performance of various algorithms on a simulated dataset for identifying disease-associated SNP combinations, this paper uses two evaluation metrics: F-measure and power. These two metrics are widely used in the field to measure the detection efficiency of algorithms. The specific calculation method is as follows:

[0125]

[0126] In the present invention, power and F-measure are defined as functions used to evaluate the effectiveness of the algorithm in detecting disease-related SNP combinations. For each model in the data set, 100 simulated data files are included, and the power of the algorithm for the model is determined by analyzing these 100 files and calculating the proportion of files that correctly identify the pathogenic model. In addition, for each simulated file, after the algorithm runs, multiple SNP combinations will be identified, which are considered to be significantly associated with the disease. Here, TP represents the number of pathogenic SNP combinations correctly identified by the algorithm, FN is the number of pathogenic SNP combinations not detected by the algorithm, and FP is the number of non-pathogenic SNP combinations that are incorrectly labeled as pathogenic. Among them, recall is the recall rate and precision is the precision rate. For each model in the data set, 100 F-measure values ​​will be obtained after the algorithm is completed. By averaging these values, the F-measure score of the algorithm on the model is obtained.

[0127] The present invention analyzes a large number of simulated data sets using algorithms such as AntEpiSeeker, DECMDR, HSMMGKG, SEE, and GFP_CPO to evaluate the performance of the newly proposed GFP_CPO algorithm. According to the parameters listed in Table 2, these experiments are intended to fairly compare the ability of different algorithms to detect SNP combinations. Among them, the AntEpiSeeker algorithm does not support the detection of 3rd-order SNP combinations. To ensure the fairness of the comparison, the parameter settings of all algorithms follow the best practice recommendations in their relevant literature. In addition, considering that most algorithms are based on the principle of swarm intelligence, the analysis of different algorithms on the same data set in the experiment uses the same group size and number of iterations.

[0128] Table 2 Execution parameters of the algorithm on the simulated dataset

[0129]

[0130] Table 3 shows the experimental results on the Epi2 100 dataset, where “p” and “f” represent detection power and F-measure, respectively. The specific calculation method is detailed in formula (1).

[0131] Experimental results show that the GFP_CPO algorithm outperforms the comparison algorithms in most models. Overall, GFP_CPO achieved an average power of 0.97 and an F-measure of 0.90 on this dataset, while the highest average power and F-measure among existing algorithms were 0.94 and 0.15, respectively.

[0132] The above results indicate that GFP_CPO has significant advantages in identifying pathogenic SNP combinations with marginal effects. Figure 2 Further visualization shows the performance differences of various algorithms on the Epi2 100 dataset (retain three decimal places).

[0133] Table 3 Comparison results of different algorithms on the Epi2 100 dataset

[0134]

[0135] The execution results on the Epi3 100 dataset are shown in Table 4. On average, the power and F-measure of the GFP_CPO algorithm on this dataset are 0.89 and 0.89, respectively. Among the other algorithms, the highest average power and F-measure are 0.77 and 0.74, respectively. Although the power of HSMMGKG is 0.77, which is only slightly lower than that of the GFP_CPO algorithm, its F-measure is only 0.07, which is a big gap compared with the GFP_CPO algorithm. It can be seen that although the HSMMGKG algorithm can effectively detect pathogenic SNP combinations, its false positive rate is too high compared with the GFP_CPO algorithm, and there is a lot of noise in the detection results.

[0136] Experimental results show that the GFP_CPO algorithm still has obvious advantages over other current algorithms in detecting pathogenic SNP combinations. Figure 3 A comparison chart of the results of different algorithms on Epi3 100 data is shown.

[0137] Table 4 Comparison results of different algorithms on Epi3 100 dataset

[0138]

[0139] The above examples prove that the present invention has the following advantages:

[0140] 1. Introducing the Porcupine Optimizer (CPO) and designing a multi-strategy search mechanism. This paper simulates the multi-stage defensive behavior of porcupines and introduces the Crested Porcupine Optimizer (CPO) as the core search framework. Combining four strategies—global guidance, coordinated perturbation, local exploration, and convergence acceleration—this method simulates complex behaviors to guide search direction, achieving a dynamic balance between global exploration and local exploitation. This significantly improves the algorithm's search efficiency and stability in high-dimensional SNP spaces, enhancing its robustness and ability to obtain optimal solutions in genome-wide association analysis tasks.

[0141] 2. No need to pre-set the order of epistasis, improving the algorithm's automatic adaptability. This method automatically limits the searchable combinatorial space by calculating the theoretical maximum order of epistasis under sample number constraints. This allows exploration and identification of SNP combinations of any order (less than or equal to the maximum order), avoiding the tedious manual setting of interaction orders required in traditional methods. This enhances the algorithm's applicability and practicality in actual complex disease research.

[0142] 3. Constructing a Genetic Frequency-Controlled Pool (GFP) module to guide search direction. This paper designs a SNP frequency statistics mechanism based on historically high-fitness individuals to construct a Genetic Frequency-Controlled Pool (GFP), which dynamically records the frequency of each site in the optimal solution. This frequency-controlled pool sampling mechanism is incorporated into individual generation, perturbation direction construction, and local perturbation processes to enhance the biological rationality and convergence guidance of the search, prevent the algorithm from falling into local optima, and enhance its ability to focus on highly correlated sites.

[0143] 4. A significance detection mechanism combining the K2 score and the G-test. This invention integrates the K2 score and the G-test to assess the significance of candidate SNP combinations. The K2 score focuses on the interpretability of the combination for the class label, while the G-test measures the degree of deviation from independence in the contingency table. The two complement each other to enhance sensitivity to higher-order interaction effects. The mean K2 value (meanK2) is also introduced as an evaluation baseline to effectively suppress false-positive combinations, improve the algorithm's ability to identify high-order (≥3) pathogenic structures, and enhance the accuracy of modeling complex disease interaction mechanisms.

Claims

1. A method for detecting SNP combinations in GWAS data based on the crested porcupine algorithm using a genetic frequency control pool, characterized by: Named as GFP-CPO algorithm, it includes the following steps, And the following steps are carried out in sequence: Step 1: Load the GWAS data into memory; save the information in binary format and use two three-dimensional arrays to manage the genotype data of the case group and the control group respectively. The first dimension of the array represents different SNP sites, the second dimension is used to distinguish various genotype variants, and the third dimension lists the individuals with the SNP genotype. Step 2: Set the maximum upper order number mo and the maximum length parameter ml of the contingency table to determine the longest SNP combination length processed by the GFP-CPO algorithm during the search process; Step 3: Initialize the population individuals and build a genetic frequency control pool module to record the frequency of each SNP in the population; Step 4: Construct an improved convergence acceleration strategy, including four individual generation strategies, namely global guidance strategy, collaborative perturbation strategy, local exploration strategy and convergence acceleration strategy; design a random strategy selection mechanism to select one of the four strategies with a medium probability before generating a new individual each time; generate individuals from the genetic frequency control pool according to frequency guidance , the direction between the individual and the current optimal solution is used as the disturbance direction; Step 5: Calculate the value of the objective function of the newly generated individual, use the K2 function as the goal of population evolution, count the occurrence frequencies of various genotype combinations under different disease states, calculate the conditional probability, and analyze the correlation between SNP combinations and diseases.

2. The method for detecting SNP combinations in GWAS data based on the genetic frequency control pool-based crested porcupine algorithm according to claim 1, characterized in that: The calculation formula for the maximum upper order number mo in step 2 is: ， Where: mo represents the maximum upper order number, log represents the logarithmic function with the natural constant as the base, m0 is the number of samples in the control group, m1 is the number of samples in the case group, and min is the function that takes the minimum value; The maximum length parameter ml of the contingency table is used to limit the number of non-zero cells in the contingency table; ， Where: ml is the maximum length of the contingency table.

3. The method for detecting SNP combinations in GWAS data based on the genetic frequency control pool crested porcupine algorithm according to claim 1, characterized in that: The population individual initialization method described in step 3 is: set the population size parameter PopSize, determine the number of SNPs contained in a single individual based on the automatically calculated maximum upper order number mo, and calculate the required total number of SNPs TotalSNPs, using the formula: ’ If the number of SNPs N is greater than or equal to TotalSNPs, the required number is randomly allocated from the SNPs to the population; if N is less than TotalSNPs, it means that the existing SNP resources cannot meet the independent needs of each individual, and repeated extraction is performed to ensure that each individual obtains the required number of SNPs; The genetic frequency control pool module performs frequency statistics on the SNP information of all individuals in the initial population, and the number of times each SNP appears in the entire population is recorded as f i ; and normalized to frequency distribution, the formula is: ， Where p i is the normalized frequency of the ith SNP, is the sampling probability in the current population, and f i is the number of times SNP i appears in the entire initial population, and j is the sum index used to traverse the numbers of all SNP sites.

4. The method for detecting SNP combinations in GWAS data based on the genetic frequency control pool crested porcupine algorithm according to claim 1, characterized in that: The global guidance strategy described in step 4 uses the difference vector between the current individual and the global optimal individual in the GFP_CPO algorithm to guide the individual's forward direction. The formula is: ， Where, is the scaling factor, set to 0.1, used to control the individual movement amplitude; is a random variable in the interval [0,1], introducing perturbation diversity; Represents the individual solution with the best fitness in the current iteration; represents the position vector of the i-th individual in the t-th iteration; The collaborative perturbation strategy generates two individuals from the genetic frequency control pool according to the frequency guidance, constructs a direction vector and drives the forward direction of the current individual, and constructs a new perturbation direction based on the frequency information. The formula is: ， Where: It represents the vector representation of the two bootstrap individuals generated by the genetic frequency control pool according to the frequency in the tth iteration; The local exploration strategy uses the relative position between the current individual and the optimal solution to construct the perturbation direction, and controls the perturbation amplitude by introducing a diffusion factor to perform directed and bounded local exploration. The formula is: ， Where: F is the diffusion factor, which is a constant that controls the amplitude of local disturbance and is a random variable in the interval [0.05, 0.5]; The convergence acceleration strategy generates individuals from the genetic frequency control pool according to the frequency guidance. , the direction between the individual and the current optimal solution is used as the disturbance direction; By introducing the convergence perturbation coefficient Control the offset amplitude and perform a fine search strategy. The formula is: ， Where: To converge the disturbance coefficient, control the degree of proximity to the frequency control direction, and take the value as a random variable in the interval [0.05, 0.1].

5. The method for detecting SNP combinations in GWAS data based on the genetic frequency control pool crested porcupine algorithm according to claim 1, characterized in that: The calculation formula of the K2 function in step 5 is: ， Where: Indicates SNP combination and disease states A measure of the correlation between The genotype set representing the SNP combination and the disease status set is {0, 1}, where 0 represents normal and 1 represents disease; Indicates a specific genotype The number of samples, Indicates that both genotype In a sick state The number of samples; When the newly generated individual Value is less than average When the value is , the algorithm performs the subsequent analysis and optimization steps; When the newly generated individual Value greater than or equal to average When the value is , the newly generated individuals do not enter the subsequent analysis and optimization steps, the algorithm updates the clustering, and continues to generate new individuals for detection; The average The value uses the mean function. If the first The calculation of the value of the current individual Value is set to average Value, if it is not the first calculation, average The value is updated using the incremental averaging formula, which is: ， Where: oldMeanK2 is the average before update Value, nofc is the value that has been previously The number of times the value is calculated, The current individual value; each time the mean function is called, nofc increases, indicating that the total number of individuals participating in the average calculation increases.

6. The method for detecting SNP combinations in GWAS data based on the genetic frequency control pool-based crested porcupine algorithm according to claim 1, characterized in that: The G-test significance test is used to analyze the correlation between the SNP combination and the disease described in step 5. The p-value of the G-test is less than 0.

05. The p-value is the probability of significance. When the p-value is less than 0.05, it means that the difference is significant enough, and it is considered that the two variables are not independent and there is a significant association. There is a statistically significant correlation between the newly generated individual and the disease status. The newly generated individual is added to the population. The calculation formula is: ， Where: g(X,Y) is the p-value of the G-test significance test of SNP combination X and sample state Y, F(X,Y) is the degree of freedom of the independence test, m x is the number of samples with SNP combination genotype x, m y is the number of samples with sample status y, m is the total number of samples, E(x,y) is the expected number of samples with SNP combination genotype x and sample status y, pvalueOfG It is a function that calculates the p-value of a statistic under the chi-square distribution.

Images