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

Abstract

The present application relates to the field of computer science and biological technology, and particularly relates to a method for detecting epistasis in GWAS data based on sparrow search algorithm, which comprises the following steps: loading GWAS data, calculating the maximum epistasis order and the maximum length of contingency table; initializing the position of sparrow population; calculating the fitness of sparrow in the population; updating the position of discoverer; updating the position of follower; generating new sparrow; screening sparrow in the population to detect epistasis; judging whether to fall into local optimum; and generating results. The above method can comprehensively evaluate SNP combinations related to disease state, improve analysis depth, automatically identify the order of epistasis interaction significantly related to disease state, and be more in line with actual demand; in iteration, new individuals are generated based on SNP weight vector, and SNP weight vector is updated in each iteration according to detection results to guide the evolution direction of the population; and the false positive problem in detecting epistasis interaction is solved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of computer science and biology, and particularly relates to a method for detecting epistasis in GWAS data based on sparrow search algorithm. BACKGROUND

[0002] In recent years, although significant progress has been made in identifying genes associated with Mendelian diseases, there are still great challenges in analyzing the genetic basis of non-Mendelian (i.e. complex diseases), which mainly comes from the existence of epistasis, which significantly increases the complexity of genetic analysis. Complex diseases are not caused by a single gene, but by the combined effects of multiple gene variations, which have a significant cumulative effect on the disease as a whole, although individually they may have a negligible effect on individuals. This cumulative effect is usually referred to as 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. Genome-wide association studies (GWAS) as a new strategy have significantly promoted the understanding of the genetic basis of common and complex diseases. GWAS covers genotyping analysis of tens of thousands of SNPs in thousands of individuals. In these studies, the detection of epistatic interaction provides a new perspective for disease genetics, helping us to better understand these diseases and at the same time providing new ways for the prevention, diagnosis and treatment of diseases.

[0004] In recent years, scholars have proposed a variety of strategies to detect epistatic interactions on GWAS data from different perspectives. SNPHarvester(a filtering-based approach for detecting epistaticinteractions) uses multiple approaches to screen SNP combinations significantly associated with diseases and identifies SNP combinations closely related to diseases(epistatic interactions) through strict statistical methods. 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 find epistatic interactions significantly associated with diseases. It is an early method for finding epistatic interactions without full 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 accurate 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 on 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, thereby improving the algorithm's efficiency in detecting epistatic interactions. FHSA-SED(Harmony Search Algorithm with twoscoring functions) uses the harmony search algorithm to detect second-order epistasis on GWAS data. By using K2 and Gini as the optimization objective function in the harmony search algorithm, it enhances the algorithm's ability to detect epistasis. DECMDR(thedifferential evolution algorithm with a classification based multifactor-dimensionality reduction) combines the differential evolution algorithm and the MDR method, using CMDR as the fitness measure to detect epistatic interactions on GWAS data.HS-MMGKG(harmony search with MDR, MI, Gini, K2, G-test) combines the harmony search algorithm and multiple metrics of SNP combinations associated with disease, and adopts a new strategy of combining p-value and MDR method to improve the accuracy of detection results. SEE(sort, exploitation and exploration) is a new swarm intelligence optimization algorithm based on sorting, exploration and utilization strategy, which uses eight functions as the objective function of group optimization at the same time, and improves the ability of the algorithm to detect epistatic interaction through the complementarity between multiple objective functions. SHEIB-AGM(stochastic approach for detecting high-order epistatic interactions using bioinformation with automatic gene matrix) is a random method based on automatic gene matrix, which randomly detects higher-order SNP combinations of epistatic interaction according to the content of gene matrix in each iteration, and updates the gene matrix according to the detection result 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, which uses K2 and Gini to guide the search process of genetic algorithm together.

[0005] Although many algorithms for detecting epistatic interactions on GWAS data have emerged in recent years, they still have deficiencies in accuracy and efficiency, mainly including the following problems: 1. Limitation of objective function: Most current methods based on swarm intelligence optimization algorithm only rely on a single objective function to evaluate the association between SNP combination and disease; however, the potential pathogenesis of complex diseases is extremely complex, and a single objective function is difficult to fully capture this association, and the strategy of optimizing the population based on multiple objective functions is not mature. 2. Preset problem of epistatic order: Most algorithms require users to specify the epistatic order considered; given the complex and variable pathogenesis of complex diseases, it is difficult for users to accurately predict the epistatic order involved in the actual disease, which increases the difficulty of use. 3. Neglect of high-order epistasis: Existing algorithms often only consider two-order epistasis when designed, ignoring higher-order (three-order and above) epistatic interactions; experimental verification is also mainly carried out on simulated data containing two-order epistatic interactions, resulting in insufficient ability of these algorithms to identify high-order epistasis. 4. Fairness problem of objective function: The widely used K2, CE, Gini and other objective functions in the field have the problem of "the finer the sample segmentation, the more biased the calculation result", which hinders the accurate measurement of the relationship between SNP combination and disease by the algorithm. 5. False positive problem: Analysis of different pathogenic models on GWAS data shows that there are differences in the significance level associated with the disease; while trying to maintain high detection ability, traditional detection methods often have difficulty in effectively reducing the false positive rate, which leads to a large amount of noise in the detection results, thereby affecting the accuracy of the algorithm in identifying epistatic interactions. SUMMARY

[0006] The purpose of the present application is to overcome the defects of the prior art and provide a method for detecting epistasis in GWAS data based on sparrow search algorithm, which guides the sparrow population to evolve in the direction of reducing multiple objective function values through repeated iteration process, and then obtains the epistatic interactions related to the disease by detecting part of the optimal sparrows of each generation.

[0007] To achieve the above purpose, the following specific technical solutions are adopted in the present application:

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

[0009] S1, load Genome-Wide Association Studies (GWAS) data, calculate the maximum epistatic order (mo) and the maximum length of contingency table (ml);

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

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

[0012] S4, update the position of the discoverer; divide the sparrow population into discoverers and followers, and the discoverers have better fitness in the population than the followers;

[0013] S5, update the position of the follower;

[0014] S6, generate new sparrows;

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

[0016] S8, judge whether it is trapped in local optimum; evaluate whether the search process reaches local optimum by analyzing the proportion of de-duplication SNPs in the population;

[0017] S9, generate results.

[0018] Further, in step S1, the maximum epistatic order is used to define the upper limit of the order of the epistatic interaction related to the disease that can be detected, to reduce the risk of method failure due to processing too long single nucleotide polymorphism (SNP) combinations, and the calculation formula of the maximum epistatic order is as follows:

[0019]

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

[0021] The maximum length of the contingency table is used to limit the size of the contingency table generated in 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, and the calculation formula of the maximum length of the contingency table is as follows:

[0022]

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

[0024] Further, in step S2, the expression of initializing the position of the sparrow population 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 ith sparrow in the population, i∈[1,n], n is the number of sparrows in the population, whose value is specified by the user; each element s i in the vector X i,j corresponds to the s i,j th SNP in the GWAS dataset, s i,j ∈[1,N], N represents the total number of SNPs in the GWAS dataset.

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

[0028]

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

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

[0031] According to the independent value of each objective function, each sparrow in the population is sorted to determine the rankK2, rankCE and rankGini values of each sparrow;

[0032] The rankK2, rankCE and rankGini values of each sparrow are added up to obtain the comprehensive rank sum rankSum.

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

[0034] Further, in step S6, the new sparrow is generated as follows:

[0035] First, n x sd sparrows are randomly selected in the population. For each selected sparrow a, if it is the optimal 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, form a new position vector (sparrow), and add it to the population; if a is not the optimal individual, the algorithm will randomly select a better finder b 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.

[0036] Further, in step S7, the sparrow in the population is screened for epistasis as follows:

[0037] Select n x sd sparrows with the highest ranking from the population;

[0038] For each sparrow position vector X, based on the K2 value, repeatedly attempt 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 association significance of R and the disease;

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

[0041] For each SNP x in X, update pss[x] = pss[x] x 0.9, the purpose of which is to reduce the algorithm's attention to SNPs that have participated in epistasis detection in the next iteration;

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

[0043]

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

[0045] Further, in step S8, whether to fall into a local optimum is judged as follows:

[0046] By analyzing the proportion of the de-duplication SNPs in the population, whether the search process reaches a local optimum is evaluated, if a local optimum is reached, the top n×sd sparrows in the population are removed, otherwise, the last n×sd sparrows in the population are removed, and the calculation formula is as follows:

[0047]

[0048] Wherein, 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 of the algorithm, which 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 epistasis order.

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

[0050] All the epistases in the result are sorted from strong to weak (from small to large) based on the G-test significance;

[0051] Suppose that ne epistases are detected, for each i∈[2, ne], the ratio of the significance of the i-th epistasis to the significance of the i-1-th epistasis is calculated, and the i value corresponding to the maximum ratio is recorded as iBiggest;

[0052] The epistases in the result before iBiggest are output to the result file as the final detected epistases, and the epistases after iBiggest are ignored as noise.

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

[0054] The method for detecting epistasis in GWAS data based on the sparrow search algorithm provided by the application adopts an optimization strategy based on multiple objective functions, can comprehensively evaluate SNP combinations related to a disease state, and improves the depth of analysis of GWAS data; can automatically identify the order of epistatic interaction significantly related to the disease state, does not need to preset the order of epistasis, and is more in line with the actual needs of GWAS research; generates new individuals based on the SNP weight vector in iteration, and updates the SNP weight vector according to the detection result in each iteration, guiding the evolution direction of the population; and solves the false positive problem when detecting epistatic interaction. BRIEF DESCRIPTION OF DRAWINGS

[0055] Figure 1 is a flowchart of the method for detecting epistasis in GWAS data based on the sparrow search algorithm provided by the embodiment of the application.

[0056] Figure 2 is a strategy diagram for sorting the sparrow population based on the rank provided by the embodiment of the application.

[0057] Figure 3 is a pseudo code diagram for updating the position of the discoverer provided by the embodiment of the application.

[0058] Figure 4 is a pseudo code diagram for updating the position of the tracker provided by the embodiment of the application.

[0059] Figure 5 is an introduction to the simulation data set provided by the embodiment of the application.

[0060] Figure 6 is the execution parameter of the algorithm on the simulation data set provided by the embodiment of the application.

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

[0062] Figure 8 is a comparison diagram of the results of different algorithms on the DNME 100 data set provided by the embodiment of the application.

[0063] Figure 9 is the comparison result of different algorithms on the DNME 100 data set provided by the embodiment of the application

[0064] Figure 10 is a comparison diagram of the results of different algorithms on the DNME 100 data set provided by the embodiment of the application.

[0065] Figure 11is a comparison result of different algorithms on the DME 1000 dataset according to an embodiment of the present application

[0066] Figure 12 is a comparison result of different algorithms on the DME 1000 dataset according to an embodiment of the present application

[0067] Figure 13 is a comparison result of different algorithms on the DME 1000 dataset according to an embodiment of the present application

[0068] Figure 14 is a comparison result of different algorithms on the DME 1000 dataset according to an embodiment of the present application

[0069] Figure 15 is a comparison result of different algorithms on the DME 1000 dataset according to an embodiment of the present application

[0070] Figure 16 is a comparison result of different algorithms on the DME 1000 dataset according to an embodiment of the present application DETAILED DESCRIPTION

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

[0072] In order to make the objectives, technical solutions and advantages of the present application clearer, further detailed descriptions will be given below in combination with the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and do not constitute a limitation on the present application.

[0073] The embodiment of the present application provides a method for detecting epistasis in GWAS data based on sparrow search algorithm, which is established on the basis of the idea of sparrow optimization algorithm, focuses on identifying significant epistatic interactions related to diseases in genome-wide association study (GWAS) data, and guides the sparrow population to evolve in the direction of reducing the values of multiple objective functions through repeated iteration process, and then obtains the epistatic interactions related to diseases by detecting part of the optimal sparrows of each generation, the overall structure and execution process of the method are as shown in Figure 1 The method specifically includes the following steps:

[0074] S1, calculate the maximum epistatic order and the maximum length of the contingency table.

[0075] The method for detecting epistasis in GWAS data provided by the application uses the maximum epistasis order (mo) to define the upper limit of the order of epistatic interaction related to the disease that it can detect. 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, and 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, and to ensure that the average number of samples for each genotype combination is maintained at the level of the natural constant e during the entire running process of the algorithm. This strategy effectively reduces the risk of evaluation function failure that may be caused by processing too long SNP combinations, thereby improving the stability and reliability of the algorithm.

[0076]

[0077] wherein 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 (ml) of the contingency table is set to limit the size of the contingency table generated by the evaluation function during the calculation process. When analyzing the correlation between SNP combinations and diseases using the evaluation function, even for SNP combinations of the same length, due to the lack of samples for certain genotype combinations, the actual length of the non-zero contingency table may differ. This difference may cause the evaluation function to produce calculation bias 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 the research needs, or it can be automatically calculated based on the number of samples in the GWAS data, and the specific calculation method is shown in formula 2.

[0079]

[0080] wherein ml represents the maximum length of the contingency table, and m0 and m1 are defined the same as in formula 1.

[0081] S2, initialize the position of the sparrow population.

[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 process of the algorithm, the position vectors of the sparrows in the population will undergo continuous optimization, 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 This represents the position vector of the i-th sparrow in the group, where i ∈ [1, n], and n is the number of sparrows in the group, the value of which is specified by the user when running the algorithm. Vector X i Each element s in i,j Corresponding to the sth in the GWAS dataset i,j There are 1 SNP, and s i,j ∈[1,N], where N represents the total number of SNPs in the GWAS dataset.

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

[0086] The method for detecting epistasis in GWAS data provided by this invention employs three objective functions to evaluate the position vector of each sparrow during population optimization. These objective functions include K2, CE, and Gini, which are widely adopted for detecting epistasis interactions in GWAS data. The detailed calculation methods for these functions are shown in Equation 4, and they measure the association between sparrow position vectors and disease from multiple dimensions. The lower the values ​​of these objective functions, the more significant the association between the corresponding position vector and the disease.

[0087]

[0088] Here, X is the vector representing the sparrow's location, and Y is the disease state of the sample. The correlation between X and Y is quantified using K2 (k2(X,Y)), CE (ce(X,Y)), and Gini (gini(X,Y)). XG represents the set of all possible genotype combinations corresponding to X. For example, for a vector X of length 2, corresponding to two SNPs in GWAS data, each SNP has three genotype states in our study: 0, 1, and 2. XG contains all possible genotype combinations, i.e., {(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 in this patent study only includes two values: 0 for normal samples and 1 for diseased samples. x m is the number of samples in the data that have a specific combination of genotypes x on the SNP corresponding to X. x,y This represents the number of samples with the combinatorial genotype x and sample state y. p(x,y) is the m x,y The ratio of m to the total number of samples in the data, p(x) is m xThe ratio of the total number of samples in the data, while p(y|x) is the ratio of m x,y The ratio of m x , which reflects the conditional probability of the sample state being y given the genotype x.

[0089] Although the three objective functions listed in formula 4 are widely adopted 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 the three objective functions tend to decrease. To address this issue, the method provided by the present application for detecting epistasis in GWAS data introduces an optimization strategy applicable to the 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, in the calculation process, we sort the SNP combination genotypes on the contingency table and only retain 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, improving the accuracy and applicability of the algorithm.

[0090] To effectively integrate the three objective functions as the optimization objectives of the population in the sparrow optimization algorithm, we adopt a ranking-based sorting strategy. This process can be described in detail by Figure 2 .

[0091] The specific steps are as follows: First, we sort each sparrow in the population according to the independent value of each objective function, thereby determining the rankK2, rankCE and rankGini values of each sparrow. Second, we accumulate the rankK2, rankCE and rankGini values of each sparrow to obtain a comprehensive rank sum rankSum. Finally, we sort the sparrows in the population 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 association between the position vector of the corresponding sparrow and the disease state in the comprehensive consideration of the three objective functions. Through this ranking-based sorting method, we can efficiently identify and select sparrows with strong association with the disease state within the multi-objective optimization framework, thereby optimizing the performance of the algorithm.

[0092] S4, update the position of the discoverer.

[0093] Based on the design principles of the sparrow optimization algorithm, the sparrow population is divided into discoverers and followers during iteration. Discoverers have better fitness than followers in the population. For the pseudocode for updating the discoverer, please refer to [reference needed]. Figure 3 The parameter `pd` is a decimal between 0 and 1, representing the proportion of finders in the sparrow flock, with a default value of 0.4. The parameter `n` is the number of sparrows in the flock. The parameter `st` is the safety threshold, a decimal between 0.5 and 1, with a default value of 0.8. The parameter `mo` is the maximum epigenetic order, specified by the user or calculated according to Formula 2.

[0094] S5, Update the position of the followers.

[0095] After updating the discoverer's position vector, the followers in the group need to be updated. The idea behind this update is to move the better followers towards the discoverer and the weaker followers in random directions, attempting to help them find better positions. For specific update methods, please refer to [reference needed]. Figure 4 The definition of the parameters is the same as in step S4.

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

[0097] After updating the follower position vectors, the sparrow search algorithm works by the idea that when individuals at the edge of the group sense a threat, they move towards the core area of ​​the group; simultaneously, individuals at the center of the group also conduct random exploration. To simulate this behavior, which is tailored to the characteristics of GWAS data, some new sparrows are generated. Specifically, the algorithm first randomly selects n×sd sparrows from the group. For each selected sparrow a, if it is the best sparrow in the current group, the algorithm randomly selects mo SNPs from all SNPs in the GWAS data based on the probabilities stored in the pss vector, forming a new position vector (sparrow) and incorporating it into the group. If a is not the best individual, the algorithm randomly selects a discoverer b that is better than a, and then selects half of the SNPs from the position vectors of a and b respectively to combine them into a new position vector, which is then added 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 application adopts a strategy of detecting 2 to mo-order epistatic interactions on mo-order SNP combinations using K2 function. The core idea can be summarized as follows: considering a SNP combination X, the K2 value of which is calculated by formula 4 and recorded as k2X. When a SNP x is removed from X, a new SNP combination R is formed, and the K2 value of R is recorded as k2R. If x is a SNP related to the disease, or x and other SNPs in X jointly affect the disease (showing epistasis), then k2R is greater than k2X; on the contrary, if x is noise, then k2R should be less than or equal to k2X.

[0100] According to this idea, the following steps are taken to detect epistasis:

[0101] 1. Select n×sd sparrows with the best performance 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, use G-test to evaluate the association significance between R and the disease according to formula 5.

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

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

[0106]

[0107] where g(X,Y) is the p-value value of G-test independence test of sparrow position vector X and sample state Y, F(X,Y) is the degree of freedom of independence test, 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 pvalue of G is a function of calculating the pvalue value of the statistic under the chi-square distribution.

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

[0109] The proportion of the de-duplicated SNPs in the population is analyzed to evaluate whether the search process reaches a local optimum, which is based on the specific calculation method of formula 6. If the algorithm reaches a local optimum, the top n x sd sparrows in the population are removed, otherwise, the bottom n x sd sparrows in the population are removed.

[0110]

[0111] wherein the spasChaos index is used to determine 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 of the algorithm, which 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 order of epistasis.

[0112] S9, generating a result.

[0113] In order to reduce false positives in the detection result, a new strategy is proposed to screen the results of epistasis detection, the purpose of which is to screen out the epistasis with relatively weak disease association in the detection result as noise, and the specific steps are as follows:

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

[0115] Suppose that ne epistasis are detected, for each i∈[2, ne], the ratio of the significance of the ith epistasis to the significance of the (i-1)th epistasis is calculated, and the value of i corresponding to the maximum ratio is recorded as iBiggest;

[0116] The epistasis in the result before iBiggest is sorted as the final detected epistasis and output to the result file, and the epistasis after iBiggest is ignored as noise.

[0117] Through this strategy, the detection accuracy of the algorithm is maintained while greatly reducing the false positives in the detection result.

[0118] The following will be described in conjunction with specific embodiments.

[0119] 1. Simulated data set

[0120] In order to verify the detection ability of the method for detecting epistasis in GWAS data based on the sparrow search algorithm proposed in the present application, comparative experiments with other detection algorithms are carried out on a total of 38 simulated data sets, which are generated based on different pathogenic models using GAMETES_2.1 software, and their brief introduction is as follows Figure 5shown.

[0121] DME 100 Dataset: contains 8 DME datasets, which are organized in 8 folders, each of which corresponds to a different DME (disease loci with marginal effects) model, under which there are 100 files, which are generated by the DME model corresponding to the folder, and the content of each file is a simulated GWAS data, which contains 1600 samples (800 cases, 800 controls), 100 SNPs (2 causal, 98 irrelevant). Experiments on this simulated dataset can verify the algorithm's ability to detect disease-related epistatic interactions (SNP combinations) with marginal effects.

[0122] DNME 100 Dataset: contains 8 DNME datasets, which are organized in 8 folders, each of which corresponds to a different DNME (disease loci without marginal effects) model, and the rest of the content is similar to the DME 100 Dataset. The causal model simulated in the DNME dataset lacks marginal effects compared to the DME dataset, so it is more difficult to detect causal SNP combinations. Experiments on this simulated dataset can verify the algorithm's ability to detect disease-related epistatic interactions without marginal effects.

[0123] DME 1000 Dataset: its content is basically the same as the DME 100 Dataset, the only difference is that each GWAS data file in this dataset contains 1000 SNPs. Experiments on this simulated dataset can verify the algorithm's ability to detect disease-related epistatic interactions with marginal effects when faced with more noisy SNPs.

[0124] DNME 1000 Dataset: its content is basically the same as the DNME 1000 Dataset, the only difference is that each GWAS data file in this dataset contains 1000 SNPs. Experiments on this simulated dataset can verify the algorithm's ability to detect disease-related epistatic interactions 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 (3 orders), and the contents are similar to the DNME 100 dataset, the difference is that the DNME3 100 model simulates a 3-order pathogenic SNP combination model, therefore, each GWAS simulation data file in the dataset contains 100 SNPs (3 pathogenic, 97 irrelevant). Experiments on this simulated dataset can further verify the algorithm's ability to detect high-order marginal effect-free disease-related epistatic interactions.

[0126] 2. Evaluation index

[0127] In order to compare the detection ability of different algorithms on simulated datasets, the present application uses F-measure and power as evaluation indexes to measure the detection ability of the algorithm, which are widely used in this field, and the calculation method is shown in formula 7.

[0128]

[0129] Wherein, power and F-measure are functions for evaluating the detection ability of the algorithm, for a model in the dataset, there are 100 simulation data files, the power of the algorithm on the model is the proportion of correctly detected pathogenic models in the files after analyzing the 100 files corresponding to the model using the algorithm. And for each simulation file, when an algorithm is executed, it will detect a number of epistatic interactions, which the algorithm considers to be significantly associated with the disease, at this time, TP is the number of real pathogenic epistatic interactions detected by the algorithm, FN is the number of real pathogenic epistatic interactions that the algorithm fails to detect, FP is the number of non-pathogenic epistatic interactions detected by the algorithm, recall is the recall rate, precision is the precision rate, and for each model in the dataset, the algorithm will obtain 100 F-measure after execution, and the average of these values is taken 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 application, the present application analyzes a large number of simulated datasets using AntEpiSeeker, DECMDR, HS-MMGKG, SEE, SHEIB-AGM, SNPHarvester, SNPRuler and Epi-SSA algorithms based on the parameters in Figure 6 Among them, AntEpiSeeker and SNPHarvester cannot be used to detect three-order epistasis. As Figure 6As shown, in order to fairly compare the ability of different algorithms to detect hypergenes on simulated datasets, the parameter settings in the experiment followed the optimal recommendations in the literature for each algorithm. Furthermore, since most of these algorithms are based on swarm intelligence algorithms to detect hypergenes, the swarm size and number of iterations are the same when different algorithms are analyzed on the same simulated dataset.

[0132] Experimental results on the DME 100 dataset are as follows: Figure 7 As shown in Equation 7, p and f represent power and F-measure, respectively, and their calculation methods are illustrated in Equation 7. This invention analyzed eight simulated datasets on this dataset using different algorithms. These eight datasets were constructed from eight different simulated pathogenic models. On most models, the Epi-SSA algorithm showed a significant advantage in detection capability. On average, on the DME 100 dataset, the Epi-SSA algorithm's power was higher than other algorithms, slightly lower than SHEIB-AGM, but its F-measure was higher than all other algorithms. This is because the Epi-SSA algorithm effectively reduced false positives in the results during the generation process, thus improving its overall performance. Figure 8 The graph shows a comparison of the results of different algorithms on the DME100 dataset.

[0133] Experimental results on the DNME 100 dataset are as follows: Figure 9 As shown, the Epi-SSA algorithm demonstrates a significant advantage in detection capability on most models, with results similar to those on the DME 100 dataset. While Epi-SSA's power is slightly lower than SHEIB-AGM, its overall capability surpasses other algorithms due to reduced false positives. Further experiments demonstrate that in the task of detecting epithelial markers, the Epi-SSA algorithm maintains high detection capability even when facing different disease models. Figure 10 The graph shows a comparison of the results of different algorithms on the DNME 100 dataset.

[0134] Experimental results on the DME 1000 dataset are as follows: Figure 11 As shown, the Epi-SSA algorithm has a very significant advantage in F-measure and is also superior to most other algorithms in power. Its experimental results show that in the detection task of the superordinate model with marginal effects, even if the number of SNPs in the GWAS data increases, the detection capability of the Epi-SSA algorithm still has a great advantage over other algorithms. Figure 12 The graph shows a comparison of the results of different algorithms on the DME 1000 dataset.

[0135] Experimental results on the DNME 1000 dataset are as follows: Figure 13As shown, the F-measure of the Epi-SSA algorithm still has a very obvious advantage, and the power is also higher than that of most other algorithms. The experimental results show that in the detection task of the upper-order model without marginal effect, 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. Figure 14 A comparison chart of results of different algorithms on the DNME 1000 dataset is shown.

[0136] The experimental results on the DNME3 100 dataset are as shown in the table. Figure 15 As shown, the F-measure of the Epi-SSA algorithm still has a very obvious advantage, and the power is also higher than that of most other algorithms. The experimental results show that in the detection task of the three-order upper-order, the advantage of the Epi-SSA algorithm over other algorithms is more obvious. Figure 16 A comparison chart of results of different algorithms on the DNME3 100 dataset is shown.

[0137] In the description of the present specification, the description of the terms "one embodiment", "some embodiments", "an example", "a specific example", or "some examples" and the like means that the specific features, structures, materials or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present application. In the present specification, the illustrative description of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any appropriate manner in any one or more embodiments or examples. In addition, different embodiments or examples described in the present specification and the features of different embodiments or examples can be combined and modified by those skilled in the art without contradiction.

[0138] Although the embodiments of the present application have been shown and described above, it should be understood that the above embodiments are exemplary and should not be construed as limiting the present application, and those skilled in the art can make changes, modifications, replacements and variations to the above embodiments within the scope of the present application.

[0139] The specific embodiments of the present application described above do not constitute a limitation on the scope of protection of the present application. Any various other corresponding changes and modifications made according to the technical concept of the present application shall be included in the scope of protection of the claims of the present application.

Claims

1. A method for detecting epistasis in GWAS data based on sparrow search algorithm, characterized in that, Comprise the following steps: S1, load the whole genome association study data, calculate the maximum order of epistasis mo and the maximum length of contingency table ml; S2, initialize the position of sparrow population; generate n vectors with length mo, representing the position of n sparrows in the population; S3, calculate the fitness of sparrows in the population, including K2, CE and Gini objective functions; S4, update the position of the discoverer; divide the sparrow population into discoverers and followers, and the discoverers have better fitness in the population than the followers; S5, update the position of the follower; S6, generate new sparrows; First, randomly select n×sd sparrows in the population, for each selected sparrow a, if it is the optimal sparrow in the current population, the algorithm will randomly select mo SNPs from all SNPs in GWAS data according to the probability stored in the pss vector to form a new position vector, that is, a sparrow, and add it to the population; if a is not the optimal individual, the algorithm will randomly select a better discoverer b from 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, that is, a sparrow, to the population; S7, screen the sparrow in the population to detect epistasis; use K2 function to detect 2 to mo order epistatic interaction on mo order SNP combination; S8, judge whether it is trapped in local optimum; evaluate whether the search process reaches local optimum by analyzing the proportion of de-duplication SNPs in the population; The judgment of whether it is trapped in local optimum is as follows: If the local optimum is reached, remove the top n×sd sparrows in the population, otherwise, remove the last n×sd sparrows in the population, the calculation formula is as follows: ； Wherein, spasChaos index is used to judge whether the population reaches local optimum state, if spasChaos<thresholdSpasChaos, it is considered that the population has been trapped in local optimum state; thresholdSpasChaos is one of the parameters of the algorithm, which can be specified by the user, the default value is 0.6; n represents the number of sparrows in the population, mo is the maximum order of epistasis; S9, generate results.

2. The method for detecting epistasis in GWAS data based on sparrow search algorithm according to claim 1, characterized in that, In step S1, the maximum order of epistasis is used to define the upper limit of the order of epistatic interaction related to the disease that can be detected, to reduce the risk of method failure due to processing too long single nucleotide polymorphism combination, the calculation formula of the maximum order of epistasis is as follows: ； In the formula, mo represents the maximum order of epistasis, m0 is the number of normal samples in GWAS data, and m1 is the number of diseased samples in GWAS data; The maximum length of contingency table is used to limit the size of the contingency table generated in the calculation process, to evaluate that the longer the length of non-zero contingency table, the more significant the correlation between SNP combination and disease, the calculation formula of the maximum length of contingency table is as follows: ； In the formula, ml represents the maximum length of contingency table, m0 is the number of normal samples in GWAS data, and m1 is the number of diseased samples in GWAS data. 3.The method for detecting epistasis in GWAS data based on sparrow search algorithm according to claim 1, wherein, In step S2, the expression of the initial position of the sparrow population is as follows: ； where X i represents the position vector of the i-th sparrow in the population, i ∈ [1, n], n is the number of sparrows in the population, whose value is specified by the user; each element s i in the vector X i,j corresponds to the s i,j -th SNP in the GWAS dataset, s i,j ∈ [1, N], N represents the total number of SNPs in the GWAS dataset. 4.The method of detecting epistasis in GWAS data based on sparrow search algorithm according to claim 1, wherein, In step S3, K2, CE and Gini objective functions are as follows: ； ； ； where X is a vector of sparrow positions, Y is the disease status of samples; the association between X and Y is quantified by K2 value k2(X,Y), CE value ce(X,Y) and Gini value gini(X,Y); XG represents the set of all possible combined genotypes corresponding to X; YG represents the set of sample states, including two values: 0 represents normal samples, 1 represents diseased samples; m x is the number of samples in the data with specific combined genotype x at the SNPs corresponding to X, m x,y is the number of samples with combined genotype x and sample state y; p(x,y) is the ratio of m x,y to the total number of samples in the data, p(x) is the ratio of m x to the total number of samples in the data, and p(y|x) is the ratio of m x,y to m x , reflecting the conditional probability of sample state y given genotype x; In order to integrate K2, CE and Gini objective functions as the optimization objective of the population in the sparrow optimization algorithm, a ranking-based sorting strategy is adopted, which is as follows: According to the independent value of each objective function, each sparrow in the population is sorted to determine the rankK2, rankCE and rankGini values of each sparrow. 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 population are sorted from small to large, and the lower the rankSum value, the more significant the relevance between the position vector of the corresponding sparrow and the disease state.

5. The method for detecting epistasis in GWAS data based on sparrow search algorithm according to claim 1, characterized in that, In step S7, the sparrows in the population are screened for epistasis as follows: Select n×sd sparrows with high ranking 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, use G-test to evaluate the association significance of R and disease; If the association significance of R and disease is less than or equal to the threshold set by the user, record R as the result; For each SNP x in X, update pss[x]=pss[x]×0.9, the purpose is to reduce the attention of the algorithm to the SNP that has participated in the epistasis detection in the next iteration. ； where g(X, Y) is the p-value of G-test for independence between the Sparrow position vector X and the 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 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 pvalue of G is a function of calculating the p-value of the chi-square distribution according to the statistic.

6. The method for detecting epistasis in GWAS data based on sparrow search algorithm according to claim 1, wherein, In step S9, the result is generated as follows: Sort all epistasis in the result from strong to weak based on the significance of G-test; Assuming that 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 maximum ratio as iBiggest. The epistasis in the result before iBiggest is sorted as the final detected epistasis and output to the result file, while the epistasis after iBiggest is ignored as noise.

Images