Method for detecting emotiality in GWAS data based on sparrow search algorithm

Abstract

The invention relates to the technical field of computer science and biology, in particular to a method for detecting the horizon in GWAS data based on a sparrow search algorithm, which comprises the following steps of: loading the GWAS data, and calculating the maximum horizon order and the maximum length of a column list; initializing the position of the sparrow group; calculating the fitness of sparrows in the group; updating the position of the discoverer; updating the position of the follower; generating a new sparrow; screening sparrows in the group, and detecting the epistatic property; judging whether local optimum exists or not; and generating a result. According to the method, the SNP combination related to the disease state can be comprehensively evaluated, and the analysis depth is improved; the order of the epistatic interaction significantly related to the disease state can be automatically identified, and the actual demand is better met; a new individual is generated based on the SNP weight vector in iteration, the SNP weight vector is updated in each iteration according to a detection result, and the evolution direction of the population is guided; the problem of false positive during detection of epistatic interaction is solved.

Description

Technical Field

[0001] The present invention relates to the fields of computer science and biological technology, and in particular to a method for detecting epistasis in GWAS data based on a sparrow search algorithm. Background Art

[0002] In recent years, despite significant progress in identifying genes associated with Mendelian genetic diseases, there are still daunting challenges in analyzing the genetic basis of non-Mendelian (i.e., complex diseases). This challenge mainly stems from the existence of epistasis, which significantly increases the complexity of genetic analysis. Complex diseases are not caused by a single gene, but are the result of the combined effects of multiple gene variants, which have a significant cumulative effect on the disease as a whole, although individually, their effects on individuals may be negligible. This cumulative effect is usually called epistatic interaction or multi-locus interaction.

[0003] With the development of high-throughput genotyping and sequencing technologies, it is now possible to detect the genotypes of millions of single nucleotide polymorphisms (SNPs) at the individual level. As an emerging strategy, Genome-Wide Association Studies (GWAS) have significantly promoted the understanding of the genetic basis of common and complex diseases. GWAS covers genotyping analysis of hundreds of thousands of SNPs in thousands of individuals. In these studies, the detection of epistatic interactions has provided new perspectives on disease genetics, helping us to understand these diseases more comprehensively and at the same time providing new avenues for disease prevention, diagnosis, and treatment.

[0004] In recent years, scholars have proposed a variety of strategies for detecting epistatic interactions in GWAS data from different perspectives. SNPHarvester (a filtering-based approach for detecting epistatic interactions) uses multiple approaches to screen out SNP combinations that are significantly associated with the disease, and uses strict statistical methods to identify SNP combinations that are closely associated with the disease (epistatic interactions). SNPRuler (a nonparametric learning approach based on a predictive rule learning algorithm) is a method based on predictive rule reasoning that can be used to discover epistatic interactions that are significantly associated with the disease. It is an earlier proposed method for finding epistatic interactions without a comprehensive search. BOOST (Operation-based Screening and Testing) is a method that can quickly detect epistatic interactions. The method mainly includes two steps: first, it uses a logical operation strategy to preliminarily screen SNP pairs; second, it performs precise statistical analysis on the screened SNP pairs to evaluate their significant association with specific diseases. AntEpiSeeker (a two-stage ant colony optimization algorithm) detects epistatic interactions in GWAS data through a two-stage ant colony optimization algorithm. MACOED (a multi-objective heuristic optimization methodology) is a multi-objective heuristic optimization method that integrates logistic regression and Bayesian network technology into the ant colony optimization algorithm to improve the efficiency of the algorithm in detecting epistatic interactions. FHSA-SED (Harmony Search Algorithm with two scoring functions) uses the harmony search algorithm to detect second-order epistasis on GWAS data. It enhances the algorithm's ability to detect epistasis by using K2 and Gini as the optimization objective function in the harmony search algorithm. DECMDR (the differential evolution algorithm with a classification based multifactor-dimensionality reduction) combines the differential evolution algorithm with the MDR method and uses CMDR as the fitness metric to detect epistatic interactions on GWAS data.HS-MMGKG (harmony search with MDR, MI, Gini, K2, G-test) combines the harmony search algorithm with multiple functions that measure the association between SNP combinations and diseases, and adopts a new strategy that combines p-value and MDR methods to improve the accuracy of the test results. SEE (sort, exploitation and exploration) is a new swarm intelligence optimization algorithm based on sorting, exploration and exploitation strategies. It uses eight functions as the objective function of group optimization at the same time. Through the complementarity between multiple objective functions, the algorithm's ability to detect epistatic interactions is improved. SHEIB-AGM (stochastic approach for detecting high-order epistatic interactions using bioinformation with automatic gene matrix) is a stochastic method based on automatic gene matrix. In each iteration, epistatic interactions on higher-order SNP combinations are randomly detected according to the content of the gene matrix, and the gene matrix is ​​updated according to the test results to ensure the detection ability of the algorithm. EpiMOGA (a multi-objective epistasis detection method based on a genetic algorithm) is a multi-objective genetic algorithm for epistasis detection. It uses K2 and Gini to jointly guide the search process of the genetic algorithm.

[0005] Although many algorithms for detecting epistatic interactions on GWAS data have emerged in recent years, they are still insufficient in accuracy and efficiency. The main problems are as follows: 1. Limitation of objective function: Most current methods based on swarm intelligence optimization algorithms rely on only a single objective function to evaluate the association between SNP combinations and diseases; however, the potential pathogenic mechanisms of complex diseases are extremely complex, and a single objective function is difficult to fully capture this association, and strategies for optimizing populations based on multi-objective functions are not yet mature. 2. Preset problem of epistasis order: Most algorithms require users to specify the epistasis order to be considered; given the complex and variable pathogenic mechanisms of complex diseases, it is difficult for users to accurately predict the epistasis order involved in actual diseases, which increases the difficulty of use. 3. Neglect of high-order epistasis: Existing algorithms often only consider second-order epistasis when designing, while ignoring higher-order (third-order and above) epistatic interactions; experimental verification is also mainly carried out on simulated data containing second-order epistatic interactions, resulting in the inability of these algorithms to identify high-order epistasis. 4. Fairness of objective function: The widely used objective functions such as K2, CE, and Gini in the field have the problem of "the finer the sample segmentation, the more biased the calculation results" when evaluating the association between SNP combinations and diseases. This bias prevents the algorithm from accurately measuring the relationship between SNP combinations and diseases. 5. False positive problem: In GWAS data, analysis of different pathogenic models shows that there are differences in the significance level of association with diseases; while traditional detection methods try to maintain high detection capabilities, they often find it difficult to effectively reduce the false positive rate, which results in a large amount of noise in the test results, which in turn affects the accuracy of the algorithm in identifying epistatic interactions. Summary of the invention

[0006] The purpose of the present invention is to overcome the defects of the existing technology and propose a method for detecting epistasis in GWAS data based on a sparrow search algorithm. Through a repeated iterative process, the sparrow population is guided to evolve in the direction of reducing multiple objective function values, and then by detecting some of the optimal sparrows in each generation, the epistatic interactions related to the disease are obtained.

[0007] To achieve the above object, the present invention adopts the following specific technical solutions:

[0008] The method for detecting epistasis in GWAS data based on the sparrow search algorithm provided by the present invention comprises the following steps:

[0009] S1. Load the genome-wide association study (GWAS) data and calculate the maximum epistasis order (mo) and the maximum length of the contingency table (ml);

[0010] S2, initialize the position of the sparrow group; generate n vectors of length mo, representing the positions of n sparrows in the group;

[0011] S3, calculate the fitness of sparrows in the group, including K2, CE and Gini objective functions;

[0012] S4, update the position of the discoverer; divide the sparrow group into discoverers and followers, and the discoverer has a better fitness in the group than the follower;

[0013] S5. Update the follower's position;

[0014] S6, generate new sparrows;

[0015] S7. Screen sparrows in the population to detect epistasis; use K2 function to detect epistatic interactions from 2 to mo order on mo order SNP combinations;

[0016] S8. Determine whether it is trapped in the local optimum; evaluate whether the search process has reached the local optimum by analyzing the proportion of deduplicated SNPs in the population;

[0017] S9. Generate results.

[0018] Furthermore, in step S1, the maximum epistasis order is used to define the upper limit of the order of epistatic interactions associated with the disease that can be detected, reducing the risk of method failure due to processing of too long single nucleotide polymorphism (SNP) combinations. The maximum epistasis order is calculated as follows:

[0019]

[0020] Where mo represents the maximum epistasis order, m0 is the number of normal samples in the GWAS data, and m1 is the number of diseased samples in the GWAS data;

[0021] The maximum length of the contingency table is used to limit the size of the contingency table generated during the calculation process to evaluate that the longer the length of the non-zero contingency table, the more significant the correlation between the SNP combination and the disease. The maximum length of the contingency table is calculated as follows:

[0022]

[0023] Where ml represents the maximum length of the contingency table, m0 is the number of normal samples in the GWAS data, and m1 is the number of diseased samples in the GWAS data.

[0024] Furthermore, in step S2, the expression for initializing the position of the sparrow group is as follows:

[0025] X i=s i,1 ,s i,2 ,......,s i,j ,......,s i,mo-1 ,s i,mo ;

[0026] In the formula, X i Represents the position vector of the i-th sparrow in the group, i∈[1,n], n is the number of sparrows in the group, and its value is specified by the user; vector X i Each element s in i,j Corresponding to the sth in the GWAS dataset i,j SNPs, i,j ∈[1,N], where N represents the total number of SNPs in the GWAS dataset.

[0027] Furthermore, in step S3, the K2, CE and Gini objective functions are as follows:

[0028]

[0029] Where X is the vector of sparrow positions, and Y is the disease status of the sample; K2 value (k2(X,Y)), CE value (ce(X,Y)) and Gini value (gini(X,Y)) are used to quantify the association between X and Y; XG represents the set of all possible combined genotypes corresponding to X; YG represents the set of sample status, including two values: 0 represents normal samples and 1 represents diseased samples; m x is the number of samples in the data with a specific combination of genotypes x on the SNP corresponding to X, m x,y is the number of samples with combined genotype x and sample state y; p(x,y) is m x,y The ratio of p(x) to the total number of samples in the data is m x The ratio of p(y|x) to the total number of samples in the data, while p(y|x) is m x,y With m x The ratio reflects the conditional probability that the sample state is y given the genotype x;

[0030] In order to integrate K2, CE and Gini objective functions as the optimization objectives of the group in the sparrow optimization algorithm, a rank-based sorting strategy is adopted, as follows:

[0031] According to the independent value of each objective function, each sparrow in the group is ranked and the rankK2, rankCE and rankGini values ​​of each sparrow are determined;

[0032] The rankK2, rankCE and rankGini values ​​of each sparrow are accumulated to obtain the comprehensive rank sum rankSum;

[0033] According to the rankSum value, the sparrows in the group are sorted from small to large. The lower the rankSum value, the more significant the correlation between the position vector of the corresponding sparrow and the disease status.

[0034] Further, in step S6, generating a new sparrow is specifically as follows:

[0035] First, n×sd sparrows are randomly selected from the population. For each selected sparrow a, if it is the best sparrow in the current population, the algorithm will randomly select mo SNPs from all SNPs in the GWAS data according to the probability stored in the pss vector to form a new position vector (sparrow) and include it in the population; if a is not the best individual, the algorithm will randomly select a discoverer b that is better than a, and then select half of the SNPs from the position vectors of a and b to combine into a new position vector, and add this new vector (sparrow) to the population.

[0036] Furthermore, in step S7, the sparrows are screened in the population to detect epistasis as follows:

[0037] Select n × sd top-ranked sparrows from the population;

[0038] For each sparrow's position vector X, based on the K2 value, repeatedly try to remove all noise SNPs in X;

[0039] Get a purified noise-free SNP combination R; if the length of R is greater than 1, use G-test to evaluate the significance of the association between R and the disease;

[0040] If the significance of the association between R and the disease is less than or equal to the threshold set by the user, R is recorded as the result;

[0041] For each SNP x in X, update pss[x] = pss[x] × 0.9, with the goal of reducing the algorithm's focus on SNPs that have already participated in epistasis detection in the next iteration;

[0042] F(X,Y)=(|XG|-1)×(|YG|-1)

[0043]

[0044] In the formula, g(X,Y) is the p-value of the G-test independence test between the sparrow position vector X and the sample state Y, F(X,Y) is the degree of freedom of the independence test, and m x is the number of samples with SNP combination genotype x, m yis 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, and pvalueOfG is a function that calculates its pvalue value under the chi-square distribution based on the statistic.

[0045] Further, in step S8, the judgment of whether it falls into a local optimum is as follows:

[0046] The search process is evaluated for reaching a local optimum by analyzing the proportion of unique SNPs in the population. If a local optimum is reached, the top n×sd sparrows in the population are removed; otherwise, the bottom n×sd sparrows in the population are removed. The calculation formula is as follows:

[0047]

[0048] Among them, the spasChaos index is used to judge whether the population reaches a local optimum state. If spasChaos < thresholdSpasChaos, it is considered that the population has fallen into a local optimum state; thresholdSpasChaos is one of the parameters for algorithm operation and can be specified by the user, with a default value of 0.6; n represents the number of sparrows in the population, and mo is the maximum epistatic order.

[0049] Further, in step S9, the generation of the result is as follows:

[0050] Sort all the epistasis in the result based on the significance of the G-test from strong to weak (from small to large);

[0051] Assume that a total of ne epistasis are detected. For each i ∈ [2, ne], calculate the ratio of the significance of the i-th epistasis to the significance of the (i - 1)-th epistasis, and record the i value corresponding to the combination with the largest ratio as iBiggest;

[0052] Output the epistasis sorted before iBiggest in the result to the result file as the finally detected epistasis, and ignore the epistasis after iBiggest as noise.

[0053] The present invention can achieve the following technical effects:

[0054] The method for detecting epistasis in GWAS data based on the sparrow search algorithm provided by the present invention adopts an optimization strategy based on multiple objective functions, which can comprehensively evaluate SNP combinations related to disease status and improve the analysis depth of GWAS data; it can automatically identify the order of epistasis interaction significantly related to disease status without the need for users to preset the epistasis order, which is more in line with the actual needs of GWAS research; in iteration, new individuals are generated based on SNP weight vectors, and the SNP weight vectors are updated in each iteration according to the detection results to guide the evolutionary direction of the group; and the false positive problem in detecting epistasis interaction is solved. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 It is a flowchart of a method for detecting epistasis in GWAS data based on a sparrow search algorithm provided according to an embodiment of the present invention.

[0056] Figure 2 It is a schematic diagram of a strategy for sorting a sparrow group based on rank provided according to an embodiment of the present invention.

[0057] Figure 3 It is a pseudo code diagram for updating the location of the finder according to an embodiment of the present invention.

[0058] Figure 4 It is a pseudo code diagram for updating the position of a tracker according to an embodiment of the present invention.

[0059] Figure 5 It is an introduction to the simulation data set provided according to an embodiment of the present invention.

[0060] Figure 6 are execution parameters of the algorithm provided according to an embodiment of the present invention on a simulated data set.

[0061] Figure 7 is the comparison result of different algorithms provided by the embodiment of the present invention on the DME 100 data set

[0062] Figure 8 4 is a comparison chart of results of different algorithms provided according to an embodiment of the present invention on the DME 100 data set.

[0063] Fig. 9 The comparison results of different algorithms provided by the embodiments of the present invention on the DNME 100 dataset are shown in FIG.

[0064] Fig.10 4 is a comparison chart of results of different algorithms provided according to an embodiment of the present invention on the DNME 100 dataset.

[0065] Fig.11The comparison results of different algorithms provided by the embodiments of the present invention on the DME 1000 data set are shown in FIG.

[0066] Fig.12 4 is a comparison chart of results of different algorithms provided according to an embodiment of the present invention on the DME 1000 data set.

[0067] Fig.13 The comparison results of different algorithms provided by the embodiments of the present invention on the DNME 1000 dataset are shown in FIG.

[0068] Fig.14 4 is a comparison chart of the results of different algorithms provided according to an embodiment of the present invention on the DNME 1000 data set.

[0069] Fig.15 1 is a comparison result of different algorithms provided according to an embodiment of the present invention on the DNME3 100 dataset.

[0070] Fig.16 3 is a comparison chart of the results of different algorithms provided according to an embodiment of the present invention on the DNME3 100 dataset. DETAILED DESCRIPTION

[0071] Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings. In the following description, the same modules are represented by the same reference numerals. In the case of the same reference numerals, their names and functions are also the same. Therefore, the detailed description thereof will not be repeated.

[0072] 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 accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and do not constitute a limitation of the present invention.

[0073] The embodiment of the present invention provides a method for detecting epistasis in GWAS data based on a sparrow search algorithm. The method is based on the idea of ​​a sparrow optimization algorithm and focuses on identifying epistatic interactions significantly associated with diseases in genome-wide association studies (GWAS) data. The method guides a sparrow population to evolve in a direction of reducing multiple objective function values ​​through a repeated iterative process, and then obtains epistatic interactions associated with diseases by detecting some optimal sparrows in each generation. The overall structure and execution flow of the method are as follows Figure 1 As shown, the specific steps include:

[0074] S1. Calculate the maximum epistasis order and the maximum length of the contingency table.

[0075] The method for detecting epistasis in GWAS data provided by the present invention uses the maximum epistasis order (mo) to define the upper limit of the order of epistasis interactions related to the disease that can be detected. The setting of this parameter can be specified by the user according to the research needs, or it can be automatically calculated based on the number of samples in the GWAS data. The specific calculation method is shown in Formula 1. The purpose of setting mo is to limit the length of the evaluation function when processing SNP combinations, ensuring that the average number of samples for each genotype combination is maintained at the level of the natural constant e during the entire operation of the algorithm. This strategy effectively reduces the risk of failure of the evaluation function that may be caused by processing too long SNP combinations, thereby improving the stability and reliability of the algorithm.

[0076]

[0077] Among them, mo represents the maximum epistasis order, m0 is the number of normal samples in the GWAS data, and m1 is the number of diseased samples in the GWAS data.

[0078] The maximum length of the contingency table (ml) is set to limit the size of the contingency table generated by the evaluation function during the calculation process. When using the evaluation function to analyze the correlation between SNP combinations and diseases, even for SNP combinations of the same length, some genotype combinations may lack samples, resulting in differences in the actual length of the non-zero contingency table. This difference may cause the evaluation function to produce calculation deviations on contingency tables of different lengths. Generally, the longer the length of the non-zero contingency table, the more significant the correlation between the SNP combination and the disease it reflects. In order to fairly evaluate this correlation, a mechanism is introduced to control the length of the contingency table during the evaluation process. The setting of this parameter can be specified by the user according to research needs, or it can be automatically calculated based on the number of samples in the GWAS data. The specific calculation method is shown in Formula 2.

[0079]

[0080] Where ml represents the maximum length of the contingency table, and the definitions of m0 and m1 are the same as in Formula 1.

[0081] S2. Initialize the position of the sparrow group.

[0082] Randomly generate n vectors of length mo, which represent the positions of n sparrows in the population. The position vector of each sparrow is defined according to Formula 3, which describes the composition of the position vector in detail. During the iteration of the algorithm, the position vectors of the sparrows in the population will undergo a continuous optimization process, which aims to identify epistatic interactions related to the disease.

[0083] X i =s i,1 ,s i,2,......,s i,j ,......,s i,mo-1 ,s i,mo (3);

[0084] Among them, X i Represents the position vector of the i-th sparrow in the group, and i∈[1,n], n is the number of sparrows in the group, and its value is specified by the user when running the algorithm. i Each element s in i,j Corresponding to the sth in the GWAS dataset i,j SNPs, and s i,j ∈[1,N], where N represents the total number of SNPs in the GWAS dataset.

[0085] S3. Calculate the fitness of the sparrows in the group.

[0086] The method for detecting epistasis in GWAS data provided by the present invention adopts three objective functions to evaluate the position vector of each sparrow in the process of optimizing the population. These objective functions include K2, CE and Gini, which are widely adopted when used to detect epistatic interactions in GWAS data. The detailed calculation methods of these functions are shown in Formula 4, and they measure the correlation between the sparrow position vector and the disease from multiple dimensions. The lower the value of these objective functions, the more significant the correlation between the corresponding position vector and the disease.

[0087]

[0088] Where X is the vector of sparrow positions, and Y is the disease status of the sample. The K2 value (k2(X,Y)), CE value (ce(X,Y)) and Gini value (gini(X,Y)) are used to quantify the association between X and Y. Where XG represents the set of all possible combined genotypes corresponding to X. For example, for a vector X of length 2, corresponding to 2 SNPs in the GWAS data, each SNP has 3 genotype states in our study: 0, 1, 2. XG contains all possible genotype combinations, namely {(0,0), (0,1), (0,2), (1,0), (1,1), (1,2), (2,0), (2,1), (2,2)}. YG represents the set of sample states, which only includes two values ​​in the study of this patent: 0 represents a normal sample and 1 represents a diseased sample. m x is the number of samples in the data with a specific combination of genotypes x on the SNP corresponding to X, m x,y is the number of samples with combined genotype x and sample state y. p(x,y) is m x,y The ratio of p(x) to the total number of samples in the data is m xThe ratio of p(y|x) to the total number of samples in the data, while p(y|x) is m x,y With m x This reflects the conditional probability that the sample state is y given the genotype x.

[0089] Although the three objective functions listed in Formula 4 are widely used in algorithms for detecting epistatic interactions, they have a common limitation: these functions are all calculated based on the contingency table between X and Y, and show significant sensitivity to the length of XG. Specifically, as the length of XG increases, the values ​​of these three objective functions tend to decrease. In order to solve this problem, the method for detecting epistasis in GWAS data provided by the present invention introduces an optimization strategy applicable to these three objective functions, which reduces the dependence on the length of XG by limiting the length of the contingency table to no more than ml. The specific operation is as follows: By analyzing Formula 4, we can find that the value of each objective function is obtained by cumulative summation of XG, and the smaller the value of the objective function, the stronger the association between X and Y. Therefore, during the calculation process, we sort the SNP combination genotypes on the contingency table and retain only the values ​​of the smallest ml-1 cells. For the remaining cells, we merge their samples to ensure that the maximum length of the contingency table does not exceed ml. This method effectively alleviates the bias of the objective function to the length of XG and improves the accuracy and applicability of the algorithm.

[0090] In order to effectively integrate the three objective functions as the optimization objectives of the group in the sparrow optimization algorithm, we adopt a rank-based sorting strategy. This process can be achieved by Figure 2 Get detailed instructions.

[0091] The specific steps are as follows: First, we sort each sparrow in the group according to the independent value of each objective function to determine the rankK2, rankCE and rankGini values ​​of each sparrow. Secondly, we add up the rankK2, rankCE and rankGini values ​​of each sparrow to obtain a comprehensive rank sum rankSum. Finally, we sort the sparrows in the group according to the rankSum value, from small to large. According to the definition of the three objective functions, the lower the rankSum value, the more significant the correlation between the position vector of the corresponding sparrow and the disease state when the three objective functions are comprehensively considered. Through this rank-based sorting method, sparrows with a strong correlation with the disease state can be efficiently identified and selected within the multi-objective optimization framework, thereby optimizing the performance of the algorithm.

[0092] S4. Update the location of the discoverer.

[0093] Based on the design idea of ​​the sparrow optimization algorithm, in the iteration, the sparrow group is divided into discoverers and followers. The discoverer has a better fitness in the group than the follower. Please refer to the pseudo code for updating the discoverer Figure 3 . The parameter pd is a decimal between 0 and 1, which is the proportion of discoverers in the sparrow population, and the default value is 0.4. The parameter n is the number of sparrows in the population. The parameter st is the safety threshold, which is a decimal between 0.5 and 1, and the default value is 0.8. The parameter mo is the maximum epistasis order, which is specified by the user or calculated according to Formula 2.

[0094] S5. Update the follower's position.

[0095] After updating the position vector of the discoverer, the followers in the group need to be updated. The idea of ​​the update is to let the better followers move towards the discoverer and the worse followers move in random directions, trying to let the followers find a better position. For the specific update method, please refer to Figure 4 , where the definition of parameters is the same as step S4 code.

[0096] S6. Generate n×sd new sparrows.

[0097] After updating the position vector of the followers, the idea of ​​the sparrow search algorithm is that when individuals at the edge of the group perceive a threat, they will move to the core area of ​​the group; at the same time, individuals in the center of the group will also conduct random exploration. In view of the characteristics of GWAS data, this behavior is simulated by generating some new sparrows. Specifically, the algorithm first randomly selects n×sd sparrows in the group. For each selected sparrow a, if it is the best sparrow in the current group, the algorithm will randomly select mo SNPs from all SNPs on the GWAS data according to the probability stored in the pss vector to form a new position vector (sparrow) and include it in the group. If a is not the best individual, the algorithm will randomly select a better discoverer b than a, and then select half of the SNPs from the position vectors of a and b to combine into a new position vector, and add this new vector (sparrow) to the group.

[0098] S7. Detect epistasis on some of the best sparrows in the population.

[0099] The method for detecting epistasis in GWAS data provided by the present invention adopts a strategy that uses the K2 function to detect epistatic interactions of 2 to mo orders on mo-order SNP combinations. The core idea can be summarized as follows: Consider a SNP combination X, whose K2 value is calculated by formula 4 and is recorded as k2X. When a SNP x is removed from X, a new SNP combination R is formed, and its K2 value is recorded as k2R. If x is a SNP associated with the disease, or x and other SNPs in X jointly affect the disease (showing epistasis), then k2R is greater than k2X; conversely, if x is noise, then k2R should be less than or equal to k2X.

[0100] Based on this concept, the following steps are taken to detect epistasis:

[0101] 1. Select the n×sd best sparrows from the population (the top n×sd sparrows in the population).

[0102] 2. For each sparrow’s position vector X, repeatedly try to remove all noise SNPs in X based on the K2 value.

[0103] 3. After the above process, a purified noise-free SNP combination R is obtained. If the length of R is greater than 1, the G-test (according to Formula 5) is used to evaluate the significance of the association between R and the disease.

[0104] 4. If the significance of the association between R and the disease is less than or equal to the threshold set by the user, R is recorded as the result.

[0105] 5. For each SNP x in X, update pss[x] = pss[x] × 0.9. The purpose is to reduce the algorithm's attention to SNPs that have already participated in epistasis detection in the next iteration.

[0106]

[0107] Among them, g(X,Y) is the p-value of the G-test independence test between the sparrow position vector X and the sample state Y, F(X,Y) is the degree of freedom of the independence test, and m x is the number of samples with SNP combination genotype x, m y is the number of samples with sample state y, m is the total number of samples, E(x,y) is the expected number of samples with SNP combination genotype x and sample state y, and pvalueOfG is a function that calculates the pvalue of the statistic under the chi-square distribution.

[0108] S8. Determine whether it is trapped in a local optimum.

[0109] The proportion of deduplicated SNPs in the population is analyzed to evaluate whether the search process reaches the local optimum, and this evaluation is based on the specific calculation method of Formula 6. If the algorithm reaches the local optimum, the top n×sd sparrows in the population are removed; otherwise, the bottom n×sd sparrows in the population are removed.

[0110]

[0111] Among them, the spasChaos index is used to judge whether the population reaches the local optimum state. If spasChaos < thresholdSpasChaos, it is considered that the population has fallen into the local optimum state; thresholdSpasChaos is one of the parameters for the algorithm to run and can be specified by the user, and the default value is 0.6; n represents the number of sparrows in the population, and mo is the maximum epistatic order.

[0112] S9. Generate results.

[0113] To reduce the false positives in the detection results, a new strategy is proposed to screen the results of epistasis detection. The purpose is to screen out the epistasis with relatively weak disease association in the detection results as noise, as follows:

[0114] All the epistasis in the results are sorted from strong to weak (from small to large) based on the significance of the G-test (Formula 5);

[0115] Suppose a total of ne epistasis are detected. For each i ∈ [2, ne], calculate the ratio of the significance of the i-th epistasis to the significance of the (i - 1)-th epistasis, and record the i value corresponding to the combination with the largest ratio as iBiggest;

[0116] The epistasis sorted before iBiggest in the results are output as the finally detected epistasis to the result file, and the epistasis after iBiggest are regarded as noise and ignored.

[0117] Through this strategy, while maintaining the detection accuracy of the algorithm, the false positives in the detection results are greatly reduced.

[0118] The following is illustrated with specific embodiments.

[0119] 1. Simulated datasets

[0120] To verify the detection ability of the method for detecting epistasis in GWAS data based on the sparrow search algorithm proposed by the present invention, a comparative experiment with other detection algorithms was conducted on a total of 38 simulated datasets. These simulated datasets were all generated based on different pathogenic models using the GAMETES_2.1 software, and their brief introduction is as Figure 5shown.

[0121] DME 100 dataset: contains 8 DME datasets, which are composed of 8 folders. Each folder corresponds to a different DME (disease loci with marginal effects) model, and contains 100 files generated by a DME model corresponding to the folder. The content of each file is a simulated GWAS data, which contains 1600 samples (800 diseased and 800 normal), 100 SNPs (2 pathogenic and 98 unrelated). Experiments on this simulated dataset can verify the ability of the algorithm to detect epistatic interactions (SNP combinations) related to diseases with marginal effects.

[0122] DNME 100 dataset: contains 8 DNME datasets, which are composed of 8 folders. Each folder corresponds to a different DNME (disease loci without marginal effects) model. The other contents are similar to the DME100 dataset. The pathogenic model simulated in the DNME dataset lacks marginal effects compared to the DME dataset, so it is more difficult to detect pathogenic SNP combinations. Experiments on this simulated dataset can verify the algorithm's ability to detect epistatic interactions related to diseases without marginal effects.

[0123] DME 1000 dataset: Its content is basically the same as DME 100 dataset, except that each GWAS data file in this dataset contains 1000 SNPs. Experiments on this simulated dataset can verify the ability of the algorithm to detect epistatic interactions with marginal effects related to diseases when faced with more noisy SNPs.

[0124] DNME 1000 dataset: Its content is basically the same as the DNME 1000 dataset, except that each GWAS data file in this dataset contains 1000 SNPs. Experiments on this simulated dataset can verify the ability of the algorithm to detect epistatic interactions related to diseases without marginal effects when faced with more noisy SNPs.

[0125] DNME3 100 dataset: contains 6 DNME datasets, which are composed of 6 folders. Each folder corresponds to a different DNME model (3rd order). Its content is similar to that of DNME 100 dataset, except that the DNME3100 model simulates a 3rd order pathogenic SNP combination model. Therefore, each GWAS simulation data file in the dataset contains 100 SNPs (3 pathogenic and 97 irrelevant). Experiments on this simulated dataset can further verify the algorithm's ability to detect high-order epistatic interactions related to diseases without marginal effects.

[0126] 2. Evaluation indicators

[0127] In order to compare the ability of different algorithms to detect epistatic interactions related to diseases on simulated data sets, the present invention uses F-measure and power as evaluation indicators to measure the detection ability of the algorithms, which are widely used in this field, and their calculation methods are shown in Formula 7.

[0128]

[0129] Among them, power and F-measure are functions used to evaluate the detection ability of the algorithm. For a model in the data set, there are 100 simulated data files. The power of the algorithm on the model is the proportion of pathogenic models in the file that are correctly detected after the algorithm is used to analyze the 100 files corresponding to the model. For each simulated file, when an algorithm is executed, several epistatic interactions will be detected. The algorithm considers them to be epistatics significantly associated with the disease. At this time, TP is the number of true pathogenic epistatics detected by the algorithm, FN is the number of true pathogenic epistatics that the algorithm failed to detect, FP is the number of non-pathogenic epistatics detected by the algorithm, recall is the recall rate, precision is the precision rate, and for each model in the data set, 100 F-measures will be obtained after the algorithm is executed. These values ​​are averaged as the F-measure of the algorithm on the model.

[0130] 3. Experimental results

[0131] In order to verify the detection ability of the Epi-SSA algorithm proposed in the present invention, the present invention uses AntEpiSeeker, DECMDR, HS-MMGKG, SEE, SHEIB-AGM, SNPHarvester, SNPRuler and Epi-SSA algorithm based on Figure 6 A large number of simulated data sets were analyzed using the parameters in , where AntEpiSeeker and SNPHarvester were unable to detect third-order epistasis. Figure 6As shown in , in order to fairly compare the ability of different algorithms to detect epistasis on simulated data sets, the parameter settings in the experiment follow the optimal recommendations in the literature of each algorithm. Moreover, since most of these algorithms are based on swarm intelligence algorithms to detect epistasis, when different algorithms are analyzed on the same simulated data set, the group size and number of iterations in the swarm intelligence algorithm are the same.

[0132] The experimental results on the DME 100 dataset are as follows: Figure 7 As shown, p and f represent power and F-measure respectively, and the calculation method is shown in Formula 7. The present invention uses different algorithms to analyze 8 simulated data sets on this data set. These 8 data sets are constructed by 8 different simulated pathogenic models. In most models, the detection ability of the Epi-SSA algorithm has obvious advantages. On average, on the DME 100 data set, the power of the Epi-SSA algorithm is higher than that of other algorithms and slightly lower than that of SHEIB-AGM, but the F-measure is higher than all other algorithms. This is because the Epi-SSA algorithm effectively reduces the false positives in the results in the process of generating results, and the comprehensive ability is improved. Figure 8 The results comparison of different algorithms on the DME100 dataset are shown.

[0133] The experimental results on the DNME 100 dataset are as follows: Fig. 9 As shown in the figure, the detection ability of the Epi-SSA algorithm has obvious advantages on most models. The results are similar to those on the DME 100 dataset. The Epi-SSA algorithm is slightly lower than SHEIB-AGM in power, but it is higher than other algorithms in comprehensive ability due to the reduction of false positives in the results. Further experiments have proved that in the task of detecting epistasis, the Epi-SSA algorithm has a high detection ability even when facing different treatment models. Fig.10 The following figure shows the comparison of the results of different algorithms on the DNME 100 dataset.

[0134] The experimental results on the DME 1000 dataset are as follows: Fig.11 As shown in the figure, the F-measure of the Epi-SSA algorithm has a very obvious advantage and is also higher than most other algorithms in terms of power. The experimental results show that in the detection task of epistasis models with marginal effects, even if the number of SNPs in the GWAS data increases, the detection ability of the Epi-SSA algorithm still has a great advantage over other algorithms. Fig.12 The results of different algorithms on the DME 1000 dataset are compared.

[0135] The experimental results on the DNME 1000 dataset are as follows: Fig.13As shown in the figure, the F-measure of the Epi-SSA algorithm still has a very obvious advantage, and its power is also higher than that of most other algorithms. The experimental results show that in the detection task of the epistasis model without marginal effects, even if the number of SNPs in the GWAS data increases, the detection ability of the Epi-SSA algorithm still has a great advantage over other algorithms. Fig.14 The results comparison of different algorithms on the DNME 1000 dataset are shown.

[0136] The experimental results on the DNME3100 dataset are as follows: Fig.15 As shown in the figure, the F-measure of the Epi-SSA algorithm still has a very obvious advantage, and its power is also higher than that of most other algorithms. The experimental results show that in the task of detecting third-order epistasis, the advantage of the Epi-SSA algorithm over other algorithms is more obvious. Fig.16 A comparison chart of the results of different algorithms on the DNME3 100 dataset is shown.

[0137] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, without contradiction.

[0138] Although the embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and are not to be construed as limitations of the present invention. A person skilled in the art may change, modify, replace and vary the above embodiments within the scope of the present invention.

[0139] The above specific implementations of the present invention do not constitute a limitation on the protection scope of the present invention. Any other corresponding changes and modifications made based on the technical concept of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. A method for detecting epistasis in GWAS data based on the sparrow search algorithm, characterized in that: The steps include: S1. Load the genome-wide association study (GWAS) data and calculate the maximum epistasis order (mo) and the maximum length of the contingency table (ml); S2, initialize the position of the sparrow group; generate n vectors of length mo, representing the positions of n sparrows in the group; S3, calculate the fitness of sparrows in the group, including K2, CE and Gini objective functions; S4, update the position of the discoverer; divide the sparrow group into discoverers and followers, and the discoverer has a better fitness in the group than the follower; S5. Update the follower's position; S6, generate new sparrows; S7. Screen sparrows in the population to detect epistasis; use K2 function to detect epistatic interactions from 2 to mo order on mo order SNP combinations; S8. Determine whether it is trapped in the local optimum; evaluate whether the search process has reached the local optimum by analyzing the proportion of deduplicated SNPs in the population; S9. Generate results.

2. The method for detecting epistasis in GWAS data based on the sparrow search algorithm according to claim 1, characterized in that: In step S1, the maximum epistasis order is used to define the upper limit of the order of epistatic interactions associated with the disease that can be detected, reducing the risk of method failure due to processing too long single nucleotide polymorphism (SNP) combinations. The maximum epistasis order is calculated as follows: Where mo represents the maximum epistasis order, m0 is the number of normal samples in the GWAS data, and m1 is the number of diseased samples in the GWAS data; The maximum length of the contingency table is used to limit the size of the contingency table generated during the calculation process to evaluate that the longer the length of the non-zero contingency table, the more significant the correlation between the SNP combination and the disease. The maximum length of the contingency table is calculated as follows: Where ml represents the maximum length of the contingency table, m0 is the number of normal samples in the GWAS data, and m1 is the number of diseased samples in the GWAS data.

3. The method for detecting epistasis in GWAS data based on the sparrow search algorithm according to claim 1, characterized in that: In step S2, the expression for initializing the position of the sparrow group is as follows: X i =s i,1 ,s i,2 ,......,s i,j ,......,s i,mo-1 ,s i,mo ; Where, X i Represents the position vector of the i-th sparrow in the group, i∈[1,n], n is the number of sparrows in the group, and its value is specified by the user; vector X i Each element s in i,j Corresponding to the sth in the GWAS dataset i,j SNPs, i,j ∈[1,N], where N represents the total number of SNPs in the GWAS dataset.

4. The method for detecting epistasis in GWAS data based on the sparrow search algorithm according to claim 1, characterized in that: In step S3, the K2, CE and Gini objective functions are as follows: Where X is the vector of sparrow positions, and Y is the disease status of the sample; K2 value (k2(X,Y)), CE value (ce(X,Y)) and Gini value (gini(X,Y)) are used to quantify the association between X and Y; XG represents the set of all possible combined genotypes corresponding to X; YG represents the set of sample status, including two values: 0 represents normal samples and 1 represents diseased samples; m x is the number of samples in the data with a specific combination of genotypes x on the SNP corresponding to X, m x,y is the number of samples with combined genotype x and sample status y; p(x,y) is m x,y The ratio of p(x) to the total number of samples in the data is m x The ratio of p(y|x) to the total number of samples in the data, while p(y|x) is m x,y With m x The ratio reflects the conditional probability that the sample state is y given the genotype x; In order to integrate K2, CE and Gini objective functions as the optimization objectives of the group in the sparrow optimization algorithm, a rank-based sorting strategy is adopted, as follows: According to the independent value of each objective function, each sparrow in the group is ranked and the rankK2, rankCE and rankGini values ​​of each sparrow are determined; The rankK2, rankCE and rankGini values ​​of each sparrow are accumulated to obtain the comprehensive rank sum rankSum; According to the rankSum value, the sparrows in the group are sorted from small to large. The lower the rankSum value, the more significant the correlation between the position vector of the corresponding sparrow and the disease status.

5. The method for detecting epistasis in GWAS data based on the sparrow search algorithm according to claim 1, characterized in that: In step S6, the new sparrow is generated as follows: First, randomly select n×sd sparrows from the population. For each selected sparrow a, if it is the optimal sparrow in the current population, then the algorithm will randomly select mo SNPs from all SNPs on the GWAS data according to the probability stored in the pss vector, form a new position vector (sparrow), and incorporate it into the population; if a is not the optimal individual, the algorithm will randomly select a discoverer b that is better than a, and then select half of the SNPs from the position vectors of a and b respectively to combine into a new position vector, and add this new vector (sparrow) to the population.

6. The method for detecting epistasis in GWAS data based on the sparrow search algorithm according to claim 1, characterized in that: In step S7, the specific method for screening sparrows in the population to detect epistasis is as follows: Select n×sd sparrows with the top rankings from the population; For the position vector X of each sparrow, based on the K2 value, repeatedly try to remove all noise SNPs in X; Obtain a purified noise-free SNP combination R; if the length of R is greater than 1, then use the G-test to evaluate the association significance between R and the disease; If the association significance between R and the disease is less than or equal to the threshold set by the user, then record R as the result; For each SNP x in X, update pss[x] = pss[x]×0.9, with the aim of reducing the attention of the algorithm to SNPs that have participated in epistasis detection in the next iteration; F(X,Y) = (|XG| - 1)×(|YG| - 1) In the formula, g(X,Y) is the p-value of the G-test independence test between the sparrow position vector X and the sample state Y, F(X,Y) is the degree of freedom of the independence test, and m x is the number of samples with SNP combination genotype x, m y is the number of samples with sample state y, m is the total number of samples, E(x,y) is the expected number of samples with SNP combination genotype x and sample state y, and pvalueOfG is a function that calculates the pvalue of the statistic under the chi-square distribution.

7. The method for detecting epistasis in GWAS data based on the sparrow search algorithm according to claim 1, characterized in that: In step S8, the specific method for determining whether it is trapped in a local optimum is as follows: Evaluate whether the search process has reached a local optimum by analyzing the proportion of unique SNPs in the population. If it has reached a local optimum, then remove the n×sd sparrows with the top rankings in the population; otherwise, remove the n×sd sparrows with the bottom rankings in the population. The calculation formula is as follows: Among them, the spasChaos index is used to determine whether the population has reached a local optimum state. If spasChaos < thresholdSpasChaos, it is considered that the population has fallen into a local optimum state; thresholdSpasChaos is one of the parameters for the algorithm to run and can be specified by the user, with a default value of 0.6; n represents the number of sparrows in the population, and mo is the maximum epistasis order.

8. The method for detecting epistasis in GWAS data based on the sparrow search algorithm according to claim 1, characterized in that: In step S9, the specific method for generating the result is as follows: Sort all the epistasis in the result based on the significance of the G-test from strong to weak (from small to large); Suppose a total of ne epistasis are detected. For each i ∈ [2, ne], calculate the ratio of the significance of the i-th epistasis to the significance of the (i - 1)-th epistasis, and record the i value corresponding to the combination with the largest ratio as iBiggest; Output the epistasis ranked before iBiggest in the result to the result file as the finally detected epistasis, and regard the epistasis after iBiggest as noise and ignore it.

Images