HSIC_CR: A Lightweight Scoring Criterion Based on Measuring the Degree of Causality for the Detection of SNP Interactions

Abstract

Recently, research on detecting SNP interactions has attracted considerable attention, which is of great significance for exploring complex diseases. The formulation of effective swarm intelligence optimization algorithms is a primary resolution to this issue. To achieve this goal, an important problem needs to be solved in advance; that is, designing and selecting lightweight scoring criteria that can be calculated in $O(m)$ time and can accurately estimate the degree of association between SNP combinations and disease status. In this study, we propose a high-accuracy scoring criterion (HSIC${}_{CR}$) by measuring the degree of causality dedicated to assessing the degree. First, we approximate two kinds of dependencies according to the structural equation of the causal relationship between epistasis SNP combination and disease status. Then, inspired by these dependencies, we put forward this scoring criterion that integrates a widely used method of measuring statistical dependencies based on kernel functions (HSIC). However, the computing time complexity of HSIC is $O(m^2)$, which is too costly to be an integral part of the scoring criterion. Since the sizes of the sample space of the disease status, SNP loci and SNP combination are small enough, we propose an efficient method of computing HSIC for variables with a small sample in $O(m)$ time. Eventually, HSIC${}_{CR}$ can be computed in $O(m)$ time in practice. Finally, we compared HSIC${}_{CR}$ with five representative high-accuracy scoring criteria that detect SNP interactions for 49 simulation disease models. The experimental results show that the accuracy of our proposed scoring criterion is, overall, state-of-the-art.

1. Introduction

Since many complex diseases are usually caused by multiple genes and multiple factors, in recent years, with the emergence of high-throughput genotypic technology, genome-wide association analysis (GWAS) has been one of the main methods used to study complex diseases. Furthermore, the identification of single-nucleotide polymorphism (SNP) interactions from GWAS data is of great importance for exploring the explanation, prevention and treatment of complex diseases [1,2]. Therefore, over the past decade, this research topic has attracted considerable attention [3,4,5,6,7,8,9].

It is well known that SNP interactions represent combinations of multiple SNPs that affect complex diseases in a linear or non-linear manner, also known as k-order epistasis SNPs. The research topic of detecting k-order epistasis SNPs is a typical case of combinatorial optimization problems in k-dimensional discrete space ($k\in \{2,3,4,5\}$ in practice), and swarm intelligence optimization (SIO) algorithms are one of the main methods used to solve the problems [9,10,11]. For this study to be successful, an important problem needs to be solved in advance; that is, designing and selecting lightweight scoring criteria that can be calculated in $O(m)$ time and can accurately estimate the degree of association between SNP combinations and disease status.

To date, few lightweight scoring criteria can accurately estimate the degree of association of SNP combinations with disease status in most disease models due to the widely varying characteristics of different disease models. As one of the primary methods used to work on this combinatorial optimization problem, SIO algorithms mostly tackle this issue by combining multiple criteria [9,10,11,12,13]; however, using too many objective functions will often make the proposed algorithm difficult to converge effectively. Therefore, picking a few high-accuracy objective functions instead of using too many objective functions can dramatically improve the performance of the used algorithms [14,15].

This paper’s goal is not to contribute toward fixing the issue entirely but to use a different methodology to propose a scoring criterion that can accurately estimate the associations in most disease models. The contributions of this paper are:

• We propose a high-accuracy scoring criterion based on measuring the degree of causality that integrates a widely used method of measuring statistical dependencies (HSIC);
• We put forward an efficient algorithm of computing HSIC on two variables with a small sample in $O(m)$ time, thus enabling us to compute HSIC${}_{CR}$ in $O(m)$ time in practice.

2. Related Works

So far, the proposed lightweight scoring criteria can be roughly divided into two categories.

The first category covers various approaches, which are so-called Bayesian scoring criteria. The Bayesian scoring criteria calculate the posterior probability distribution, proceeding from a prior belief on the possible DAG models, conditional on the data [16]. The K2-Score is an efficient Bayesian scoring criterion that obtains priors under the assumption that all DAG models are equally likely. Other such scoring criteria representatives include the Bayesian Dirichlet equivalent (BDe) scoring criterion and the Bayesian Dirichlet equivalent uniform (BDeu) scoring criterion [4,17,18,19].

The second category is usually known as the information-theoretic scoring criteria. Mutual information (Mi) is a lightweight method but has preferences for certain disease models [6,20]. The JS divergence is a symmetrized divergence measure, derived from the Kullback–Leibler (KL) divergence, which is an asymmetric divergence measure of two probability distributions [21]. This approach can be utilized to evaluate the SNP genotype deviation between control samples and case samples. Lately, joint entropy (JE) and normalized distance with joint entropy (ND-JE) have been proposed as criteria for guiding harmony search algorithms to discover clues for exploring the epistasis of SNP combinations [9].

There are also a few approaches that do not fall into any of the above two main categories. For example, the LR is a composite indicator that reflects both sensitivity and specificity and can be used for a related measure to find the likelihood difference between a disease-causing SNP combination and an SNP combination that is not involved in the disease process [22,23].

In statistics, G-test is a significant test method of natural ratio or maximum likelihood. In recent years, scholars have tended to use the G-test independence test instead of the chi-square independence test recommended in the past. In genome association analysis, G-test has been extensively used. Different from other scoring criteria, G-test will provide its p-value when measuring the relationship between an SNP combination and sample state, which can indicate whether the SNP combination has a significant relationship with the sample state [24].

Published research has found that the results were different when employing different scoring criteria, The K2-Score has been widely used to evaluate the association. This measure has a high capacity for detecting SNP interactions and is superior in discriminating certain disease models with low marginal effects. However, for the interaction model with low minor allele frequencies (MAFs) and low genetic heritability (h), the K2-Score has a low performance in detecting high-order SNP interactions. The ND-JE is proposed based on the properties of the disease-causing SNP combination models without marginal effects, so this metric is more suitable for evaluating diseases with this type of mode. The LR-score aims to discover the relationship between likelihood differences in functional SNP combinations and non-functional SNP combinations. The method is well adapted to unbalanced datasets of cases and controls. In practice, the use of the G-test as a single evaluation criterion for detection is found to be inadequate, as there are often many SNP combinations with G-test values close to 0 [6,9,11].

The scoring criterion proposed in this work is based on the theory of causality. It is distinct from the theoretical approach taken by the current existing criteria. From the perspective of a comparison with correlation, causality strictly distinguishes “cause” variables and “result” variables, and plays an irreplaceable role in revealing the mechanism of the occurrence of things and guiding intervention behaviors [25]. Thus, the proposed criterion is a useful for and complementary to the current existing criteria.

3. Methodology

3.1. Concepts and Terms

In this work, x = {$x_1$, $x_2$, …, $x_n$} represents a set of n SNP loci, and X = $\left[\begin{array}{cccc}{x}_{11}& x_{21}& .& x_{1n}\\ .& .& .& .\\ .& .& .& .\\ x_{m1}& x_{m2}& .& x_{mn}\end{array}\right]$ is a set of m samples of x; $Y={\{y_1,y_2,\dots ,y_m\}}^t$ denotes a set of m samples of disease status y. For $\forall 1⩽i⩽n$, $1⩽j⩽m$, $i,j\in \mathbb{Z}$, $x_{ij}\in$ {0, 1, 2}, $x_{ij}$ is equal to 0, 1 and 2, which implies that it is the homozygous major allele (AA), heterozygous allele (AT), and homozygous minor allele (TT), respectively; $y_j=0$ for control and $y_j=1$ for case; D = (X,Y) is a dataset with m samples.

Definition 1

(k-order epistasis SNP combination). Let $S_k$ = {{$x_{i_1}$, $x_{i_2}$, …, $x_{i_k}$}} be a collection of a set with k SNP loci ($1<k<n$). $f(D)$: $S_k$ $\to {\mathbb{R}}^{+}$ is a score function used for measuring the association between any k SNP loci and disease status y based on a dataset D. If x has the k-order epistasis SNP combination on y (denoted as s${}_k^{\ast },s_k^{\ast }\in S_k$), and $f(D)$ is a correct score function (or scoring criterion), then, for $\forall s_k\in S_k$, s${}_k\ne s_k^{\ast }$, $f(D)$(s${}_k^{\ast }$) < $f(D)$(s${}_k$) or $f(D)$(s${}_k^{\ast }$) > $f(D)$(s${}_k$).

3.2. Causal Relationship

According to how the data are generated, the structural equation of the causal relationship between epistasis SNP combination and disease status can be modeled as [26]:

$$y=f(s_k^{\ast })+e_y,$$

(1)

where $e_y$ is the noise variable.

From the above equation, we can find that $s_k^{\ast }$ and $e_y$ are independent (denoted as $s_k^{\ast }\perp e_y$). In other words, among all $s_k\in S_k$, $s_k^{\ast }$ and $e_y$ have the lowest degree of dependency. However, it is unrealistic to measure the dependence degree of any $s_k$ and $e_y$ as the evaluation criterion for epistasis detection, because it requires too high a computational cost to obtain data generated by $e_y$ based on the regression method.

Thus, this paper herein let $e_y$ be the constant 0, i.e., y $\approx f(s_k^{\ast })$, which approximately introduces two kinds of dependencies as described in Figure 1 and Figure 2, respectively [25]. Obviously, the dependence between $s^{\ast }$ and y is direct (denoted as $s_k^{\ast }\cancel{\perp }y$); and the other is derived from the v-structure, i.e., $x_{i_1}$ and $x_{i_2}$ are dependent given a value of y and $s_{k(i_1{i}_2)}^{\ast }$ (denoted as $x_{i_1}\cancel{\perp }{x}_{i_2}|y=y_i,s_{k(i_1{i}_2)}^{\ast }=c_j$). Let $s_k^{\ast }=\{x_{i_1},x_{i_2},\dots ,x_{i_{k-1}},x_{i_k}\}$, $s_{k(i_1{i}_2)}^{\ast }$ represent $s_k^{\ast }\backslash \{x_{i_1},x_{i_2}\}$. In particular, the set $s_{k(i_1{i}_2)}^{\ast }$ is empty when k is equal to 2.

3.3. Scoring Criterion

These two kinds of dependencies described above inspire us to raise this scoring criterion, which integrates a widely used method of measuring statistical dependencies based on kernel functions (HSIC).

For $\forall s_k\in S_k$, then, for $\forall x_{i_1},x_{i_2}\in s_k$, given $y=p$ ($p\in \{0,1\}$) and $s_{k(i_1{i}_2)}=q$ ($q\in \{c_1,c_2,\dots ,c_{3^{k-2}}\}$), let $D_{i_1{i}_2}^{pq}=(X_{i_1{i}_2}^{pq},Y_{i_1{i}_2}^{pq})$ be a slice of D on $x_{i_1},x_{i_2}$ under the constraint; $m_{i_1{i}_2}^{pq}$ is the number of rows of the data slice ($m_{i_1{i}_2}^q=m_{i_1{i}_2}^{0q}+m_{i_1{i}_2}^{1q}$); HSIC$(X,Y)$ is used to measure the degree of statistical dependence of two random variables (x and y) based on dataset (X, Y); the scoring criterion can be computed by the following Equations (2)–(6).

$$\begin{array}{cc}\hfill b_{i_1{i}_2}^q& =\frac{m_{i_1{i}_2}^{0q}}{m_{i_1{i}_2}^q}\mathrm{HSIC}(X_{i_1}^{0q},X_{i_2}^{0q})+\frac{m_{i_1{i}_2}^{1q}}{m_{i_1{i}_2}^q}\mathrm{HSIC}(X_{i_1}^{1q},X_{i_2}^{1q})\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill b^q& =\sum _{i_1\ne i_2,x_{i_1},x_{i_2}\in s_k}(\frac{m_{i_1{i}_2}^q}{{\sum }_{i_a\ne i_b,x_{i_a},x_{i_b}\in s_k}{m}_{i_a{i}_b}^q}{b}_{i_1{i}_2}^q)\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\overline{m}}^q& =(\sum _{i_a\ne i_b,x_{i_a},x_{i_b}\in s_k}{m}_{i_a{i}_b}^q)/c_k^2\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\mathrm{HSIC}}_{CR}^q& =\frac{m}{m+{\overline{m}}^q}\mathrm{HSIC}(X_{i_1{i}_2\dots ,i_k},Y_{i_1{i}_2\dots ,i_k})+\frac{{\overline{m}}^q}{m+{\overline{m}}^q}{b}^q\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\mathrm{HSIC}}_{CR}(X_{i_1{i}_2\dots ,i_k},Y_{i_1{i}_2\dots ,i_k}:data)& =\sum _{q=1}^{q=3^{k-2}}{\mathrm{HSIC}}_{CR}^q\hfill \end{array}$$

(6)

To facilitate the reader’s understanding, we now define following notations:

1.

The value of $b_{i_1{i}_2}^q$ is a linear weighted sum of $\mathrm{HSIC}(X_{i_1}^{0q},X_{i_2}^{0q})$ and $\mathrm{HSIC}(X_{i_1}^{1q},X_{i_2}^{1q})$ based on respective sample sizes;

2.

For all $x_{i_1},x_{i_2}\in s_k$, the value of $b^q$ is a linear weighted sum of all $b_{i_1{i}_2}^q$ since there are $c_k^2$ v-structures given $s_{k(i_1{i}_2)}=q$;

3.

The value of ${\mathrm{HSIC}}_{CR}^q$ is a linear weighted sum of $\mathrm{HSIC}(X_{i_1{i}_2\dots ,i_k},Y_{i_1{i}_2\dots ,i_k})$ and $b^q$ based on sample size of $D_{i_1{i}_2\dots ,i_k}$ and average sample size of all $D_{i_a{i}_b}^q$, which is a component and basis of ${\mathrm{HSIC}}_{CR}$;

4.

In particular, $b_{i_1{i}_2}^q=b^q$ as is ${\mathrm{HSIC}}_{CR}^q={\mathrm{HSIC}}_{CR}$ when $k=2$;

5.

For robustness purposes, let $b_{i_1{i}_2}^q$ = 0 if and only if the denominator of the weighted factor term is 0, like $b^q$;

6.

The effort to calculate scoring criterion is reduced to calculate $\mathrm{HSIC}(X_{i_1{i}_2\dots ,i_k},Y_{i_1{i}_2\dots ,i_k})$ once, and is reduced up to $C_k^2\ast 3^{k-2}$ times to calculate the type of problem $b_{i_1{i}_2}^q$ (fortunately, $k\in \{2,3,4,5\}$ in practice).

Thus, our estimate of $s_k^{\ast }$ can eventually be obtained by solving the problem

$$\begin{array}{c}\hfill \underset{s_k\in S_k}{maxf(D)(s_k)}={\mathrm{HSIC}}_{CR}(X_{i_1{i}_2\dots ,i_k},Y_{i_1{i}_2\dots ,i_k})\end{array}$$

(7)

3.4. Method for Measuring Statistical Dependence

3.4.1. HSIC

HSIC is a measuring statistical dependence criterion proposed by other authors [27,28] based on the eigenspectrum of covariance operators in reproducing kernel Hilbert spaces (RKHSs), denoted by $\mathrm{HSIC}\left(P_{xy}\right)$ as follows:

$$\begin{array}{cc}\hfill \mathrm{HSIC}\left(P_{xy}\right)& =E_{x,x{}^{\prime },y,y{}^{\prime }}\left[k(x,x{}^{\prime })l(y,y{}^{\prime })\right]+\hfill \\ \hfill & E_{x,x{}^{\prime }}\left[k(x,x{}^{\prime })\right]E_{y,y{}^{\prime }}\left[l(y,y{}^{\prime })\right]-\hfill \\ \hfill & 2E_{x,y}\left[E_{x^{{}^{\prime }}}\left[k(x,x{}^{\prime })\right]E_{y^{{}^{\prime }}}\left[l(y,y{}^{\prime })\right]\right],\hfill \end{array}$$

(8)

where $k(,)$ and $l(,)$ are two kernel functions.

Let $\mathbb{X}$ and $\mathbb{Y}$ be the separable sample spaces of random variables x and y, respectively, assuming that ($\mathbb{X},\mathsf{\Gamma }$) and ($\mathbb{Y},\mathsf{\Lambda }$) are furnished with probability measures $p_x,p_y$, respectively ($\mathsf{\Gamma }$ being the Borel sets on $\mathbb{X}$, and $\mathsf{\Lambda }$ the Borel sets on $\mathbb{Y}$); $p_{xy}$ is a joint measure over $(\mathbb{X}\times \mathbb{Y},\mathsf{\Lambda }\times \mathsf{\Gamma })$; $\mathrm{HSIC}\left(P_{xy}\right)\ge 0$ (the higher the degree of dependence of x and y, the greater the value) and $\mathrm{HSIC}\left(P_{xy}\right)$ is zero if and only if x and y are independent.

In order to show that HSIC is a practical criterion for measuring independence or the degree of dependence given a finite number of observations, it consists of an empirical estimator with $O(m^{-1})$ expectation bias, denoted by $\mathrm{HSIC}\left(D\right)$, formulated as follows:

$$\begin{array}{cc}\mathrm{HSIC}\left(D\right)& ={(m-1)}^{-2}trace(KHLH),\end{array}$$

(9)

where D:=$\{(x_1,y_1),\dots ,(x_m,y_m)\}⊑\mathbb{X}\times \mathbb{Y}$, $K,H,L\in {\mathbb{R}}^{m\times m}$,$K_{i,j}:=k(x_i,x_j)$, $L_{i,j}:=l(y_i,y_j)$, $H_{i,j}:={\delta }_{ij}-m^{-1}$.

An advantage of HSIC compared with other kernel-based independence criteria is that it can be computed in $O(m^2)$ time. However, such computational costs are too high as an integral part of the scoring criterion. Fortunately, the sample space of the disease status, SNP loci and k-order SNP combination is finite discrete. Thus, immediately below, we put forward an efficient HSIC calculation method for variables with a small sample, which can be approximately calculated in $O(m)$ time. Thus, as $k\in \{2,3,4,5\}$, we can compute ${\mathrm{HSIC}}_{CR}$ in $O(c_k^2\times m)\approx O(m)$ time in practice.

3.4.2. Efficient Computation

Proposition 1

(Efficient computation). Let x and y be two random discrete variables with p and q states, respectively, where $p^2\times q^2<m$, or $p^2\times q^2\approx m$. Then, we can compute $trace(KHLH)$ in $O(m)$ time.

Proof. 

Let e be a column vector with a length of m, $e=\left[\begin{array}{c}1\\ .\\ .\\ 1\end{array}\right],$ $L=\left[\begin{array}{c}{l}_1{}^t\\ .\\ .\\ l_m{}^t\end{array}\right]$, and I be an identity matrix with a size of $m\ast m$; we have $H=I-\frac{1}{m}ee{}^t$ and $LH=L-\frac{1}{m}Lee{}^t$.

As $\frac{1}{m}Le=\left[\begin{array}{c}\frac{1}{m}{l}_1{}^{t}e\\ .\\ .\\ \frac{1}{m}{l}_m{}^{t}e\end{array}\right]$, which implies that each i-th element of $\frac{1}{m}Le$ is the mean of the corresponding row elements of L, we have $L-\frac{1}{m}Lee{}^t=\left[\begin{array}{c}{l}_1{}^t\\ .\\ .\\ l_m{}^t\end{array}\right]-\left[\begin{array}{c}\frac{1}{m}\overline{l_1}\\ .\\ .\\ \frac{1}{m}\overline{l_m}\end{array}\right]e{}^t$, where each $\overline{l_i}$ is the sum of the corresponding row elements of L.

Let $\overline{L}=\left[\begin{array}{c}\frac{1}{m}\overline{l_1}\\ .\\ .\\ \frac{1}{m}\overline{l_m}\end{array}\right]e{}^t$ be an $m\times m$ matrix (${\overline{L}}_{ij}=\frac{1}{m}\overline{l_i}$) and $\overline{K}=\left[\begin{array}{c}\frac{1}{m}\overline{k_1}\\ .\\ .\\ \frac{1}{m}\overline{k_m}\end{array}\right]e{}^t$ be an $m\times m$ matrix (${\overline{K}}_{ij}=\frac{1}{m}\overline{k_i}$); we have $trace(KHLH)=trace((K-\overline{K})(L-\overline{L}))$.

Let P be an $m\times m$ row transformation matrix; we have $trace(PKHPLH)=trace(P(K-\overline{K})P(L-\overline{L}))=trace((K-\overline{K})(L-\overline{L}))$=$trace(KHLH)$.

Thus, without a loss of generality, we can assume that: D=$\{(c_1,y_{i_1}),\dots (c_1,y_{i_2-1}),$ $(c_2,y_{i_2}),\dots (c_2,y_{i_3-1}),\dots (c_p,y_{i_p}),\dots (c_p,y_{i_m})\}$, i.e., the number of observed instances of x with the value of $c_j$ (denoted as $xt(c_j)$) is $i_{j+1}-i_j$($i_{p+1}=i_m+1$); K can be viewed as a $p\times p$ partitioned matrix. Let $K^{jl}$ be the $j\times l$th block having $xt(c_j)\times xt(c_l)$ elements, all having the same value ( $k(c_j,c_l)$). Let $\widehat{K}=K-\overline{K}$, where all elements in ${\widehat{K}}^{jl}$ have the same value equal to $k(c_j,c_l)-\frac{{\sum }_{n=1}^{p}xt(c_n)\times k(c_j,c_n)}{m}$ (denoted as $\widehat{K}(j,l)$).

Let $yx(j,i)$ be the number of observed instances with the value $(c_i,d_j)$ ($1\le j\le q$, $d_j$ is the j-th state of y); the definition of $\widehat{L}$ and $\overline{L}$ is similar to that of $\widehat{K}$ and $\overline{K}$. We also view $\widehat{L}$ as a $p\times p$ partitioned matrix, where each $L^{jl}$ has the same number of rows and columns as $K^{jl}$; for $\forall p,1\le p\le m$, $\frac{1}{m}{\overline{l}}_p=\frac{{\sum }_{v=1}^{q}yt(d_v)\times k(d_h,d_v)}{m}$ (denoted as $\overline{l}(d_h)$), where $y_p=d_h$.

As $trace(\widehat{K}\widehat{L})={\sum }_{i=1}^p{\sum }_{j=1}^p<{\widehat{K}}^{ij},{({\widehat{L}}^{ji})}^T>$ (<.,.> denoted as inner product operator) and $<{\widehat{K}}^{ij},{({\widehat{L}}^{ji})}^T>=\widehat{K}(i,j)\times {\sum }_{u=1}^{q}yx(u,j)\times {\sum }_{v=1}^q(yx(v,i)\times (l(d_u,d_v)-\overline{l}(d_u))$, we can obtain that the computational complexity of $trace(\widehat{K}\widehat{L})$ is $O(p^2{q}^2)$.

As described above, we can know that the total computational complexity of $xt$, $yt$ and $yx$ is $O(m)$, and that those of $\overline{K}$ and $\overline{L}$ are $O(p^2)$ and $O(q^2)$, respectively.

Hence, we have that the total computational complexity of $\mathrm{HSIC}$ is $O(m)$, where $p^2\times q^2<m$, or $p^2\times q^2\approx m$.

The proof is complete.    □

In fact, the proof above gives the simplified process of efficient computation to $\mathrm{HSIC}$. The detailed processes are shown in Algorithms 1–5. Algorithm 1 is the main process of the method, consisting of three functions:

• (X, Y) is m observations of a tuple of x and y with p and q states, respectively;
• $kernels$ includes two kernel functions used to calculate $K_{i,j}$ and $L_{i,j}$, the parameters of which are $delta(1)$ and $delta(2)$, respectively;
• $GetInfo$ (see Algorithm 2) is used to calculate $xt$, $yt$ and $yx$;
• For all $1\le j,l\le p$, $KH$ (see Algorithm 3) is used to calculate $\widehat{K}(j,l)$;
• $Trace$ (see Algorithm 4) is used to calculate $trace(\widehat{K}\widehat{L})$;
• $RowAverage$ (see Algorithm 5) is used to calculate $\overline{K}$ and $\overline{L}$.

Algorithm 1 Calculate $value=\mathrm{HSIC}(X,Y,p,q,m,kernels,deltas)$

Require: $\left|X\right|=\left|Y\right|=m$ 1: $[xt,yt,yx]\Leftarrow GetInfo(X,Y,p,q,m)$ 2: $\widehat{K}\Leftarrow KH(xt,p,kernels(1),deltas(1),m)$ 3: $value\Leftarrow Trace(\widehat{K},yt,yx,p,q,kernels(2),deltas(2),m)$

Algorithm 2 Calculate $[xt,yt,yx]$= $GetInfo(X,Y,p,q,m)$

1: $xt\Leftarrow zeros(1,p)$   2: $yt\Leftarrow zeros(1,q)$   3: $yx\Leftarrow zeros(q,p)$   4: $col\Leftarrow 1$   5: while$col\le m$  do   6:     $stx\Leftarrow X(col)$   7:     $sty\Leftarrow Y(col)$   8:     $xt(stx)\Leftarrow x(stx)+1$   9:     $yt(sty)\Leftarrow y(sty)+1$ 10:     $yx(sty,stx)\Leftarrow yx(sty,stx)+1$ 11:     $col\Leftarrow col+1$ 12: end while

Algorithm 3 Calculate $\widehat{K}=KH(xt,p,kernel,delta,m)$

1: $\overline{K}\Leftarrow RowAverage(p,k,delta,m,xt)$   2: $\widehat{K}\Leftarrow zeros(p,p)$   3: $i\Leftarrow 1$   4: while$i\le p$  do   5:     $j\Leftarrow 1$   6:     while $j\le p$ do   7:         $\widehat{K}(i,j)\Leftarrow kernel(i,j,delta)-\overline{K}(i)$   8:         $j\Leftarrow j+1$   9:     end while 10:     $i\Leftarrow i+1$ 11: end while

Algorithm 4 Calculate $value=Trace(\widehat{K},yt,yx,p,q,kernel,delta,m)$

1: $\overline{L}\Leftarrow RowAverage(q,kernel,delta,m,yt)$   2: $value\Leftarrow 0$   3: $i\Leftarrow 1$   4: while$i\le p$  do   5:     $j\Leftarrow 1$   6:     while $j\le p$ do   7:         $u\Leftarrow 1$   8:         $t\Leftarrow 0$   9:         while $u\le q$ do 10:            $v\Leftarrow 1$ 11:            $s\Leftarrow 0$ 12:            while $v\le q$ do 13:                $s\Leftarrow s+yx(v,i)\times (kernel(u,v,delta)-\overline{L}(u))$ 14:                $v\Leftarrow v+1$ 15:            end while 16:            $t\Leftarrow t+yx(u,j)\times s$ 17:            $u\Leftarrow u+1$ 18:         end while 19:         $value\leftarrow value+\widehat{K}(i,j)\times t$ 20:         $j\Leftarrow j+1$ 21:     end while 22:     $i\Leftarrow i+1$ 23: end while 24: $value\Leftarrow value/((m-1)\ast (m-1))$

Algorithm 5 Calculate $\overline{M}=RowAverage(p,kernel,delta,m,y)$

1: $\overline{M}\Leftarrow zeros(1,p)$   2: $i\Leftarrow 1$   3: while  $i\le p$  do   4:     $j\Leftarrow 1$   5:     while $j\le p$ do   6:         $\overline{M}(i)\Leftarrow \overline{M}(i)+y(j)\ast kerenl(i,j,delta)$   7:         $j\Leftarrow j+1$   8:     end while   9:     $\overline{M}(i)\Leftarrow \overline{M}(i)/m$ 10:     $i\Leftarrow i+1$ 11: end while

4. Experiments

We employed representation of the data in a matrix $X_{i_j}^{pq}\in {\{0,1,2\}}^{m_{i_j}^{pq}\times 1}$ to calculate $\mathrm{HSIC}(X_{i_1}^{pq},X_{i_2}^{pq})$ by using a Gaussian kernel (${\sigma }^2=0.1$). In addition, we mapped $X_{i_1{i}_2\dots ,i_k}\in {\{0,1,2\}}^{m\times k}$ onto $X_{i_1{i}_2\dots ,i_k}\in {\{0,1\}}^{m\times 3k}$ to compute $\mathrm{HSIC}(X_{i_1{i}_2\dots ,i_k},Y_{i_1{i}_2\dots ,i_k})$ by also using a Gaussian kernel (${\sigma }^2=1$), i.e., $0\mapsto (1,0,0)$, $1\mapsto (0,1,0)$ and $2\mapsto (0,0,1)$. The advantages and disadvantages of the two representations have been explained by the other authors [29].

4.1. Evaluation Criterion

The evaluation criterion that we adopted in the experiments is by [9]:

$$\begin{array}{cc}Power& =\frac{♯S}{♯T},\end{array}$$

(10)

where $♯S$ is the number of found disease-causing SNP combinations (the epistasis SNPs score the highest) and $♯T$ is the number of datasets. Each dataset includes one disease-causing SNP combination. $Power$ is a measure of the accuracies of scoring criteria from genome data.

4.2. Simulated Datasets

For any data set, the worst-case scenario for checking the correctness of the scoring criteria is extensive testing of all SNP combinations. It is too computationally expensive for k = 4 and k = 5 cases. Therefore, tests were only conducted for k = 2 and k = 3.

4.2.1. Disease Models with k = 2

For k = 2, we used thirty-five disease models without marginal effects (DNME1–35) and six disease models with marginal effects (DME1–6). The models were designed based on interaction structures with different diseases, MAFs, prevalence (p) and h (the parameter settings are described in the supplementary files). Each data set contains 1000 SNPs and includes pairs of interacting SNPs (M0P0 and M1P1) generated according to the disease model setting, while other SNPs are generated using MAFs uniformly selected in [0.05, 0.5). For each model, we generated two simulated 100 data sets using the software GAMETES2.1 [30] with sample sizes of 400 (200 controls and 200 cases) and with sample sizes of 800 (400 control and 400 cases) [31].

Disease Models without Marginal Effects

We divided all DNMEs into seven subgroups for analysis according to the different combined values of h and MAF (DNME1–5 MAF = 0.2, h = 0.2; DNME6–10 MAF = 0.4, h = 0.2; DNME11–15 MAF = 0.2, h = 0.1; DNME16–20 MAF = 0.4, h = 0.1; DNME21–25 MAF = 0.2, h = 0.05; DNME25–30 MAF = 0.4, h = 0.05; DNME1–5 MAF = 0.2, h = 0.025).

The analysis results of subgroups of DNME1–35, each of which has 400 samples, are shown in Figure 3:

• Except for Mi, using tests on DNME1–10, the accuracy of all scoring criteria is close to 100%;
• All criteria are not very accurate using tests on DNME21–25 and DNME31–35;
• Mi has an extremely poor accuracy on all subgroup tests;
• LR has the highest accuracy using tests on DNME11–15 and DNME16–20, close to 100%, but is only a little more accurate than Mi on DNME21–25, DNME26–30 and DNME31–35 tests;
• The accuracy rates of both ND-JE and G-test rank in the middle overall, but G-test has the highest accuracy on the DNME26–30 test;
• The accuracy rate of K2-Score ranks second on DNME11–15 and DNME21–25 tests, third on the DNME26–30 test and slightly worse than ND-JE, HS${}_{CR}$ and G-test on the DNME16–20 test, but is only a little more accurate than Mi on the most difficult model (DNME31–35) test;
• HSIC${}_{CR}$ has the highest accuracy on the two most difficult model subgroup (DNME21–25 and DNME31–35) tests, especially on DNME31–35, where the accuracy is much higher than other criteria, and the overall accuracy on other model subgroups tests is similar to the other four criteria.

When the size of samples increased from 400 to 800, the accuracy of all criteria was greatly improved. The analysis results of subgroups of DNME1–35, each of which has 800 samples, are shown in Figure 4:

• Except for Mi, the accuracy of all criteria is close to 100% excluding tests on the two most difficult model subgroups (DNME21–25 and DNME31–35);
• Although the accuracy rate of Mi can be significantly improved with the increase in the size of samples, it is still relatively poor overall;
• With the number of samples increasing, there is still no change in the overall ranking, but the accuracy of the K2-Score on the DNME31–35 test rises to second;
• HSIC${}_{CR}$ has the highest accuracy on the two most difficult model subgroup tests.

Table 1. The number of times, out of 3500 data sets generated by 35 models without marginal effects, where $k=2$, that each scoring criterion identified epistasis SNPs of snp1000 for sample sizes of 400 and 800. The fourth column gives the total accuracy over all sample sizes. The last column gives the accuracy over all sample sizes in the most difficult subgroup models. The scoring criteria are listed in descending order of total accuracy.

Scoring Criterion	400 Samples	800 Samples	Total (%)	DNME31–35 (%)
HSIC${}_{CR}$	2535	3220	5755 (82.2%)	445 (44.5%)
K2-Score	2479	3186	5665 (80.9%)	336 (33.6%)
G-test	2443	3169	5612 (80.2%)	314 (31.4%)
ND-JE	2440	3163	5603 (80.0%)	301 (30.1%)
LR	2437	3158	5595 (79.9%)	297 (29.7%)
Mi	494	1971	2465 (35.2%)	192 (19.2%)

reveals the total average accuracy. From Table 1, we can find that Mi has a poor average accuracy; HSIC${}_{CR}$ has the best average accuracy regardless of the model’s sample scale of 400 or 800; although HSIC${}_{CR}$ is only slightly higher than the other four criteria on the total average accuracy, and the average accuracy on the most difficult model subgroup test is much better than other criteria.

Disease Models with Marginal Effects

We tested six DMEs for analysis according to MAF = 0.1 and the different combined values of heritability and prevalence (DME1 h = 0.031 and p = 0.050; DME2 h = 0.014 and p = 0.050; DME3 h = 0.01 and p = 0.050; DME4 h = 0.016 and p = 0.046; DME5 h = 0.009 and p = 0.026; DME6 h = 0.008 and p = 0.017).

The analysis results of DME1–6, each of which has 400 samples, are shown in Figure 5:

• The accuracy of all scoring criteria is close to 100% tested on DME1, except for Mi;
• Mi has extremely poor accuracy on all six models tests;
• Except for Mi, the accuracy rate of LR is worse than the other four criteria on DME2–6 tests, except that the accuracy on the DME3 test is nearly the same as that of HSIC${}_{CR}$;
• The accuracy rate of ND-JE ranks third on DME1 and DME3 tests, and fourth on the other four models tests;
• The accuracy rate of G-test ranks first on the DME1 test, second on DME2 and DME3 tests and third on the other four models tests;
• HSIC${}_{CR}$ has the highest accuracy rate on DME4 and DME6 tests, its accuracy rate on the DME5 test is slightly worse than LR and the accuracy rate on DME1–2 tests ranks third, whereas the accuracy rate on the DME3 test is a little better than Mi;
• K2-Score has the highest accuracy rate on DME1–3 and DME5 tests, its accuracy rate on DME4 and DME6 tests ranks second and it significantly outperforms the others on the most difficult model (DME3) test (although its accuracy rate in DME3 is below 50%).

When the size of samples increased from 400 to 800, the accuracy of all criteria was greatly improved. The analysis results of DME1–6, each of which has 800 samples, are shown in Figure 6:

• Although the accuracy of Mi can be significantly improved with the increase in the size of samples, it is still relatively poor overall;
• The accuracy rates of the other five scoring criteria all exceed 95% tested by the models, except on DME3;
• K2-Score has the highest accuracy rate on the most difficult model test (over 90%), the accuracy rate of G-test ranks second (over 80 %) and the accuracy rates of ND-JE, HSIC${}_{CR}$ and LR are not good enough, at just over $70\%$.

Table 2. The number of times, out of 600 data sets generated by six models with marginal effects, where $k=2$, that each scoring criterion identified epistasis SNPs of snp1000 for sample sizes of 400 and 800. The fourth column gives the total accuracy over all sample sizes. The last column gives the accuracy over all sample sizes in the most difficult model. The scoring criteria are listed in descending order of total accuracy.

Scoring Criterion	400 Samples	800 Samples	Total (%)	DME3 (%)
K2-Score	419	586	1005 (83.8%)	159 (79.5%)
HSIC${}_{CR}$	379	570	949 (79.1%)	86 (43%)
G-test	353	578	931 (77.6%)	108 (54%)
ND-JE	286	567	853 (71.1%)	92 (46%)
LR	273	559	832 (69.3%)	86 (43%)
Mi	55	349	404 (33.7%)	48 (24%)

reveals the total average accuracy. From Table 2, we can find that Mi has a poor average accuracy; the K2-Score has the best average accuracy rate regardless of the model’s sample scale of 400 or 800, where the main reason is that its accuracy rate on the DME3 test is much better than the other five scoring criteria; although the average accuracy of HSIC${}_{CR}$ tested on DME3 is not good enough, it ranks second in the overall average accuracy rate.

4.2.2. Disease Models with k = 3

For $k=3$, the data sets are generated by eight third-order epistasis pathogenic models (DM1–8), which are modeled by GAMETES2.1 according to the combinations of different MAFs ([0.2, 0.4]) and different heritability ([0.025, 0.05, 0.1, 0.2]) (DM1 MAF = 0.2, h = 0.025; DM2 MAF = 0.2, h = 0.05; DM3 MAF = 0.2, h = 0.1; DM4 MAF = 0.2, h = 0.2; DM5 MAF = 0.4, h = 0.025; DM6 MAF = 0.4, h = 0.05; DM7 MAF = 0.4, h = 0.1; DM8 MAF = 0.4 h = 0.2). The ♯quantiles of each combination is five. Every quantile of each pathogenic model corresponds to 100 simulated data files. Each file contains 100 SNPs and 1600 samples (800 normal, 800 diseased), and includes three interacting SNPs (M0P0, M1P1 and M2P2) generated according to the disease model settings, while other SNPs were generated using MAFs uniformly selected in [0.05, 0.5]. Therefore, the total number of the data sets is 4000 [6]. Detailed parameter settings are described in the supplementary file.

The analysis results of DM1–8, each of which has 1600 samples, are shown in Figure 7:

• The accuracy of all scoring criteria is close to 100% tested on DM2, DM3–4 and DM7–8;
• The accuracy of all scoring criteria is close to 80% on the DM1 test;
• The K2-Score has a poor accuracy rate on the most difficult model (DM5) test, whose accuracy is just close to $10\%$;
• The accuracy rates of the scoring criteria on the DM6 test are good enough, and the accuracy rates of the scoring criteria are close to 100%, except the K2-Score;
• HSIC${}_{CR}$ has the highest accuracy rate on the DM5 test, and its accuracy rate is the only one that exceeds 60% among all scoring criteria.

Table 3. The number of times, out of 4000 data sets generated by eight models, where $k=3$, that each scoring criterion identified epistasis SNPs of snp100 for 1600 samples. The second column gives the total accuracy over a sample size of 1600. The last column gives the accuracy over a sample size of 1600 in the most difficult model. The scoring criteria are listed in descending order of total accuracy.

Scoring Criterion	Total (%)	DM5 (%)
HSIC${}_{CR}$	3710 (92.8%)	316 (63.2%)
LR	3702 (92.6%)	300 (60%)
ND-JE	3700 (92.5%)	290 (58%)
Mi	3696 (92.4%)	289 (57.8%)
G-test	3677 (91.9%)	257 (51.4%)
K2-Score	3402 (85.1%)	63 (12.6%)

reveals the total average accuracy. From Table 3, we can find that the K2-Score has a poor average accuracy on the most difficult model test; the total average accuracy rates for all scoring criteria are good enough, with all criteria except the K2-Score achieving over 90% accuracy. Although HSIC has the best total average accuracy, its accuracy is not significantly better than other four criteria; however, tested on the most difficult model, the accuracy rate of HSIC${}_{CR}$ outperforms LR by 3.2%, and significantly outperforms the other four criteria, especially the K2-Score.

4.2.3. The Running Time Analysis

To demonstrate that our proposed method can be used as a lightweight scoring criterion, we proved in the previous section that its time complexity is $O(m)$. Furthermore, we calculated the average running time per dataset (unit in seconds) for the two-order with a sample size of 800 and the three-order in the simulation experiments, and we found that the average running time per dataset of our proposed method is between the other five lightweight methods (see

Table 4. Average running time (s) for the six scoring criteria per dataset for both the two-order tests and the three-order tests in the simulation experiments.

Scoring Criterion	2-Order (s)	3-Order (s)
HSIC${}_{CR}$	120.2426	80.3045
LR	77.3055	55.4438
ND-JE	134.922	76.5377
Mi	78.3	55.6138
G-test	77.6589	56.5221
K2-Score	125.0655	81.4237

). This further demonstrates the applicability of our proposed method as a lightweight scoring criterion.

In the experiment, all scoring criteria were implemented based on Matlab, and all tests were run on the environment of Windows 10 64 desktop computer with 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80 GHz, and 16.0 GB memory.

4.3. Case Study: A Real Chronic Dialysis Data

A real data set of 193 cases and 704 controls was selected from the mitochondrial D-loop region of chronic dialysis patients who were observed in a study by other authors [32]. The genotypes and locations of 77 SNPs are presented in

Table 5. Positions of chronic dialysis-associated 77 SNPS in mitochondrial d-loop region. ${}^a$ Left and right letters are major and minor genotypes, respectively. The number is the SNP position in the mitochondrial D-loop region.

SNP		D-Loop Position
1∼5	A16051G ${}^a$	T16086C	T16092M	T16093C	C16108T
6∼10	C16111T	T16126C	G16129A	T16136C	T16140C
11∼15	G16145A	C16148T	T16157C	A16162G	A16164G
16∼20	C16167T	T16172C	T16209C	T16217C	C16218T
21∼25	T16223C	A16227G	C16234T	A16235G	T16243M
26∼30	C16248T	T16249C	C16256T	C16257W	C16260T
31∼35	C16261T	C16266D	A16272G	G16274A	C16278T
36∼40	C16290T	C16291T	C16295T	C16297T	C16298T
41∼45	C16304T	A16309G	T16311C	A16316G	G16319A
46∼50	T16324C	C16327T	A16335G	C16355T	T16356C
51∼55	T16357C	T16362C	G16390A	A16399G	A16463G
56∼60	C16519T	A93G	G103A	T146M	C150T
61∼65	C151T	T152C	A153G	G185A	A189G
66∼70	C194T	T195C	T199C	A200G	T204C
71∼75	G207A	A210G	T217C	A234G	A235G
76∼77	T317C	C461T

[33].

The 77 SNPs contained in the subset of the chronic dialysis data set were used in the case study, which aims to give our readers more specificity regarding our proposed scoring criterion. First, for this dataset, we performed a full-space two-order SNP combination detection, meaning that the HSIC${}_{CR}$ values were evaluated for 2926 ($C_{77}^2$) possible combinations. Then, we selected the top ten HISC${}_{CR}$-valued combinations as candidate two-order epistasis SNP combinations to be raised for medical researchers, and the 10 candidate combinations are presented in

Table 6. The top ten highest HISC${}_{CR}$-valued two-order SNPs were used as ten candidate combinations.

Rank	Combination	HSIC${}_{\mathbf{CR}}$
1	41,   21	0.035922
2	52,   21	0.033105
3	41,   17	0.019069
4	56,   21	0.018961
5	68,   39	0.018545
6	21,   19	0.017254
7	60,   21	0.014506
8	17,   8	0.011645
9	17,   14	0.0097405
10	75,   36	0.0095467

.

5. Conclusions

In this paper, we verified with rigorous mathematical proof that HSIC${}_{CR}$ can be computed in $O(m)$ time. Moreover, we compared HSIC${}_{CR}$ with five representative scoring criteria for 49 simulation disease models. The experimental results show that: Mi has a poor accuracy on two-order disease models; the K2-Score has a poor accuracy on three-order difficult disease models; the accuracy rates of LR are not good enough on two-order disease models tests; HSIC${}_{CR}$, G-test and ND-JE have a high accuracy on all three classes of disease models tests; the accuracy rates of HSIC${}_{CR}$ rank first on two-order disease models without marginal effects tests and three-order disease model tests, and rank second on two-order disease models with marginal effects tests, although its advantage is not significant.

The advantages of HSIC${}_{CR}$ are: the methodology used is different from other scoring criteria, which makes it more complementary to other scoring criteria; it has a high accuracy on most disease models.

In the future, we will further investigate proposing efficient SIO algorithms to solve this problem by combining HSIC${}_{CR}$ and other effective lightweight criteria that already exist as weighted single or multi-objective functions. In addition, we will work with several local medical research institutions to use their real disease case-control study data to mine for disease-related SNP combinations by using our proposed approach. This will ultimately provide new guidance for drug development in complex diseases.