Multiobjective differential evolution-based multifactor dimensionality reduction for detecting gene–gene interactions

Epistasis within disease-related genes (gene–gene interactions) was determined through contingency table measures based on multifactor dimensionality reduction (MDR) using single-nucleotide polymorphisms (SNPs). Most MDR-based methods use the single contingency table measure to detect gene–gene interactions; however, some gene–gene interactions may require identification through multiple contingency table measures. In this study, a multiobjective differential evolution method (called MODEMDR) was proposed to merge the various contingency table measures based on MDR to detect significant gene–gene interactions. Two contingency table measures, namely the correct classification rate and normalized mutual information, were selected to design the fitness functions in MODEMDR. The characteristics of multiobjective optimization enable MODEMDR to use multiple measures to efficiently and synchronously detect significant gene–gene interactions within a reasonable time frame. Epistatic models with and without marginal effects under various parameter settings (heritability and minor allele frequencies) were used to assess existing methods by comparing the detection success rates of gene–gene interactions. The results of the simulation datasets show that MODEMDR is superior to existing methods. Moreover, a large dataset obtained from the Wellcome Trust Case Control Consortium was used to assess MODEMDR. MODEMDR exhibited efficiency in identifying significant gene–gene interactions in genome-wide association studies.

problems 20 , in which n (n > 1) objectives are considered to synchronously search for optimal solutions 21 ; for example, maximization objectives can be formulated as maximize(f 1 (X), …, f i (X)), where X ∈  X , i is the number of objectives,  X is the feasible solution set, and f(X) is an objective function. In maximization problems, solution X 1 dominates solution X 2 if f j (X 1 ) > f j (X 2 ) for all indices j ∈ (1, …, n). Pareto optimal solution sets (Pareto sets) represent a powerful technique for collecting good solutions not dominated by one another. These good solutions are the results of MODE.
Several contingency table measures, such as chi-square, likelihood ratio, normalized mutual information (NMI), and et al., have been applied to score model quality in MDR 22,23 , and these measures can be regarded as various objectives in MODE. Currently, MDR-based methods focus only on a single measure to determine genegene interactions. Various simulation dataset types have been adopted to evaluate which contingency table measures can significantly improve MDR performance 22 , revealing that MDR performance could be measured based on the correct classification rate (CCR) 10 or NMI 22 . However, no optimal measure for determining gene-gene interactions involving various dataset types has yet been found. Each measure may fit specific dataset types; however, deriving data distributions from real datasets is difficult, especially for complex diseases. Therefore, developing a method that can synchronously consider multiple measures to detect gene-gene interactions is essential.
In this study, a multiobjective DE (hereafter MODEMDR) was proposed to merge various contingency table measures based on MDR and detect significant gene-gene interactions. Two objectives involving the aforementioned two measures of CCR and NMI were selected for MODEMDR. Several epistatic models with and without marginal effects and with various parameter settings (heritability (h 2 ) and minor allele frequencies (MAF)) were selected to generate high-dimensional simulation datasets. In addition, a large real dataset was obtained from the Wellcome Trust Case Control Consortium (WTCCC) 24 . The results of the simulation and real datasets indicated that MODEMDR can effectively detect gene-gene interactions.

Results
Simulation data experiments. The goal of the simulation datasets was to successfully detect the specific two-locus SNP combination (target) in each artifact epistasis model. Epistatic models with and without marginal effects were simulated to compare the epistatic interaction identification ability of SNPRuler 11 , MDR 25 , single measure DE MDR (DEMDR), and MODEMDR.
Comparison between MODEMDR and existing methods on disease loci with marginal effects. The eight epistatic models with marginal effects were used to evaluate the performance of SNPRuler, MDR, DEMDR (CCR), and MODEMDR. Models 1-6 were obtained from Namkung et al. 23 and models 7 and 8 were obtained from Bush et al. 22 . These models reflect the strength of genetic effects and were proposed according to the interaction structure, MAF, and prevalence. The details of the multilocus penetrances of the eight models are shown in Table 1 in the supplementary file. The penetrances of the eight models were computed under the Hardy-Weinberg equilibrium (HWE) assumption for each SNP. In each model, 100 datasets were simulated under identical settings with uniform MAF of [0.05, 0.5). The detection success rate was computed as the proportion of the generated datasets, in which a target of epistatic interaction was detected. GAMETES software was used to simulate the simulation datasets 26 .
In the eight models, MDR, DEMDR, and MODEMDR outperformed SNPRuler in the large samples ( Fig. 1; 1,000 cases and 1,000 controls), in which MODEMDR outperformed MDR and DEMDR in models 7 and 8. Regarding the small samples (200 cases and 200 controls), SNPRuler, MDR, and DEMDR had difficulties identifying the specific two-locus SNP combinations in the epistatic models with marginal effects. Clearly in the small samples, MODEMDR outperformed MDR, DEMDR, and SNPRuler in the eight epistatic models with marginal effects. The generated datasets of eight epistatic models with marginal effects were used to compare DEMDR (CCR) (P), DEMDR (NMI) (N), and MODEMDR (two objectives merging CCR and NMI) (B). DEMDR (CCR) achieved higher detection success rates than DEMDR (NMI) in all epistatic models with marginal effects (Fig. 2). Moreover, MODEMDR outperformed DEMDR (CCR) and DEMDR (NMI), indicating that multiple contingency table measures are superior to single contingency table measures in detecting epistatic interactions with marginal effects. MODE effectively improves MDR with respect to performing evaluations to facilitate the identification of significant gene-gene interactions.
Comparison between MODEMDR and existing methods on disease loci without marginal effects. A total of 60 two-locus epistatic models were obtained from Wan et al. 11 and used to assess the performance of SNPRuler, MDR, DEMDR (CCR), and MODEMDR. These models are pure epistatic models (i.e., they have no marginal effects). The multilocus penetrances are shown in Supplementary Table 2. The parameter settings (h 2 and MAF) were selected to generate simulation data by using GAMETES software 26 . The h 2 controlled the phenotypic variation of the 60 models and ranged from 0.025 to 0.4. The MAF ranged from 0.2 to 0.4. For each epistatic model, the 100 datasets consisting of 1,000 SNPs, 200 cases, and 200 controls were generated. The detection success rate was calculated as the proportion of the 100 datasets in which the specific disease-associated two-locus SNP combination was detected.
In the 60 models, MODEMDR outperformed SNPRuler, MDR, and DEMDR in detecting epistatic interactions without marginal effects (Fig. 3). The results of Wilcoxon signed-rank testing (Table 1) showed that MODEMDR achieved the highest R + (number of victories), lowest R − (number of losses), and a p value of < 0.05, indicating that MODEMDR is significantly superior to the other methods. In the epistatic models with MAF = 0.2 or 0.4 and h 2 ≥ 0.2, all detection success rates of SNPRuler, MDR, DEMDR, and MODEMDR were ≥ 80%, which degraded as h 2 was decreased. When MAF = 0.2 and h 2 ≤ 0.05, DEMDR and MODEMDR achieved detection success rates of approximately 30% and 40%, respectively. By contrast, SNPRuler and MDR almost completely lost their detection abilities. MODEMDR achieved the highest detection success rates for all settings, especially h 2 ≤ 0.01 (Fig. 4). All the test results show that MODEMDR outperformed SNPRuler, MDR, and DEMDR in the epistatic models with no marginal effects.
The generated datasets of models 31-60 were used to compare the DEMDR (CCR) (D), DEMDR (NMI) (N), and MODEMDR (two objectives merging CCR and NMI) (O) in detecting epistatic interactions without marginal effects. Detection success rates were calculated as the proportion of the 100 datasets in which the specific  Under each setting, the detection success rate was calculated as the proportion of 100 datasets in which a specific disease-associated epistatic interaction was detected. Each dataset contained 1,000 SNPs. The gray bars represent the detection success rate for 1,000 cases and 1,000 controls. The black bars represent the detection success rate for 200 cases and 200 controls. The absence of bars indicates a detection success rate of zero.
disease-associated two-locus SNP combination was identified. DEMDR (CCR) achieved a higher detection success rate than DEMDR (NMI) in all epistatic models without marginal effects (Fig. 5) Results of WTCCC data. To evaluate the ability of MODEMDR to handle large datasets, a large dataset was obtained from the WTCCC 24 , consisting of 500,569 SNPs, including 1,988 cases of coronary artery disease (CAD) and 1,500 controls obtained from people living in Great Britain who self-identified as white Europeans.  Table 2. The gene names in Table 2 were obtained from dbSNP at the National Center for Biotechnology Information. The designation "UNKNOWN" in the table refers to SNP not being located on a gene. The minimum and maximum numbers of detected significant epistatic interactions in each chromosome were 1 and 6, respectively. The p value was estimated through an χ 2 test using the raw datasets to determine the significance level for epistatic interaction between two SNPs 11 . All epistatic interactions detected by MODEMDR in 24 chromosomes yielded a p value of <0.0001, indicating a strong significant interaction between two SNPs. When the CCR was larger than 0.5, the frequency of chance can be significantly reduced 27 , indicating that our results identified significant epistatic interactions. The high NMI values indicate that uncertainty was reduced in the model in a true state. The CCR and NMI values show the multiobjective optimization property. The epistatic interaction with the highest NMI value (rs41399650, rs41397248) was not the epistatic interaction with the highest CCR value (rs41399650, rs17163057) in chromosome 1. Further detection of significant epistatic interactions may provide an etiological understanding of epistasis in systems biology 28 . The MODEMDR for CCR measures was located between 0.588 and 0.959 and the mean CCR was 0.750 (standard deviation (SD) = 0.096). The NMI measures were located between 0.033 and 0.759 and the mean NMI was 0.267 (SD = 0.182). Notably, the epistatic interaction of SNPs rs16926425 and rs7299571 (chromosome 12) obtained the highest CCR (0.959) and NMI values (0.759). The ten detected epistatic interactions indicate the beneficial measures of NMI > 0.4 and CCR > 0.8 (marked by stars in Table 2). The details of all epistatic interactions are shown in Supplementary Fig. 1. In all figures, the (black) left bar in a class represents the frequency of cases and the (white) right bar represents the frequency of controls. Gray classes indicate being in the high-risk group.  Table 1. Comparison of SNPRuler, MDR, DEMDR, and MODEMDR across 60 epistasis models using Wilcoxon signed-rank testing. R − : number of epistasis models when MODEMDR lost another algorithm; R + : number of epistasis models when MODEMDR won another algorithm; R = : number of epistasis models when two algorithms tied; N: number of R − , R + , and R = . P < 0.05 indicates a significant difference between two algorithms. The running times of chromosomes in the WTCCC dataset are shown in Table 2

Discussion
MODE enables MDR to use multiple measures to detect gene-gene interactions. Although the CCR in MDR-based methods is a powerful measure for determining such interactions, it could fail to determine interactions in some epistatic models (e.g., models 31-60 without marginal effects (Fig. 5)). Furthermore, the NMI could not always determine specific targets within models 31-60 without marginal effects. Therefore, both measures were considered for synchronous use to effectively determine the targets (Table 2). MODEMDR can effectively detect gene-gene interactions because the MODE fits the joint effect property 29 , which consists of the main effect, overall effect, and high-order interaction effect. The main effect refers to any effect that could serve as a guide to identifying the correct multilocus interaction. The overall effect refers to an effect that commonly appears among n risk factors. The high-order interaction effect refers to the least proper subset of the loci that interacts epistatically. SNPs strongly associated with diseases or cancers are often likely to be significant factors in high-locus interactions. High CCR and NMI values in MODEMDR indicate a more significant risk of n-factor effects. In the MODEMDR selection operation, promising SNPs can be retained for the next generation. These SNPs are subsequently combined through mutation and recombination operations to produce better SNP combinations, enabling MODEMDR to detect the significant epistatic interactions.
In MDR, combinations of high-dimensional factors can be reduced by assigning multilocus genotypes to high-or low-risk groups, enabling gene-gene interaction quality to be measured through two-way contingency table analysis 10 . CCR is the measure most commonly applied in MDR-based methods 30 . Bush et al. (2008) compared the ten general measures in the text classification field to evaluate the degree of improvement in the ability of MDR to detect gene-gene interactions. CCR and NMI were suggested as being able to improve MDR identification in the simulation 22 . The results of the present study exhibited the most successful gene-gene interaction identification when the NMI and CCR were used to synchronously determine significant gene-gene interactions.
The WTCCC dataset is well-known in GWAS analyses, in which large SNPs in chromosomes are collected. MODEMDR efficiently identifies gene-gene interactions from combinatorially explosive search spaces (running time in Table 2), and uses the rational performing time (population size × generation size) to calculate MDR measures, enables it to handle GWAS analysis. MODEMDR has the advantages of MDR because the fitness functions of multiple objectives are designed based on MDR. The advantages of MODEMDR include the following: (i) suitability for application in small sample datasets, (ii) suitability for application in unbalanced datasets, (iii) ability to describe the loci genotype combinations associated with high and low risk of disease, (iv) the model-free method, (v) ability to detect a higher-order gene-gene interactions, and (vi) the nonparametric method.
MODEMDR was applied in this study for synchronous consideration of the multiple measures used to detect significant gene-gene interactions. To our knowledge, MODEMDR is the first MDR-based method that accounts for multiple measures. The experimental results demonstrate that multiple measures engender an identification performance superior to that of the MDR-based single measure method. In the WTCCC analysis, MODEMDR successfully handled the large-scale dataset in terms of speed and identification of significant gene-gene interactions. Furthermore, MODEMDR provides a multiobjective method for identifying gene-gene interactions. Improvements could be made by using more combinations among various measures in a two-way contingency table. .
where X is the decision vector and f is the objective vector. For X i , all objectives f that are not dominated by any other vector X j (j = 1, …, k | i≠j) where k is the population size are called nondominated points. For gene-gene interaction identification, we defined "gene-gene interaction" (i.e., solution) as a decision vector and "measures" as the corresponding objective vector. Here, CCR 10 and NMI 22 were defined as f 1 and f 2 , respectively. Therefore, in this study, the objective was defined as follows: where X is the solution space and X i ∈ X.

MODEMDR.
In MODEMDR, the MDR operation process is modified to apply MDR as a fitness function in MODE. In addition, a balance strategy is introduced in data preprocessing within cross-validation in MDR to improve the accuracy of fitness evaluation. The balance strategy can effectively increase the CCR in the training and testing. In MODE operations, target vector X, mutant vector V, and trial vector U are used to seek the optimal multiobjective set. Data preprocessing. Data preprocessing uses the balance strategy in MDR to handle the balanced k-fold CV subsets that are divided by the original dataset for objective evaluation. For k-fold CV operation, the balanced k-fold CV subsets are generated through the following five steps of the balance strategy: 1.
Step 1. Divide the samples into case sets (cases) and control sets (controls). 2.
Step 2. Randomly shuffle the case and control samples. 3.
Step 3. Count the total numbers of cases and controls. 4.
Step 4. Compute the ratio between cases and controls. 5.
Step 5. Assign the case and control samples to subset j (j = 1, …, k) according to the ratio, where j is the CV index and k is the total number of CV subsets. 1.
Step 1. Comparison between an X ∈ population and s j for all indices j ∈ (1, …, i) in S. If X is not dominated by any s j , X is added into S. 2.
Step 2. Comparison between an s j for index j ∈ (1, …, i) and s k for indices k and k ∈ (1, …, i | k ≠ j) in S. If s j is dominated by any s k , s j is discarded.
Target vector definition. Let X i,g = (x 1,i,g , .., x d,i,g ) be the ith target vector in the population for the gth generation in the d-dimensional search space. A target vector is a gene-gene interaction in which the parameters are the SNP indices and are all different in a target vector. Given that y-SNPs and d-order gene-gene interaction identification are considered in case-control studies, the target vector X i is represented as follows: , )), where g is the gth generation. For initialization (i.e., g = 0), the parameter x j (j = 1, …, d) in the target vector X i is randomly generated by (2): where upper and lower represent the upper and lower boundaries of the indices of independent variables, respectively. The rand j,i is the random number generator, which returns a uniformly distributed random value from within the range [0, 1).

Mutation. Each target vector generates a mutant vector V i,g+1
, which is a vector sum of the weighted difference between two vectors and a third vector, expressed as follows: In (3), n is the population size; r 1 , r 2 , and r 3 ∈ (1, …, n) are the random indices of the storage (Pareto operation) or population; g is the gth generation; X r1,g , X r2,g , and X r3,g are the selected three target vectors from the storage or population, where all selected target vectors are different; and F is a real and constant factor ∈ [0, 2) that controls the amplification of the differential variation (X r2,g − X r3,g ). In (4), S ri,g is the r i th target vector [r i ∈ (1, …, n)] in the storage and PV is a mutation constant ∈ [0, 1) that controls the probability of the vector selected from either the storage or population.

Recombination.
Recombination operation can increase the vector diversity in the population. The trial vector U i,g+1 is expressed as (5) and the parameters of the trial vector are computed by (6), which incorporates the mutant vector V i,g+1 and current target vector X i,g at the ith target vector, expressed as follows:  In (5), i is the trial vector index in the population, d is the dimension size, and g is the gth generation. In (6), j is the index of the dimension in the mutant vector V i,g and target vector X i,g , where the two is represent the indices of the mutant vector and target vector in the population; randb(j) is the jth evaluation of a uniform random number generator with the outcome ∈ [0, 1); CR is the crossover constant ∈ [0, 1); and rnbr(i) is a randomly chosen index ∈ (1,…, d) that ensures that U i,g+1 obtains at least one parameter from V i,g+1 .
Boundary constraints. Boundary constraints can ensure that trial vectors are feasible combinations. Equation (7) guarantees that trial vector parameters do not violate boundary constraints with random values generated by (2), expressed as follows: where j is an index of the dimension in the trial vector U i,g , i is the index of the trial vector in the population, g is the gth generation, upper and lower are the upper and lower bounds of the indices of independent variables, respectively, and ∃!u j,i,g represents a variable at the jth parameter only existing in the ith trial vector for the gth generation.
Selection. Selection operation determines whether the target vector X i,g is dominated by the trial vector U i,g ; in other words, f j (U i,g ) > f j (X i,g ) for all indices j ∈ (1, 2), where j is the index of the objective function. If the trial vector U i,g+1 dominates the target vector X i,g , X i,g+1 is set to U i,g+1 , otherwise X i,g is retained as X i,g+1 . In (1) ., x d )) contain d 3 multifactor cells, each of which contains the total quantities of cases and controls for the corresponding genotype combination.
Step 1. Determine high or low risk within multifactor cells by using the training data. Each multifactor cell is deemed high or low risk by evaluating the ratio between total quantities of cases and controls in that cell. A cell is deemed high-risk if ratio ≤ 1 and low risk otherwise. In the training data, θ  a represents a ratio value and is computed by (8) to provide a more accurate ratio to determine whether a cell is high or low risk. Thus, accurate objective evaluations can be improved when the total quantities of cases and controls are unbalanced. Equation (8)  where n ab is the total number of samples within the ath multifactor cell in the b outcome risk status in the training data, and n +b represents the total number of samples in the b outcome risk status, where b = 1 for cases and 0 for controls.
Step 2. Determine high or low risk within multifactor cells by using the testing data.
To use the testing data to determine whether multifactor cells are high or low risk, the ratio θ a is computed by (9), expressed as follows: where t ab is the number of samples within the ath multifactor cell in the b outcome risk status in the testing data, where b = 1 for cases and 0 for controls. Both n +0 and n +1 are the same as in (8).
Step 3. Evaluate the true positive (TP), false positive (FP), false negative (FN), and true negative (TN) values by comparing the level of risk in multifactor cells as determined by the training and testing data. A comparison of the risk level of a single cell as determined by training and testing data can be used to compute the TP, FP, FN, and TN. Thus, all multifactor cells can be reduced to four dimensions (TP, FP, FN, and TN). Equation (10) expresses the evaluation functions of the four dimensions as follows: where t ab is the number of samples within the ath multifactor cell in the b outcome risk status; n +b is the total number of samples in the b outcome risk status, where b = 1 for cases and 0 for controls; TP is the number of correctly classified samples in the testing data within the high-risk range as determined by training data; FP is the number of incorrectly classified samples in the testing data within the low-risk range as determined by the training data; FN is the number of incorrectly classified samples in the testing data within the high-risk range as determined by the training data; TN is the number of correctly classified samples in the testing data within the low-risk range as determined by the training data.
Step 4. Evaluate the fitness functions of objectives. Objective 1: Objective 1 is the CCR (11), which is used to determine the proportion of correctly classified individuals. The CCR is computed using the TP ratio for cases and TN ratio for controls, where the maximum value indicates the optimal solution. Equation (11) (11) where TP, FP, FN, and TN are computed using (10). Objective 2: Bush et al. used the NMI to evaluate MDR. NMI is a measure of information transmission based on Shannon entropy, interpreted as the proportion of information contained in the row variable transferred or transmitted to the column variable; more concisely, it is the amount by which the model reduces our uncertainty about the true state 22 . In the 2 × 2 contingency table, the row entropy H(x), column entropy H(y), and conditional entropy H(y|x) are defined as (12), (13), and (14), respectively, and expressed as follows: where TP, FP, FN, and TN are computed using (10), with the maximum value indicating the optimal solution.
Step 5. Repeat steps 1-4 until all CV folds have been completed.
Step 6. Compute the averages of the CCR and NMI values in all CV folds.
Illustrative example of MODEMDR. The supplemental material in this paper provides an example of how MODEMDR works.
Parameter settings. The SNPRuler parameter is set to the default settings. The parameter "updateRatio" is set to 0.2, which is the step size used for updating a rule. MDR, DEMDR, and MODEMDR use the five-fold CV test. DEMDR and MODEMDR have the following four common parameters: population size (pop-size), generation size (gen-size), scaling factor (F), and crossover constant (CR). For the