ACOCMPMI: An Ant Colony Optimization Algorithm Based on Composite Multiscale Part Mutual Information for Detecting Epistatic Interactions

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In recent years, numerous single nucleotide polymorphisms (SNPs) associated with complex diseases have been successfully detected through genome-wide association studies (GWAS) [1]. However, the explanatory power of individual SNPs is limited in some complex diseases, such as cancer [2] and Alzheimer’s disease [3]. Epistatic interactions, broadly defined as nonlinear interactions between SNPs, have emerged as a key mechanism to overcome these limitations. Therefore, the precise detection of epistatic interactions has become a focal point of research [4–6].

Epistatic interaction detection focuses on two key aspects: interaction quantification measures and search strategies. Interaction can be likened to a specific type of association, predominantly manifesting as nonlinear direct associations. Quantifying these interactions relies on various association measures. Traditional statistical measures, such as logistic regression [7, 8], chi-square statistic [9], distance covariance [10], and Pearson’s correlation coefficient [11], are limited to quantifying nonlinear direct associations among target variables. Measures based on information entropy, which do not strictly depend on specific association forms, have gained significant attention in recent years. The mutual information (MI) and conditional mutual information (CMI) are commonly employed for quantifying nonlinear interactions among variables [12–15]. However, they may lead to overestimation and underestimation problems [16]. To precisely quantify nonlinear direct interactions, several measures have emerged, including maximum information coefficient (MIC) [17], conditional mutual inclusive information (CMI2) [18], part mutual information (PMI) [16], partial association (PA) [19] and multiscale part mutual information (MPMI) [20]. Notably, MPMI demonstrates higher accuracy compared with other measures and has not been applied to SNP data. Therefore, this study adopts MPMI and its variant to quantify interaction between SNPs.

The search strategy can be broadly categorized into three groups: exhaustive search, stochastic search, and heuristic search. Exhaustive search methods typically attempt to evaluate all possible SNP combinations within a dataset. However, the high dimensionality of GWAS data imposes a heavy computational burden on exhaustive methods [21]. Stochastic search methods are limited in the number of features they can handle [22]. Heuristic search transforms the epistatic interaction detection problem into an optimization problem. Heuristic search mainly focuses on metaheuristic optimization algorithms, such as the firefly algorithm [23], tree seed algorithm [24], tunicate swarm algorithm [25], side-blotched lizard algorithm [26], African vultures optimization algorithm (AVOA) [27], ant colony optimization (ACO) algorithm [28], symbiotic organisms search algorithm [29], spotted hyena optimizer algorithm [30], yellow saddle goatfish behavior optimization model [31], and grey wolf optimizer [32]. In this study, the ACO algorithm ($\text{ACO}\ast$) is employed for searching epistatic interactions, and an improved version of the $\text{ACO}\ast$ is presented. The $\text{ACO}\ast$ has been widely used in this field [33] and is considered one of the most promising methods among these metaheuristic optimization algorithms.

The main contributions of this work are as follows.

• A composite version of MPMI, termed CMPMI, is proposed. CMPMI is specifically designed for detecting nonlinear direct interactions in SNP datasets.

• Memory and filtering strategies are integrated into the $\text{ACO}\ast$ to improve the accuracy of epistatic interaction detection.

• Epistatic interactions are detected in a two-stage framework. In the first stage, an improved $\text{ACO}\ast$ combined with CMPMI is used to generate a candidate SNP set. In the second stage, an exhaustive search strategy and a Bayesian network (BN) score are adopted to further identify epistatic interactions within the candidate set.

2. Related Works

Various methods have been proposed to detect epistatic interactions. For instance, multifactor dimensionality reduction (MDR) [34], backward genotype-trait association (BGTA) [35], Boolean operation-based screening and testing (BOOST) [8], factored spectrally transformed linear mixed models (FaST-LMM) [36], and tree-based epistasis association mapping (TEAM) [37] are epistatic interaction detection methods based on exhaustive search strategies. Bayesian epistasis association mapping (BEAM) [22] and epistatic module detection (EpiMODE) [38] employ stochastic search strategies. BEAM integrates the Bayesian partitioning model with Markov chain Monte Carlo to assess and identify disease-associated SNPs and epistatic interactions. EpiMODE utilizes a Bayesian marker partition model alongside a Gibbs sampling strategy to detect epistatic interactions. For heuristic search methods, CINOEDV is designed to detect and visualize epistatic interactions of various orders, leveraging the particle swarm optimization algorithm and co-information measure [39]. AntEpiSeeker uses a two-stage $\text{ACO}\ast$ for identifying epistatic interactions in large datasets [40]. Similarly, MACOED is a multiobjective ACO supervised heuristic method for epistasis detection [41], and IACO applies an improved $\text{ACO}\ast$ to search for epistatic interactions [42]. MTHSA-DHEI is proposed for detecting high-order epistatic interactions based on a multitasking harmony search algorithm [43]. Building on this framework, MTHS-EE-DHEI is introduced as an enhanced variant that incorporates explicit encoding into the multitasking harmony search algorithm to further optimize the epistasis detection [44].

3. Materials and Methods

3.1. MPMI

MPMI is an innovative measure designed to quantify direct associations between target variables [20]. Unlike traditional measures, it is not confined to specific interaction forms during quantification. Furthermore, its higher accuracy and superior statistical power render it a significant advancement in this field. The MPMI between $X$ and $Y$ given $Z$ is defined as $$\begin{array}{}\text{(1)}& \text{MPMI}X;YZ=e^{2\text{MI}X;Z}\cdot e^{2\text{MI}Y;Z}\cdot \text{CMI}X;YZ+Dpxz{p}^{\ast }xz+Dpyz{p}^{\ast }yz,\end{array}$$where both $X$ and $Z$ represent SNPs and $Y$ represents the phenotype. $x$ and $z$ are genotypes of SNPs $X$ and $Z$, and $y$ is the class label of $Y$. $\text{MI}X;Z$ is the MI between $X$ and $Z$, and $\text{MI}Y;Z$ is the MI between $Y$ and $Z$. Both of them are defined as $$\begin{array}{}\text{(2)}& \text{MI}X;Z=\sum _{x,z}px,z\log \frac{px,z}{pxpz},\text{(3)}& \text{MI}Y;Z=\sum _{\text{y},z}py,z\log \frac{py,z}{pypz},\end{array}$$where $px,z$ is the joint probability distribution of $x$ and $z$, $py,z$ is the joint probability distribution of $y$ and $z$, and $\text{CMI}X;YZ$ is the CMI between $X$ and $Y$ given $Z$, which is defined as $$\begin{array}{}\text{(4)}& \text{CMI}X;YZ=\sum _{x,y,z}px,y,z\log \frac{px,y,zpz}{px,zpy,z},\end{array}$$where $px,y,z$ is the joint probability distribution of $x$, $y$, and $z$.

In addition, both $Dpxz{p}^{\ast }xz$ and $Dpyz{p}^{\ast }yz$ are the extended Kullback–Leibler divergences [16]. They are defined as $$\begin{array}{}\text{(5)}& Dpxz{p}^{\ast }xz=\sum _{x,y,z}px,y,z\log \frac{pxz}{{\sum }_{y}pxz,ypy},\text{(6)}& Dpyz{p}^{\ast }yz=\sum _{x,y,z}px,y,z\log \frac{pyz}{{\sum }_{x}pyz,xpx},\end{array}$$where $pxz$ and $pyz$ are the probability distributions of $x$ and $y$ conditioned on $z$, respectively, and $pxz,y$ is the probability distribution of $x$ conditioned on both $z$ and $y$.

3.2. $\text{ACO}\ast$

The $\text{ACO}\ast$ is a classical swarm intelligence optimization algorithm designed to solve complex combinatorial optimization problems by simulating the cooperative behavior of ant colonies [45]. The basic idea of the $\text{ACO}\ast$ is to map feasible solutions of optimization problems to paths traversed by ants. Ants tend to release more pheromones along shorter paths during their traversal. Meanwhile, pheromones guide ants in selecting subsequent paths. Ultimately, through positive feedback, all ants converge on the optimal path, which corresponds to the optimal solution of the optimization problem. The basic $\text{ACO}\ast$ primarily involves two core strategies: path selection and pheromone update.

Ants navigate paths based on a combination of pheromones and heuristic information. Typically, the probability of an ant selecting the next position from a given current position during an iteration is defined as $$\begin{array}{}\text{(7)}& P_k^{\text{ij}}t=\begin{array}{cc}\frac{{\tau }_{\text{ij}}{t}^{\alpha }{\eta }_{\text{ij}}^{\beta }}{{\sum }_{u\in M_{k}t}{\tau }_{\text{iu}}{t}^{\alpha }{\eta }_{\text{iu}}^{\beta }}\hfill & i\in M_{k}t\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},\end{array}$$where ${\tau }_{\text{ij}}t$ is the pheromone of path $i⟶j$ in iteration $t$. Similarly, ${\eta }_{\text{ij}}$ represents the heuristic information of path $i⟶j$. $\alpha$ and $\beta$ are the weight coefficients of pheromone and heuristic information, respectively, both of which are usually set to 1. $M_{k}t$ represents the set of positions that are not detected by ant $k$ in iteration $t$.

In iteration $t+1$, the pheromone of path $i⟶j$ is defined as $$\begin{array}{}\text{(8)}& {\tau }_{\text{ij}}t+1=1-\rho {\tau }_{\text{ij}}t+\Delta {\tau }_{\text{ij}}t,\text{(9)}& \Delta {\tau }_{\text{ij}}t=\sum _{k=1}^m\Delta {\tau }_{\text{ij}}^{k}t,\text{(10)}& \Delta {\tau }_{\text{ij}}^{k}t=\begin{array}{cc}\frac{Q}{S_{k}t},\hfill & \text{if}\text{ant}k\text{passing}\text{path}\text{ij}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},\end{array}$$where ${\tau }_{\text{ij}}t$ is the pheromone of path $i⟶j$ in iteration $t$, $\rho$ is an evaporation coefficient, ${\tau }_{\text{ij}}t$ is the pheromone variation of path $i⟶j$ in iteration $t$, $Q$ is a user-defined constant, and $S_{k}t$ is the path length of ant $k$ in iteration $t$.

3.3. BNs

A BN is a network structure based on a directed acyclic graph, used to represent dependencies among observed variables. In this network, nodes represent either SNPs or phenotypes, and edges connecting nodes signify causal dependencies. The K2 score, based on BN, is widely used to quantify causal dependencies between two variables.

The K2 score is derived from the Bayesian score. The Bayesian score computes the posterior probability $PMD$ of the BN model $M$ given the data $D$, which can be written as $$\begin{array}{}\text{(11)}& PMD=\frac{PDMPM}{PD},\end{array}$$where $PDM$ is the class-conditional density and $PD$ and $PM$ are the probabilities of the data $D$ and the model $M$, respectively. Building upon prior studies [41, 46, 47], in the context of a case-control study, if all variables in the directed acyclic graph are discrete, we can derive $$\begin{array}{}\text{(12)}& PDM=\prod _{i=1}^I\frac{\Gamma {\sum }_{j=1}^J{\alpha }_{ij}}{\Gamma r_i+{\sum }_{j=1}^J{\alpha }_{ij}}\prod _{j=1}^J\frac{\Gamma r_{ij}+{\alpha }_{ij}}{\Gamma {\alpha }_{ij}}\end{array}$$where $I$ is the number of combinations of SNP nodes with different genotypes, $r_i$ is the case number of SNP nodes taking the $i_{\text{th}}$ combination numbers, and $r_{\text{ij}}$ is the number of cases with phenotypes taking the $j_{\text{th}}$ state while its parents take the $i_{\text{th}}$ combination. $J$ is the state number of phenotypes. ${\alpha }_{\text{ij}}$ is the prior belief about case numbers with model nodes taking the $i_{th}$ combination and $j_{\text{th}}$ state, which is a hyperparameter when the model satisfies the Dirichlet distribution. If ${\alpha }_{\text{ij}}=1$, $PM$ and $PD$ are constants, then, $$\begin{array}{}\text{(13)}& PMD\propto \prod _{i=1}^I\frac{J-1!}{r_i+J-1!}\prod _{j=1}^J{r}_{ij}!.\end{array}$$

The Bayesian score can be transformed into the K2 score. Subsequently, the logarithmic form of the K2 score can be derived. $$\begin{array}{}\text{(14)}& K2\text{score}=\prod _{i=1}^I\frac{J-1!}{r_i+J-1!}\prod _{j=1}^J{r}_{ij}!,\text{(15)}& K2{\text{score}}_{\log }=\sum _{i=1}^I\sum _{b=1}^{r_i+1}\log b-\sum _{j=1}^J\sum _{d=1}^{r_{ij}}\log d.\end{array}$$

3.4. ACOCMPMI

Figure 1 is the flow chart of ACOCMPMI. It can be seen that the ACOCMPMI mainly consists of two parts: Stage 1 (CMPMI + improved ACO) and Stage 2 (exhaustion search + BN). Among them, Stage 1 is the highlight of ACOCMPMI.

[figure(s) omitted; refer to PDF]

3.5. CMPMI

The MPMI possesses several properties that prove beneficial for the investigation of epistatic interaction detection. For instance, (1) $\text{MPMI}X;YZ\ge 0$; (2) $\text{MPMI}X;YZ=0$ if and only if there is no direct association between $X$ and $Y$ under the condition of $Z$; (3) when both $X$ and $Y$ are independent of $Z$, $\text{MPMI}X;YZ=\text{CMI}X;YZ=\text{MI}X;Y$; (4) $\text{MPMI}X;YZ$ is a vsconstant when $X$ is directly associated with $Y$ regardless of the influence intensity of $Z$ on their association; (5) for the target variables $X$ and $Y$, $\text{MPMI}X;YZ=\text{MPMI}Y;XZ$.

It is seen that the MPMI can be regarded as “asymmetric” for the two-order SNP combination $X,Z$ and the phenotype $Y$, which is inconsistent with the basic principle of association. Hence, to capture symmetric information in the detection of two-order epistatic interactions, we define the CMPMI as $$\begin{array}{}\text{(16)}& \text{CMPMI}=\frac{\text{MPMI}X;YZ+\text{MPMI}Z;YX}{2}.\end{array}$$

CMPMI is essentially the mean form of MPMI between the involved target variables, indicating the integration of association information related to SNPs and phenotypes. Furthermore, CMPMI incorporates the interconnectedness of SNP combinations, making it symmetric in terms of describing associations.

3.6. An Improved $\text{ACO}\ast$

Given that the basic $\text{ACO}\ast$ exhibits low convergence speed and faces challenges with local minima problems [28, 48–50], we developed an improved $\text{ACO}\ast$.

To avoid getting trapped in local optima, it is crucial to expand the search space for ants. Based on the original path selection strategy, incorporating suitable random strategies can guide ants out of cyclic paths, thereby providing them with a more diverse set of path selections [51]. The corresponding formulas for path selection can be written as $$\begin{array}{}\text{(17)}& p_k^{i}t=\begin{array}{c}R\hspace{1em}q\le q_0\hfill \\ S\hspace{1em}q>q_0\hfill \end{array},\text{(18)}& S=\begin{array}{c}1\hspace{1em}i=\mathrm{rand}{M}_{k}t\hfill \\ 0\hspace{1em}\text{otherwise}\hfill \end{array},\end{array}$$where $p_k^{i}t$ is the probability that ant $k$ selects SNP $i$ in iteration $t$, $R$ is the original path selection strategy, $q$ is a randomly generated value satisfying a uniform distribution, and $q_0$ is the user specified threshold that is set to the reciprocal of the number of iterations.

For pheromone updating, the original updating strategy is adopted. Thus, $\Delta {\tau }_{\text{ij}}^{k}t$ can be written as $$\begin{array}{}\text{(19)}& \Delta {\tau }_i^{k}t=\begin{array}{c}CMPMIS\hspace{1em}k\in {M^{\ast }}_{i}t\hfill \\ 0\hspace{1em}otherwise\hfill \end{array},\end{array}$$where $\Delta {\tau }_i^{k}t$ is the pheromone variation of SNP $i$ selected by ant $k$ at iteration $t$, ${M^{\ast }}_{i}t$ represents the set of ants that select SNP $i$ at iteration $t$, $\text{CMPMI}S$ represents the CMPMI value of SNP combination $S$.

The memory-based strategy can retain superior solutions generated in each iteration, enhancing the overall convergence of the algorithm [51, 52]. Specifically, for each iteration, solutions captured by ants are sorted in descending order based on their CMPMI values. Subsequently, a turning point can be determined. $$\begin{array}{}\text{(20)}& f=\arg \underset{g=3}{\overset{m}{\max }}\begin{array}{c}\text{CMPMI}{S}_g-\text{CMPMI}{S}_{g-1}\\ -\text{CMPMI}{S}_{g-1}-\text{CMPMI}{S}_{g-2}\end{array},\end{array}$$where $\text{CMPMI}{S}_g$ is the CMPMI value of the SNP combination $S_g$, $g$ represents the ant. In each iteration, SNP combinations before the turning point are regarded as candidate solutions, and their corresponding fitness values are stored.

To further expedite convergence, a filtering operation based on the memory strategy is incorporated into ACOCMPMI. Within the candidate solution set obtained from each iteration, $\min \text{CMPMI}$ is utilized as the filter criterion. For subsequent iterations, those SNP combinations with CMPMI values greater than the filter criterion are retained and stored in the candidate solution set.

4. Results and Discussion

4.1. Evaluation Metrics

In the experiments, three evaluation metrics, including detection power, $F$-measure, and running time, are employed to assess the performance of compared methods.

Detection power is a widely used and effective metric for assessing the performance of methods for detecting epistatic interactions [39] and is defined as $$\begin{array}{}\text{(21)}& \text{power}=\frac{D_T}{D},\end{array}$$where $D_T$ is the number of datasets that epistatic interaction models in them are successfully detected and $D$ is the total number of datasets. Besides, the $F$-measure is defined as $$\begin{array}{}\text{(22)}& \text{Precision}=\frac{\text{TP}}{\text{TP}+\text{FP}},\text{(23)}& \text{Recall}=\frac{\text{TP}}{\text{TP}+\text{FN}},\text{(24)}& F‐\text{measure}=\frac{2\text{precision}\times \text{recall}}{\text{precision}+\text{recall}}=\frac{2\text{TP}}{2\text{TP}+\text{FP}+\text{FN}},\end{array}$$where true positives (TPs) represent that the detected SNP combinations are truly associated with the phenotype, false positives (FPs) represent that the detected SNP combinations are not associated with the phenotype, and false negatives (FNs) represent that the undetected SNP combinations are indeed associated with the phenotype.

4.2. Simulation Datasets

There are 11 epistatic interaction models to evaluate the performance of compared methods, where Models 1–8 are models displaying marginal effects (DMEs), and Models 9–11 are models displaying no marginal effects (DNMEs). Table 1 lists details of these models, in which MAF represents minor allele frequency, $\text{AA}$ is the homozygous common genotype, $\text{Aa}$ is the heterozygous genotype, and $\text{aa}$ is the homozygous minor genotype [49, 53]. Using these models, the simulator EpiSIM was applied to generate datasets of different scales [54]. For small-scale datasets, each model was used to generate 100 datasets, in which the sample number is 4000 and the SNP number is 100. For large-scale datasets, each model was used to generate 50 datasets, in which the sample number is 4000 and the SNP number is 1000.

Table 1

Details of epistatic interaction models.

Models	MAF(a)	MAF(b)	AABB	AABb	AAbb	AaBB	AaBb	Aabb	aaBB	aaBb	aabb
Model 1	0.2	0.2	0.087	0.087	0.087	0.087	0.146	0.190	0.087	0.190	0.247
Model 2	0.5	0.5	0.009	0.009	0.009	0.013	0.006	0.006	0.013	0.006	0.006
Model 3	0.5	0.5	0.092	0.092	0.092	0.092	0.319	0.319	0.092	0.319	0.319
Model 4	0.2	0.2	0.084	0.084	0.084	0.084	0.210	0.210	0.084	0.210	0.210
Model 5	0.5	0.5	0.052	0.052	0.052	0.052	0.137	0.137	0.052	0.137	0.137
Model 6	0.5	0.5	0.072	0.164	0.164	0.164	0.072	0.072	0.164	0.072	0.072
Model 7	0.5	0.5	0.067	0.155	0.155	0.155	0.067	0.067	0.155	0.067	0.067
Model 8	0.3	0.3	0.486	0.960	0.538	0.947	0.004	0.811	0.640	0.606	0.909
Model 9	0.2	0.5	0.103	0.063	0.124	0.098	0.086	0.069	0.021	0.147	0.059
Model 10	0.5	0.5	0.000	0.000	0.000	0.000	0.050	0.000	0.100	0.000	0.000
Model 11	0.3	0.3	0.000	0.020	0.000	0.020	0.000	0.020	0.000	0.020	0.000

4.3. Results on Simulation Datasets

For small-scale datasets, the ant number and the iteration number are set to 200 and 70, respectively, while for large-scale datasets, the ant number and the iteration number are set to 2000 and 100, respectively. The detection power of ACOCMPMI with different iteration numbers is precomputed for all models, and those iteration numbers close to the optimal convergence point are selected as the iteration parameters, as illustrated in Figure 2.

[figure(s) omitted; refer to PDF]

For small-scale datasets, detection power and $F$-measure of compared methods are presented in Figure 3. In terms of detection power, most methods perform well and detect almost all epistatic interactions in various datasets. Specifically, ACOCMPMI demonstrates high and stable detection power in DMEs, comparable to FDHE-IW and MACOED. Notably, FDHE-IW is a method specifically designed for detecting DMEs [55]. ACOCMPMI exhibits lower detection power than MACOED in small-scale DNME datasets, contrasting with its superior performance in large-scale datasets. Although AntEpiSeeker performs effectively in Model 4–7 datasets, it fails to detect epistatic interactions in Model 2, 8, and 11 datasets, implying that AntEpiSeeker may be inconsistent and exhibit model preference. SIPSO shows similar performance to AntEpiSeeker but with greater stability. However, SIPSO struggles to adapt to DNMEs. In terms of $F$-measure, ACOCMPMI significantly outperforms most compared methods in DMEs, though its performance is inferior to MACOED in DNMEs.

[figure(s) omitted; refer to PDF]

For large-scale datasets, detection power and F-measure of compared methods are presented in Figure 4. In terms of detection power, ACOCMPMI outperforms all compared methods in almost all datasets except Model 1–2 datasets. Performance of ACOCMPMI ranks second only to SIPSO in Model 1 datasets and to FDHE-IW in Model 2 datasets, respectively, further demonstrating the stability of its detection capability. AntEpiSeeker and MACOED show detection power ranging from 0.1 to 0.5 in most models, which is significantly lower than the detection power of ACOCMPMI. SIPSO performs effectively in datasets of Models 1, 5, and 9–11, but fails to identify over 60% of epistatic interactions in other datasets. epiACO and FDHE-IW exhibit detection power comparable to ACOCMPMI. epiACO performs well in most models since both it and ACOCMPMI use the $\text{ACO}\ast$ and the information theory-based quantification measure. In terms of $F$-measure, ACOCMPMI has higher values than those of compared methods in almost all datasets except Model 1–2 datasets. Although epiACO and FDHE-IW are as effective as ACOCMPMI in identifying epistatic interactions in most models, their $F$-measure values vary widely among models, implying that both have weaker stability than ACOCMPMI. SIPSO, AntEpiSeeker, and MACOED generally have low $F$-measure values in most models, which is consistent with their performance in detection power.

[figure(s) omitted; refer to PDF]

Running times of compared methods in different datasets are shown in Figure 5. It is seen that in small-scale datasets, ACOCMPMI has similar running times to those of both epiACO and SIPSO in various models. Running times of AntEpiSeeker in all models are relatively stable, though it takes more time than ACOCMPMI, epiACO, and SIPSO. MACOED shows significantly varying running times across models, implying that it is sensitive to model type. FDHE-IW requires unacceptable running times in all models. For large-scale datasets, in DMEs and DNMEs, ACOCMPMI has a clear advantage in terms of running time. Unlike FDHE-IW, which has the worst running times in small-scale datasets, MACOED becomes the most time-consuming method in large-scale datasets. Though SIPSO and epiACO have acceptable running times, their detection power is low.

[figure(s) omitted; refer to PDF]

To demonstrate that the improved $\text{ACO}\ast$ in ACOCMPMI is effective for searching epistatic interactions, $\text{ACO}\ast$ is compared with AVOA in small-scale datasets, using CMPMI as their fitness function, in terms of detection power, $F$-measure, and running time, as shown in Figure 6. It is seen that even when facing the recently developed meta-heuristic algorithm AVOA, $\text{ACO}\ast$ still has an advantage in search performance. In general, the random strategy and memory-filter strategy incorporated into the basic $\text{ACO}\ast$ improve its detection capability without increasing running time.

[figure(s) omitted; refer to PDF]

4.4. Case Study

ACOCMPMI is applied to a real AMD dataset to detect two-order epistatic interactions. The AMD dataset contains 103,611 SNPs with 50 controls and 96 cases and has become a widely used benchmark dataset [39, 53]. ACOCMPMI runs four times on this AMD dataset, using ants and iterations as (10,000, 500), (10,000, 1000), (20,000, 250), and (20,000, 1000), respectively, to capture more epistatic interactions. Table 2 lists the Top 15 detected epistatic interactions associated with AMD.

Table 2

Top 15 detected epistatic interactions associated with AMD.

SNP 1			SNP 2			Fitness value	$\mathit{p}$ value	Times
Name	Gene	Chr	Name	Gene	Chr
rs3775652	INPP4B	4	rs380390	CFH	1	123.48	0.0175	3
rs3775652	INPP4B	4	rs725518	RRM1	11	121.99	0.0149	3
rs380390	CFH	1	rs725518	RRM1	11	121.05	0.0039	3
rs380390	CFH	1	rs54816	RRM1	11	120.19	0.0064	3
rs3775652	INPP4B	4	rs54816	RRM1	11	119.90	0.0081	3
rs7863587	/	9	rs380390	CFH	1	119.63	0.0415	2
rs4772270	PCCA	13	rs380390	CFH	1	118.82	0.0265	1
rs6480996	/	10	rs380390	CFH	1	118.11	0.0052	1
rs2019727	CFH	1	rs380390	CFH	1	118.05	0.0082	1
rs380390	CFH	1	rs365299	/	1	117.10	0.0190	1
rs3775652	INPP4B	4	rs4772270	PCCA	13	115.98	0.0458	1
rs7863587	/	9	rs3775652	INPP4B	4	115.49	0.0394	1
rs3775652	INPP4B	4	rs3775650	INPP4B	4	113.31	0.0234	1
rs3775650	INPP4B	4	rs4772270	PCCA	13	112.69	0.0046	1
rs7863587	/	9	rs725518	RRM1	11	111.21	0.0285	1

rs380390 is a G/A/T/C single-nucleotide variation in the CFH gene on human chromosome 1, and rs2019727, also located in CFH, is considered to be significantly associated with AMD in several studies [56–61]. rs3775652 is a C/T single-nucleotide variation located in the INPP4B gene on chromosome 4, and rs725518 is an A/G single-nucleotide variation in the RRM1 gene on chromosome 11, both of which have been detected as AMD-related SNPs [62, 63]. rs4772270 is a G/A/T/C single-nucleotide variation in the PCCA gene on chromosome 13, which has also been reported to be associated with AMD [55, 62, 63]. More recently, rs7863587 was reported to be highly associated with AMD [64]. Although further experiments and clinical studies are needed to confirm real epistatic interactions with AMD, we hope that these findings of ACOCMPMI can provide some clues for the pathological study of AMD.

5. Conclusions and Future Works

Epistatic interaction detection plays a pivotal role in understanding the genetic mechanisms underlying complex diseases. The effectiveness of epistatic interaction detection methods primarily depends on their interaction quantification measures and search strategies. Therefore, both are significant challenges for epistatic interaction detection. In this study, ACOCMPMI, a two-stage $\text{ACO}\ast$ based on composite MPMI is proposed for detecting epistatic interactions. In the first stage, CMPMI is introduced to quantify epistatic interactions, and an improved $\text{ACO}\ast$, incorporating filter and memory strategies, is employed to search for epistatic interactions. In the second stage, an exhaustive strategy and a BN score, that is, K2 score, are adopted to further identify epistatic interactions within the candidate SNP set obtained from the first stage. ACOCMPMI is compared with five state-of-the-art methods, including epiACO, FDHE-IW, AntEpiSeeker, SIPSO, and MACOED, using simulation data based on 11 epistatic interaction models. Furthermore, ACOCMPMI is applied to detect epistatic interactions in a real dataset related to AMD. The experimental results show that ACOCMPMI is an alternative method for epistatic interaction detection. The time complexity of ACOCMPMI is $O\text{NT}+\text{n}{\text{m}}^2$, where $N$, $T$, $n$, and $m$ are numbers of ants, iterations, SNPs, and samples, respectively.

However, there are still several limitations in ACOCMPMI, which inspire us to continue working. First, how to adjust parameter settings to adapt to different scales of input SNP datasets should be further discussed. Second, the practical applicability and scalability of ACOCMPMI require a more detailed analysis. Although some of the identified SNPs have been validated, it remains unclear whether their two-order combinations are indeed causal factors of AMD. Furthermore, the current version of ACOCMPMI focuses on capturing two-order epistatic interactions. In reality, complex diseases are often caused by epistatic interactions with different orders, especially higher orders. Therefore, its future version should be developed to detect higher order epistatic interactions.

Author Contributions

Yan Sun and Jing Wang contributed equally to this work.

Funding

This work was supported by the National Natural Science Foundation of China (62472250, 62473179, and 62172254).

References

[1] I. Dinu, S. Mahasirimongkol, Q. Liu, H. Yanai, N. Sharaf Eldin, E. Kreiter, X. Wu, S. Jabbari, K. Tokunaga, Y. Yasui, "SNP-SNP Interactions Discovered by Logic Regression Explain Crohn's Disease Genetics," PLoS One, vol. 7 no. 10,DOI: 10.1371/journal.pone.0043035, 2012.

[2] M. D. Ritchie, "Using Biological Knowledge to Uncover the Mystery in the Search for Epistasis in Genome-Wide Association Studies," Annals of Human Genetics, vol. 75 no. 1, pp. 172-182, DOI: 10.1111/j.1469-1809.2010.00630.x, 2011.

[3] T. Zhao, Y. Hu, T. Zang, Y. Wang, "Integrate GWAS, eQTL, and mQTL Data to Identify Alzheimer’s Disease-Related Genes," Frontiers in Genetics, vol. 10,DOI: 10.3389/fgene.2019.01021, 2019.

[4] Z. R. Sailer, M. J. Harms, "Detecting High-Order Epistasis in Nonlinear Genotype-Phenotype Maps," Genetics, vol. 205 no. 3, pp. 1079-1088, DOI: 10.1534/genetics.116.195214, 2017.

[5] Y. Jiang, J. C. Reif, "Efficient Algorithms for Calculating Epistatic Genomic Relationship Matrices," Genetics, vol. 216 no. 3, pp. 651-669, DOI: 10.1534/genetics.120.303459, 2020.

[6] A. J. Morrison, D. R. Wonderlick, M. J. Harms, "Ensemble Epistasis: Thermodynamic Origins of Nonadditivity Between Mutations," Genetics, vol. 219 no. 1,DOI: 10.1093/genetics/iyab105, 2021.

[7] L. Wienbrandt, J. C. Kässens, M. Hübenthal, D. Ellinghaus, "1000× Faster Than PLINK: Combined FPGA and GPU Accelerators for Logistic Regression-Based Detection of Epistasis," Journal of Computational Science, vol. 30, pp. 183-193, DOI: 10.1016/j.jocs.2018.12.013, 2019.

[8] X. Wan, C. Yang, Q. Yang, H. Xue, X. Fan, N. L. Tang, W. Yu, "BOOST: A Fast Approach to Detecting Gene-Gene Interactions in Genome-Wide Case-Control Studies," The American Journal of Human Genetics, vol. 87 no. 3, pp. 325-340, DOI: 10.1016/j.ajhg.2010.07.021, 2010.

[9] M. Xie, J. Li, T. Jiang, "Detecting Genome-Wide Epistases Based on the Clustering of Relatively Frequent Items," Bioinformatics, vol. 28 no. 1,DOI: 10.1093/bioinformatics/btr603, 2012.

[10] M. Pitsillou, K. Fokianos, dCovTS: Distance Covariance/Correlation for Time Series, 2016.

[11] Q. Ding, J. Shang, Y. Sun, G. Liu, F. Li, X. Yuan, J.-X. Liu, "NIPMI: A Network Method Based on Interaction Part Mutual Information to Detect Characteristic Genes From Integrated Data on Multi-Cancers," IEEE Access, vol. 7, pp. 135845-135854, DOI: 10.1109/ACCESS.2019.2941520, 2019.

[12] X. Zhang, X.-M. Zhao, K. He, L. Lu, Y. Cao, J. Liu, J.-K. Hao, Z.-P. Liu, L. Chen, "Inferring Gene Regulatory Networks From Gene Expression Data by Path Consistency Algorithm Based on Conditional Mutual Information," Bioinformatics, vol. 28 no. 1, pp. 98-104, DOI: 10.1093/bioinformatics/btr626, 2012.

[13] J. A. Kontio, M. J. Rinta-Aho, M. J. Sillanpää, "Estimating Linear and Nonlinear Gene Coexpression Networks by Semiparametric Neighborhood Selection," Genetics, vol. 215 no. 3, pp. 597-607, DOI: 10.1534/genetics.120.303186, 2020.

[14] D. G. Hernández Lahme, I. Samengo, "Estimating the Mutual Information Between Two Discrete, Asymmetric Variables With Limited Samples," Entropy, vol. 21 no. 6,DOI: 10.3390/e21060623, 2019.

[15] X. Cao, G. Yu, J. Liu, L. Jia, J. Wang, "Clustermi: Detecting High-Order SNP Interactions Based on Clustering and Mutual Information," International Journal of Molecular Sciences, vol. 19 no. 8,DOI: 10.3390/ijms19082267, 2018.

[16] J. Zhao, Y. Zhou, X. Zhang, L. Chen, "Part Mutual Information for Quantifying Direct Associations in Networks," Proceedings of the National Academy of Sciences, vol. 113 no. 18, pp. 5130-5135, DOI: 10.1073/pnas.1522586113, 2016.

[17] D. N. Reshef, Y. A. Reshef, H. K. Finucane, S. R. Grossman, G. McVean, P. J. Turnbaugh, E. S. Lander, M. Mitzenmacher, P. C. Sabeti, "Detecting Novel Associations in Large Data Sets," Science, vol. 334 no. 6062, pp. 1518-1524, DOI: 10.1126/science.1205438, 2011.

[18] X. Zhang, J. Zhao, J.-K. Hao, X.-M. Zhao, L. Chen, "Conditional Mutual Inclusive Information Enables Accurate Quantification of Associations in Gene Regulatory Networks," Nucleic Acids Research, vol. 43 no. 5, pp. e31-e31, DOI: 10.1093/nar/gku1315, 2015.

[19] J. Shi, J. Zhao, X. Liu, L. Chen, T. Li, "Quantifying Direct Dependencies in Biological Networks by Multiscale Association Analysis," IEEE/ACM Transactions on Computational Biology and Bioinformatics, vol. 17 no. 2, pp. 449-458, DOI: 10.1109/TCBB.2018.2846648, 2020.

[20] J. Shang, J. Wang, Y. Sun, F. Li, J.-X. Liu, H. Zhang, "Multiscale Part Mutual Information for Quantifying Nonlinear Direct Associations in Networks," Bioinformatics, vol. 37 no. 18, pp. 2920-2929, DOI: 10.1093/bioinformatics/btab182, 2021.

[21] A. Upton, O. Trelles, J. A. Cornejo-García, J. R. Perkins, "Review: High-Performance Computing to Detect Epistasis in Genome Scale Data Sets," Briefings in Bioinformatics, vol. 17 no. 3, pp. 368-379, DOI: 10.1093/bib/bbv058, 2016.

[22] Y. Zhang, J. S. Liu, "Bayesian Inference of Epistatic Interactions in Case-Control Studies," Nature Genetics, vol. 39 no. 9, pp. 1167-1173, 2007.

[23] M. Abedi, F. S. Gharehchopogh, "An Improved Opposition Based Learning Firefly Algorithm With Dragonfly Algorithm for Solving Continuous Optimization Problems," Intelligent Data Analysis, vol. 24 no. 2, pp. 309-338, DOI: 10.3233/IDA-194485, 2020.

[24] F. S. Gharehchopogh, "Advances in Tree Seed Algorithm: A Comprehensive Survey," Archives of Computational Methods in Engineering, vol. 29 no. 5, pp. 3281-3304, DOI: 10.1007/s11831-021-09698-0, 2022.

[25] F. S. Gharehchopogh, "An Improved Tunicate Swarm Algorithm With Best-Random Mutation Strategy for Global Optimization Problems," Journal of Bionic Engineering, vol. 19 no. 4, pp. 1177-1202, DOI: 10.1007/s42235-022-00185-1, 2022.

[26] O. Maciel, E. Cuevas, M. A. Navarro, D. Zaldívar, S. Hinojosa, "Side-Blotched Lizard Algorithm: A Polymorphic Population Approach," Applied Soft Computing, vol. 88,DOI: 10.1016/j.asoc.2019.106039, 2020.

[27] B. Abdollahzadeh, F. S. Gharehchopogh, S. Mirjalili, "African Vultures Optimization Algorithm: A New Nature-Inspired Metaheuristic Algorithm for Global Optimization Problems," Computers & Industrial Engineering, vol. 158,DOI: 10.1016/j.cie.2021.107408, 2021.

[28] J. Shang, X. Wang, X. Wu, Y. Sun, Q. Ding, J.-X. Liu, H. Zhang, "A Review of Ant Colony Optimization Based Methods for Detecting Epistatic Interactions," IEEE Access, vol. 7, pp. 13497-13509, DOI: 10.1109/ACCESS.2019.2894676, 2019.

[29] H. Mohammadzadeh, F. S. Gharehchopogh, "Feature Selection With Binary Symbiotic Organisms Search Algorithm for Email Spam Detection," International Journal of Information Technology & Decision Making, vol. 20 no. 1, pp. 469-515, DOI: 10.1142/S0219622020500546, 2021.

[30] S. Ghafori, F. S. Gharehchopogh, "Advances in Spotted Hyena Optimizer: A Comprehensive Survey," Archives of Computational Methods in Engineering, vol. 29 no. 3, pp. 1569-1590, DOI: 10.1007/s11831-021-09624-4, 2022.

[31] D. Zaldívar, B. Morales, A. Rodríguez, A. Valdivia-G, E. Cuevas, M. Pérez-Cisneros, "A Novel Bio-Inspired Optimization Model Based on Yellow Saddle Goatfish Behavior," Biosystems, vol. 174,DOI: 10.1016/j.biosystems.2018.09.007, 2018.

[32] A. Rodríguez, O. Camarena, E. Cuevas, I. Aranguren, A. Valdivia-G, B. Morales-Castañeda, D. Zaldívar, M. Pérez-Cisneros, "Group-Based Synchronous-Asynchronous Grey Wolf Optimizer," Applied Mathematical Modelling, vol. 93, pp. 226-243, DOI: 10.1016/j.apm.2020.12.016, 2021.

[33] S. Tuo, C. Li, F. Liu, Y. Zhu, T. Chen, Z. Feng, H. Liu, A. Li, "A Novel Multitasking Ant Colony Optimization Method for Detecting Multiorder SNP Interactions," Interdisciplinary Sciences: Computational LIfe Sciences, vol. 14 no. 4, pp. 814-832, DOI: 10.1007/s12539-022-00530-2, 2022.

[34] M. D. Ritchie, L. W. Hahn, N. Roodi, L. R. Bailey, W. D. Dupont, F. F. Parl, J. H. Moore, "Multifactor-Dimensionality Reduction Reveals High-Order Interactions Among Estrogen-Metabolism Genes in Sporadic Breast Cancer," The American Journal of Human Genetics, vol. 69 no. 1, pp. 138-147, 2001.

[35] T. Zheng, H. Wang, S.-H. Lo, "Backward Genotype-Trait Association (BGTA)-Based Dissection of Complex Traits in Case-Control Designs," Human Heredity, vol. 62 no. 4, pp. 196-212, 2006.

[36] C. Lippert, J. Listgarten, R. I. Davidson, J. Baxter, H. Poon, C. M. Kadie, D. Heckerman, "An Exhaustive Epistatic SNP Association Analysis on Expanded Wellcome Trust Data," Scientific Reports, vol. 3 no. 1,DOI: 10.1038/srep01099, 2013.

[37] X. Zhang, S. Huang, F. Zou, W. Wang, "TEAM: Efficient Two-Locus Epistasis Tests in Human Genome-Wide Association Study," Bioinformatics, vol. 26 no. 12, pp. i217-i227, DOI: 10.1093/bioinformatics/btq186, 2010.

[38] W. Tang, X. Wu, R. Jiang, Y. Li, "Epistatic Module Detection for Case-Control Studies: A Bayesian Model With a Gibbs Sampling Strategy," PLoS Genetics, vol. 5 no. 5,DOI: 10.1371/journal.pgen.1000464, 2009.

[39] J. Shang, Y. Sun, J.-X. Liu, J. Xia, J. Zhang, C.-H. Zheng, "CINOEDV: A Co-Information Based Method for Detecting and Visualizing N-Order Epistatic Interactions," BMC Bioinformatics, vol. 17 no. 1,DOI: 10.1186/s12859-016-1076-8, 2016.

[40] Y. Wang, X. Liu, K. Robbins, R. Rekaya, "AntEpiSeeker: Detecting Epistatic Interactions for Case-Control Studies Using a Two-Stage Ant Colony Optimization Algorithm," BMC Research Notes, vol. 3, 2010.

[41] P.-J. Jing, H.-B. Shen, "MACOED: A Multi-Objective Ant Colony Optimization Algorithm for SNP Epistasis Detection in Genome-Wide Association Studies," Bioinformatics, vol. 31 no. 5, pp. 634-641, DOI: 10.1093/bioinformatics/btu702, 2015.

[42] Y. Sun, J. Shang, J. Liu, S. Li, "An Improved Ant Colony Optimization Algorithm for the Detection of SNP-SNP Interactions," Intelligent Computing Methodologies: 12th International Conference, ICIC 2016, Lanzhou, China, August 2-5, 2016, Proceedings, Part III 12, pp. 21-32, 2016.

[43] S. Tuo, C. Li, F. Liu, A. Li, L. He, Z. W. Geem, J. Shang, H. Liu, Y. Zhu, Z. Feng, T. R. Chen, "MTHSA-DHEI: Multitasking Harmony Search Algorithm for Detecting High-Order SNP Epistatic Interactions," Complex & Intelligent Systems, vol. 9 no. 1, pp. 637-658, DOI: 10.1007/s40747-022-00813-7, 2023.

[44] S. Tuo, J. Jiang, "A Novel Detection Method for High-Order SNP Epistatic Interactions Based on Explicit-Encoding-Based Multitasking Harmony Search," Interdisciplinary Sciences: Computational LIfe Sciences, vol. 16 no. 3, pp. 688-711, DOI: 10.1007/s12539-024-00621-2, 2024.

[45] J. Shang, J. Zhang, X. Lei, Y. Zhang, B. Chen, "Incorporating Heuristic Information Into Ant Colony Optimization for Epistasis Detection," Genes & Genomics, vol. 34 no. 3, pp. 321-327, DOI: 10.1007/s13258-012-0003-2, 2012.

[46] C. Bielza, P. Larranaga, "Discrete Bayesian Network Classifiers," ACM Computing Surveys (CSUR), vol. 47 no. 1,DOI: 10.1145/2576868, 2014.

[47] G. F. Cooper, E. Herskovits, "A Bayesian Method for the Induction of Probabilistic Networks From Data," Machine Learning, vol. 9 no. 4, pp. 309-347, DOI: 10.1007/BF00994110, 1992.

[48] S. Chowdhury, M. Marufuzzaman, H. Tunc, L. Bian, W. Bullington, "A Modified Ant Colony Optimization Algorithm to Solve a Dynamic Traveling Salesman Problem: A Case Study With Drones for Wildlife Surveillance," Journal of Computational Design and Engineering, vol. 6 no. 3, pp. 368-386, DOI: 10.1016/j.jcde.2018.10.004, 2019.

[49] F. Dahan, K. El Hindi, H. Mathkour, H. AlSalman, "Dynamic Flying Ant Colony Optimization (DFACO) for Solving the Traveling Salesman Problem," Sensors, vol. 19 no. 8,DOI: 10.3390/s19081837, 2019.

[50] B. Guan, Y. Zhao, "Self-Adjusting Ant Colony Optimization Based on Information Entropy for Detecting Epistatic Interactions," Genes, vol. 10 no. 2,DOI: 10.3390/genes10020114, 2019.

[51] Y. Sun, J. Shang, J.-X. Liu, S. Li, C.-H. Zheng, "epiACO-a Method for Identifying Epistasis Based on Ant Colony Optimization Algorithm," Biodata Mining, vol. 10, 2017.

[52] Y. Sun, X. Wang, J. Shang, J.-X. Liu, C.-H. Zheng, X. Lei, "Introducing Heuristic Information Into Ant Colony Optimization Algorithm for Identifying Epistasis," IEEE/ACM Transactions on Computational Biology and Bioinformatics, vol. 17 no. 4, pp. 1253-1261, DOI: 10.1109/TCBB.2018.2879673, 2020.

[53] P. Mitchell, G. Liew, B. Gopinath, T. Y. Wong, "Age-Related Macular Degeneration," Lancet, vol. 392 no. 10153, pp. 1147-1159, DOI: 10.1016/S0140-6736(18)31550-2, 2018.

[54] J. Shang, J. Zhang, X. Lei, W. Zhao, Y. Dong, "EpiSIM: Simulation of Multiple Epistasis, Linkage Disequilibrium Patterns and Haplotype Blocks for Genome-Wide Interaction Analysis," Genes & Genomics, vol. 35 no. 3, pp. 305-316, DOI: 10.1007/s13258-013-0081-9, 2013.

[55] S. Tuo, "FDHE-IW: A Fast Approach for Detecting High-Order Epistasis in Genome-Wide Case-Control Studies," Genes, vol. 9 no. 9,DOI: 10.3390/genes9090435, 2018.

[56] W.-Y. Lin, W.-C. Lee, "Incorporating Prior Knowledge to Facilitate Discoveries in a Genome-Wide Association Study on Age-Related Macular Degeneration," BMC Research Notes, vol. 3 no. 1, pp. 25-26, DOI: 10.1186/1756-0500-3-26, 2010.

[57] J. Tuo, R. J. Ross, G. F. Reed, Q. Yan, J. J. Wang, C. M. Bojanowski, E. Y. Chew, X. Feng, T. W. Olsen, F. L. Ferris, "The HtrA1 Promoter Polymorphism, Smoking, and Age-Related Macular Degeneration in Multiple Case-Control Samples," Ophthalmology, vol. 115 no. 11, pp. 1891-1898, DOI: 10.1016/j.ophtha.2008.05.021, 2008.

[58] J. Tian, W. Yu, X. Qin, K. Fang, Q. Chen, J. Hou, J. Li, D. Chen, Y. Hu, X. Li, "Association of Genetic Polymorphisms and Age-Related Macular Degeneration in Chinese Population," Investigative Ophthalmology & Visual Science, vol. 53 no. 7, pp. 4262-4269, DOI: 10.1167/iovs.11-8542, 2012.

[59] P. Gili, L. Lloreda Martín, J.-C. Martín-Rodrigo, N. Kim-Yeon, L. Modamio-Gardeta, J. L. Fernández-García, A. B. Rebolledo-Poves, E. Gómez-Blazquez, R. Pazos-Rodriguez, E. Pérez-Fernández, M. Velasco, "Gene Polymorphisms Associated With an Increased Risk of Exudative Age-Related Macular Degeneration in a Spanish Population," European Journal of Ophthalmology, vol. 32 no. 1, pp. 651-657, DOI: 10.1177/11206721211002698, 2022.

[60] B. Feigl, C. P. Morris, B. Brown, A. J. Zele, "Relationship Among CFH and ARMS2 Genotypes, Macular Pigment Optical Density, and Neuroretinal Function in Persons Without Age-Related Macular Degeneration," Archives of Ophthalmology, vol. 130 no. 11, pp. 1402-1409, DOI: 10.1001/archophthalmol.2012.1940, 2012.

[61] M. Budzinskaia, T. Pogoda, É. Generozov, E. Chikun, I. Shchegoleva, É. Kazarian, N. Galoian, "Influence of Genetic Mutations on Clinical Presentation of Subretinal Neovascularization. Report 1: The Impact of CFH and Il-8 Genes Polymorphism," Vestnik Oftalmologii, vol. 127 no. 4, 2011.

[62] S. Tuo, J. Zhang, X. Yuan, Y. Zhang, Z. Liu, "FHSA-SED: Two-Locus Model Detection for Genome-Wide Association Study With Harmony Search Algorithm," PLoS One, vol. 11 no. 3,DOI: 10.1371/journal.pone.0150669, 2016.

[63] Y. Guo, Z. Zhong, C. Yang, J. Hu, Y. Jiang, Z. Liang, H. Gao, J. Liu, "Epi-GTBN: An Approach of Epistasis Mining Based on Genetic Tabu Algorithm and Bayesian Network," BMC Bioinformatics, vol. 20 no. 1,DOI: 10.1186/s12859-019-3022-z, 2019.

[64] R. Jiang, W. Tang, X. Wu, W. Fu, "A Random Forest Approach to the Detection of Epistatic Interactions in Case-Control Studies," BMC Bioinformatics, vol. 10 no. S1,DOI: 10.1186/1471-2105-10-S1-S65, 2009.