FORESEE: Fully Outsourced secuRe gEnome Study basEd on homomorphic Encryption

Abstract

Background

The increasing availability of genome data motivates massive research studies in personalized treatment and precision medicine. Public cloud services provide a flexible way to mitigate the storage and computation burden in conducting genome-wide association studies (GWAS). However, data privacy has been widely concerned when sharing the sensitive information in a cloud environment.

Methods

We presented a novel framework (FORESEE: Fully Outsourced secuRe gEnome Study basEd on homomorphic Encryption) to fully outsource GWAS (i.e., chi-square statistic computation) using homomorphic encryption. The proposed framework enables secure divisions over encrypted data. We introduced two division protocols (i.e., secure errorless division and secure approximation division) with a trade-off between complexity and accuracy in computing chi-square statistics.

Results

The proposed framework was evaluated for the task of chi-square statistic computation with two case-control datasets from the 2015 iDASH genome privacy protection challenge. Experimental results show that the performance of FORESEE can be significantly improved through algorithmic optimization and parallel computation. Remarkably, the secure approximation division provides significant performance gain, but without missing any significance SNPs in the chi-square association test using the aforementioned datasets.

Conclusions

Unlike many existing HME based studies, in which final results need to be computed by the data owner due to the lack of the secure division operation, the proposed FORESEE framework support complete outsourcing to the cloud and output the final encrypted chi-square statistics.

Introduction

Owing to the community effort on big data, biomedical science moves focus towards data-driven methodologies [1], which rely on collecting, integrating and analyzing large scale data. For biomedical studies, especially the genome analysis, the required storage and computational capacities may easily exceed the available resources in a single institution. Recently, cloud computing [2] emerges as a flexible alternative to support cost-effective biomedical research with big data. Researchers can rely on a cloud environment to easily scale up their studies with large scale data. However, the adopt of cloud computing in biomedical studies also yields more and more concerns about the potential data privacy risk in comparison with the local computing environment. As genome data are extremely sensitive, the storage of raw genome in a cloud may increase the disclosure risk.

The recently announced NIH policy [3] allows NIH funded studies to utilize public clouds to facilitate data analysis. However, the researchers instead of the cloud providers are responsible for the data security and privacy. Many existing attacks [4–6] also demonstrate the vulnerability of de-identified genome data. Thus, it is important to protect the privacy of genome data [7–9]. The rapid improvements of the data protection techniques make it possible to perform certain computations over encrypted data [10, 11] based on homomorphic encryption.

In [12], Gentry proposed the first fully homomorphic encryption scheme to enable both addition and multiplication operations over encrypted data. Brakerski et al. [13, 14] improved homomorphic encryption scheme based on learning with errors (LWE). Lauter et al. [15] presented several secure statistical algorithms for genetic association studies based on homomorphic encryption. Besides, Togan et al. [16] studied the integer comparison problem over homomorphic encrypted data. Recently, Graepel et al. [17] and Naehrig et al. [18] also showed that certain machine learning algorithms can be implemented using HME. Wang et al. [24] proposed a novel homomorphic encryption based framework to securely computing on exact logistic regression. Cheon et al. [19] developed a protocol for HME-based edit distance calculation that employed the greedy algorithm to obtain the upper bound of exact edit distance. Zhang et al. [25] improved homomorphic edit distance computation by combining path-finding algorithm and integer comparison.

In this paper, we propose the FORESEE framework to achieve secured and fully outsourced chi-square statistics computation in a public cloud. We assume that the cloud faithfully follows the protocol but may be curious of information from the received data, which is the so-called semi-honest adversary model [20]. The proposed FORESEE framework enables secure division operation over the homomorphic encrypted data and allows the cloud to directly release the study results. To be concrete, the contribution of this paper is two-fold.

• We develop a secure errorless division protocol, where a one-to-one mapping function is constructed for the floating numbers in computation and the study results can be accurately decrypted with a lookup table.

• We present a secure approximation division protocol to balance the complexity and accuracy with well-designed secure integer division in secure computation. In implementation, binary tree product and group-based computation are adopted to reduce circuit depth and the number of homomorphic multiplications.

For validation, experimental results show that the proposed FORESEE framework can identify all the significant SNPs based on the chi-square statistics with a moderate complexity using multiple slots for parallel computation.

Method

For clarity, in the rest of this paper, we use bold symbols to represent vector and matrix variables and normal symbols for scalar variables. Without specification, $\stackrel{⌢}{\Delta }$ is reserved for the encrypted version of variable or function Δ and log (·) stands for the logarithm with base 2.

Secure outsourcing GWAS

In this paper, we focus on the task of secure outsourcing GWAS in the 2015 iDASH challenge [21]. Given the genotypes from two groups over a number of single nucleotide polymorphisms (SNPs), we aim to securely calculate the chi-square statistics for the SNPs between the given case-control groups. The chi-square statistic) χ² is used by chi-square test to statistically assess whether there is significant association between the genetic variants and disease status. Typically, χ² is obtained by cumulating the normalized squared deviations between the observed and expected frequency distribution of alleles.

$${\chi }^2={\displaystyle \sum _i}{\displaystyle \sum _j}\frac{{\left(O_{i,j}-E_{i,j}\right)}^2}{E_{i,j}}$$

(1)

Here, O _i,j and E _i,j are the observed and expected allele counts for allele j, e.g. j = 1 for allele 'A' and j = 2 for allele 'a' in (see Table 1) from the case (i = 1) or control (i = 2) group, respectively.

Table 1 Observed allele counts for SNP, where O_1,1 and O_1,2 are the number of alleles A and a in the case group, O_2,1 and O_2,1 are the corresponding counts in the control group, N₁ and N₂ are the total allele counts for the case and control group, respectively.

Let us denote $N_1=O_{1,1}+O_{1,2}$ and $N_2=O_{2,1}+O_{2,2}$ the total number of alleles in the case and control groups, respectively. In general, E _i,j is computed by $\left(\left(O_{1,j}+O_{2,j}\right)\cdot N_i\right)/\left(N_1+N_2\right)$ for i = 1, 2 and j = 1, 2. If we assume that the case-control groups have the same number of n patients, we can obtain $N_1=N_2=2n$. Thus, Equation (1) can be simplified by

$${\chi }^2=\frac{4n\cdot {\left(O_{1,1}-O_{2,1}\right)}^2}{\left(O_{1,1}+O_{2,1}\right)\left[4n-\left(O_{1,1}+O_{2,1}\right)\right]}$$

(2)

Equation (2) indicates that, in addition to homomorphic additions and multiplications, the χ² statistic computation over encrypted dataset requires one secure division for fully outsourced GWAS, which is not supported in many existing HME-based schemes [15, 17, 22]. For example, if the numerator and denominator in Equation (2) are released directly due to the lack of secure division operation, one can easily infer the underlying allele counts (i.e., O_1,1and O_2,1) by solving a system of equations. To address the problem, we propose the FORESEE framework to enable secure division operation for the χ² statistic computation on an untrusted cloud.

The proposed framework

Figure 1 illustrates the proposed FORESEE framework, which allows secured and fully outsourced chi-square statistics computation in a public cloud and enable flexible release of study results. Using homomorphic encryption, the data owner can encrypt observed allele counts and directly upload to the public cloud. Consequently, the chi-square statistics can be securely computed according to Equation (2) based on homomorphic computation. Contrary to many existing HME-based schemes [15, 17, 22], the proposed framework develops two protocols for secure division operations over encrypted data, so that the final results are not necessarily computed by the data owner. As a result, authorized users are able to access the encrypted study results when granted the private key for decryption. Remarkably, the secrecy of uploaded sensitive information and released study results can be guaranteed under the proposed framework, as the trusted party would not interact with the untrusted public cloud. Thus, the proposed scheme enables secure outsourcing of the chi-square statistic computation to public cloud services, by which individuals or single institutions could contribute to the chi-square statistic computation in GWAS in a secure manner.

Figure 1

Conceptual diagram for the proposed FORESEE framework.

In the FORESEE framework, we develop two protocols for secure division operations, namely, secure errorless division and secure approximation division. The secure errorless protocol makes a secure one-to-one mapping from floating numbers to a set of encrypted positive integers. Consequently, authorized users can decrypt the study results with a lookup table. To achieve errorless division, the proposed protocol requires a deep circuit.

To balance the accuracy and complexity in chi-square statistic computation, the secure approximation division protocol is proposed as an alternative solution. Using secure integer division, the protocol approximates the study results with a tunable error rate. To improve its efficiency, binary tree product and group-based computation are designed to reduce circuit depth and the number of homomorphic multiplications.

In the following subsections, we will elaborate both protocols developed for the FORESEE framework.

Secure errorless division protocol

In this section, we propose the secure errorless division protocol when both dividend and divisor are small (e.g., less than 100). Considering that secure division operation is not available in existing HME-based schemes [15, 17, 22], we construct a one-to-one mapping function from floating numbers to a set of encrypted positive integers. Thus, the study results can be accurately decrypted with a lookup table corresponding to the one-on-one mapping function.

Secure mapping for division outcomes

To map the study result (in floating numbers), we construct a function with an integer output that uniquely corresponds to the division outcomes given a dividend and divisor. Let us denote $m\in \left[0,\bar{m}\right]$ and $w\in \left[1,\bar{w}\right]$ the dividend and divisor, respectively. Here, the upper bounds $\bar{m}$ and $\bar{w}$ of m and w should be predefined, so that the lookup table for decryption can be synchronized for all the authorized users. Consequently, we construct a two-dimensional function $\mathcal{F}\left(m,w\right)$ that returns the positive integer u _m,w corresponding to an index of the division result of m/w in floating number.

$$\mathcal{F}\left(m,w\right)=u_{m,w}$$

(3)

In the ciphertext domain, u _m,w can be determined by the polynomials of m and w related to the ciphertext modulus p. According to the Fermat Theory, we can construct a simplified function with less number of homomorphic multiplications. Given the prime p > mw , the secure mapping function is

$$\widehat{\mathcal{F}}\left(\widehat{m},ŵ\right)\equiv \widehat{m}{ŵ}^{p-2}\left(modp\right).$$

(4)

In Proposition 1, we demonstrated that the secure mapping proposed in Equation (4) is a one-to-one mapping from floating outcomes of m/w to a set of encrypted positive integers.

Proposition 1 Given arbitrary positive integers m₁, m₂, w₁, and w₂ taking their values $\left[1,\left|\sqrt{p}\right|\right]$, they satisfy

$$\frac{m_1}{w_1}=\frac{m_2}{w_2},$$

(5)

if and only if $\mathcal{F}\left(m_1,w_1\right)\equiv \mathcal{F}\left(m_2,w_2\right)\left(modp\right)$, where $\left|\sqrt{p}\right|$ is the round function that returns the maximum integer not greater than $\sqrt{p}$.

Proof. Please refer to Appendix I.

Proposition 1 implies that Equation (4) can map any pairs of $\left(\widehat{m},ŵ\right)$ with the same irreducible fraction to the same outcome $\widehat{\mathcal{F}}\left(\widehat{m}*,ŵ*\right)$, where m^* and w^* are the integer numerator and denominator, respectively that have no other common divisors. For example, given the ciphertext modulus $p=101,\widehat{\mathcal{F}}\left(\widehat{2},\widehat{1}\right)$ would be $\widehat{2}$ for the pairs $\left(\widehat{2},\widehat{1}\right),\left(\widehat{4},\widehat{2}\right)$ and $\left(\widehat{8},\widehat{4}\right)$. This fact means that the encrypted outcome can be securely released, as the authorized users can only obtain the accurate irreducible fraction m^*/w^*, but cannot infer the exact value of $\left(\widehat{m},ŵ\right)$.

Algorithm 1: Secure errorless division

0: Inputs: encrypted variable $\widehat{m},ŵ$, upper bound $\bar{m},\bar{w}$, the ciphertext modulus p.

1: Let ${{ŝ}_0}^{*}=ŵ,{ŝ}^{*}=\widehat{m}$.

2: Let $u*=⌊\mathsf{\text{log}}\left(p-2\right)⌋$

3: Decompose p − 2 as $p-2={\Sigma }_{i=0}^{\mathcal{h}}{2}^{{\mathcal{v}}_i}$, where $\mathcal{h }$ is the number of nonzero bits in the binary representation of p − 2 and ${\mathcal{v}}_i$. is the position of i-th nonzero bit.

4: For each i = 1,2,⋯, u^*

5:   ${ŝ}_i^{*}={ŝ}_{i-1}^{*}*{ŝ}_{i-1}^{*}$

6: end for

7: For each i = 0,1,⋯,$\mathcal{h }$

8:   ${\widehat{s}}^{*}={\widehat{s}}^{*}*{ŝ}_{{\mathcal{v}}_i}^{*}$

9: end for

10: Outputs: ${ŝ}^{*}$

During decryption, users can find the accurate study result $\mathcal{ r}$ with a lookup table, which consists of all possible irreducible fractions within ranges $m\in \left[0,\bar{m}\right]$ and $w\in \left[1,\bar{w}\right]$. Here, we provide two examples, where $\bar{m}=\bar{w}=10$ and p is set to 101 as the smallest prime greater than $\bar{m}\bar{w}=100$. It is worth mentioning that we can obtain the study result $\mathcal{ r}$ in floating number in Example 2. This fact verifies the accuracy of the proposed secure errorless division.

Example 1 The authorized users would obtain s^* = 2 by decrypting ${ŝ}^{*}=\widehat{2}$. The pair of co-prime integers (m, w) corresponding to $\widehat{\mathcal{F}}\left(\widehat{m},ŵ\right)\equiv \widehat{2}$ (mod 101) is (2,1). Thus, $\mathcal{ r}$ = m/w = 2.

Example 2 When ${ŝ}^{*}=\widehat{35},\left(m,w\right)=\left(4,3\right)$ as $\widehat{\mathcal{F}}\left(\widehat{m},ŵ\right)=\widehat{4}\cdot {\widehat{3}}^{99}\equiv \widehat{35}\left(mod101\right)$. As a result, $\mathcal{ r}$ = 4/3.

Secure approximation division protocol

In this subsection, we aim to develop the secure approximation division protocol. Since n (i.e., the number of patients in case or control group) is assumed to be a known integer, we denote A and B the dividend and divider of $\frac{{\left(O_{1,1}-O_{2,1}\right)}^2}{\left(O_{1,1}+O_{2,1}\right)\left[4n-\left(O_{1,1}+O_{2,1}\right)\right]}$ in Equation (2), respectively. Thus, the chi-square statistic can be rewritten as

$${\chi }^2=4n*A*\left(\frac{1}{B}\right)$$

(6)

where $A={\left(O_{1,1}-O_{2,1}\right)}^2$ is a nonnegative integer, and $B=\left(O_{1,1}+O_{2,1}\right)\left[4n-\left(O_{1,1}+O_{2,1}\right)\right]$ is a positive integer. Thus, given encrypted counts ${\widehat{O}}_{1,1}$ and ${\widehat{O}}_{2,1}$, $\widehat{A}$ and $\widehat{B}$ can be obtained with homomorphic multiplications and additions. Since the fraction team 1/B, with the value less than one, cannot be evaluated in the ciphertext domain, we scale it up by multiplying a positive integer ℳ. Therefore, the χ² statistic can be approximated by

$${\chi }^2=\frac{4n\cdot decrypt\left(\widehat{A}⌊\frac{\hat{M}}{\widehat{B}}⌋\right)}{M},$$

(7)

where $⌊M/B_i⌋$ is the round function that returns the maximum integer not greater than M/B _i , e.g., $⌊7/3⌋=2$ and $⌊10/15⌋=0$. Here, ℳ is a public information and should be large enough, as the upper bound of relative error is determined by $1/\mathsf{\text{min}}\left(\frac{\mathcal{M}}{B}\right)\times 100\%=\left(\frac{400n^2}{\mathcal{M}}\right)\%$.

Usually, we set $\mathcal{M}=\mathsf{\text{min}}\left(p-1,⌊\frac{p-1}{\mathsf{\text{max}}\left(\frac{A}{B}\right)}⌋\right)$, where p is the ciphertext modulus. According to Equation (7), we develop the secure approximation division protocol based on secure integer division.

Secure integer division

In this subsection, we describe the secure integer division protocol to achieve secure To compute $⌊\frac{\hat{\mathcal{M}}}{B}⌋$ o in Equation (7), we first introduce a vector Twith its 6-th element defined by

$$t_i=⌊\frac{\mathcal{M}}{B_i}⌋i\in \left[1,2n\right]$$

(8)

where B _i = i * (4n− i), i = 1,2,...,2n are the possible values of B in chi-square statistic computation (see equation (6)). Consequently, given ℳ, we define a function f(x) which satisfies

$$f\left(B_i\right)=⌊\frac{\mathcal{M}}{B_i}⌋i\in \left[1,2n\right]$$

(9)

In our implementation, a one-dimensional function is formulated using Lagrange interpolating polynomial with x ∈ {B₁,..., B_2n}

$$f\left(x\right)={\displaystyle \sum _{i=1}^{2n}}{t}_i\frac{{\prod }_{\underset{l\ne i}{1\le l\le 2n,}}\left(x-B_l\right)}{{\prod }_{\underset{l\ne i}{1\le l\le 2n,}}\left(B_i-B_l\right)}$$

(10)

Since division is intractable for homomorphic encrypted data, we need to derive a surrogate function for Equation (10) that can be implemented based on homomorphic multiplications and additions. For simplicity, we denote u _i . the divisor for x = B _i . in Equation (10).

$$u_i={\displaystyle \prod _{1\le l\le 2n,l\ne i}}\left(B_i-B_l\right)$$

(11)

Consequently, we can construct a surrogate function for Equation (10) by numerically finding a set of integers v _i . with 1 ≤ v _i . ≤ p − 1 for 1 ≤ i ≤ 2n, that satisfy

$$u_i{v}_i\equiv 1\left(modp\right)$$

(12)

Here, p is the cipheretext modulus (i.e., a prime under double-CRT representation in the BGV scheme). Thus, we demonstrate the existence of {v _i } in Proposition 2 to guarantee the computational tractability of f(x) in the ciphertext domain.

Proposition 2 For each u _i = 1,2,⋯,2n, given p > ℳ, at least one v _i can be found to satisfy (12).

Proof. Please refer to Appendix II.

Substituting $1/{\Pi }_{1\le l\le 2n,l\ne i}\left(B_i-B_l\right)$ with v _i in Equation (10), we can reformulate f(x) with multiplications instead.

$$f\left(x\right)={\displaystyle \sum _{i=1}^{2n}}\left[t_i{v}_i{\displaystyle \prod _{\underset{l\ne i}{1\le i\le 2n,}}}\left(x-B_l\right)\right]$$

(13)

We transform Equation (13) into the combination of polynomials of x by expanding the products and combining the coefficients.

$$f\left(x\right)={\displaystyle \sum _{i=0}^{2n-1}}{h^{\prime }}_i{x}^i$$

(14)

Here, ${h^{\prime }}_i$ is the coefficient for the i-th order of x (i.e., xⁱ.) after polynomial expansion, which includes $v_i{B}_l$ and t _i . In the ciphertext domain, we can construct the function $\widehat{f}\left(\widehat{x}\right)$ for secure integer division $⌊\hat{\mathcal{M}}/\widehat{x}⌋$.

$$\widehat{f}\left(\widehat{x}\right)\equiv {\displaystyle \sum _{i=0}^{2n-1}}{ĥ}_i{\widehat{x}}^i\left(modp\right)$$

(15)

where $\widehat{x}\in \left\{{\widehat{B}}_1,{\widehat{B}}_2,\dots ,{\widehat{B}}_{2n}\right\}$ are finite positive encrypted integers, and ${ĥ}_i\in \left[\widehat{0},\widehat{p}-\widehat{1}\right]$ is obtained by encrypting $h_i\equiv {h^{\prime }}_i\left(modp\right)$. We set h _i = 0 with i > 2n − 1.

Implementation optimization

The secure integer division can be optimized to further reduce the cumulative circuit depths (CCD) and number of homomorphic multiplications (HMs). To achieve this goal, we adopt group-based computation and binary tree product to generate $\widehat{f}\left(\widehat{x}\right)$ in implementation.

To reduce the number of HMs, a group-based computation is adopted to calculate $\widehat{f}\left(\widehat{x}\right)$. The key idea of the proposed group-based optimization is to first compute a set of ${ĥ}_{c\cdot d+i}{\widehat{x}}^i$ with i ∈ [0, d], where d is number of elements in each group, and c = 0,..., C is the group index with the total number of groups $C=⌊\left(2n-1\right)/d⌋+1$. After grouping, we get the following equation with a reduced number of HMs.

$$\widehat{f}\left(\widehat{x}\right)\equiv {\displaystyle \sum _{c=0}^{C-1}}\left[{\widehat{x}}^{c.d}\left({\displaystyle \sum _{i=0}^{d-1}}{ĥ}_{c\cdot d+i}{\widehat{x}}^i\right)\right]\left(modp\right)$$

(16)

Algorithm 2 describes the generation of $\hat{X}=\left(\widehat{1},\widehat{x},\dots ,{\widehat{x}}^d\right)$ using binary tree product. The number of HMs and CCD required to calculate $\hat{X}$ can be reduced to d − 1 and $⌊\mathsf{\text{log}}\left(d-1\right)⌋+1$, respectively.

Algorithm 2: Binary tree product for generating $\hat{X}$

0: Inputs: encrypted variable $\widehat{x}$, the maximum power d

1: For i = 2,3,⋯, d

2:   Let $l_1=2^{⌊\mathsf{\text{log}}\left(i-1\right)⌋}$.

3:   Let $l_2=i-l_1.$

4:   ${\widehat{x}}^i={\widehat{x}}^{l_1}\cdot {\widehat{x}}^{l_2}.$

5: end for

6: Outputs: $\widehat{X}=\left(\widehat{1},\widehat{x},\dots ,{\widehat{x}}^d\right)$

An additional optimization can be applied in equation (16) by replacing the multiplication ${ĥ}_{c\cdot d+i}{\widehat{x}}^i$ as the summation over a total number of ${ĥ}_{c\cdot d+i}$ additions of ${\widehat{x}}^i$ to reduce the number of HMs.

Since the time cost of HMs is larger than HAs, we determine d by minimizing the number of HMs. As shown in Table 9 the total number of HMs required for secure integer division is 2C + d − 3. The number of groups C and the number of elements in each group d are selected to minimize the number of HMs $F\left(d\right)=d+2⌊\left(2n-1\right)/d⌋-3$. Given an integer $n,F\left(d\right)\approx d+\frac{2\left(2n-1\right)}{d}-3$ can obtain its minimum$2\sqrt{4n-2}-3$, only when $d=\sqrt{2\left(2n-1\right)}$. Since d is an integer, it is estimated by $⌊\sqrt{2\left(2n-1\right)}⌋$ to minimize F(d). Thus, C can be estimated by $⌊\left(2n-1\right)/d⌋+1$ accordingly. Using the optimal d and C, secure integer division can be achieved based on the encrypted function $\widehat{f}\left(\widehat{x}\right)$ in Equation (15). Algorithm 3 elaborates the secure integer division. In line 2, in order to obtain $\widehat{X^{\prime }}$, the inputs of Algorithm 2 are set to${\widehat{x}}^d$ and C − 1, respectively.

Algorithm 3: Secure integer division

0: Inputs: encrypted variable $\widehat{x}$, group size d , the number of groups C, the ciphertext modulus p, the polynomial parameters h _i , i = 0,1,...,2n − 1

1: Compute $\widehat{X}=\left(\widehat{1},\widehat{x},\dots ,{\widehat{x}}^d\right)$ according to Algorithm 2

2: Compute $\widehat{X}=\left(\widehat{1},{\widehat{x}}^d,\dots ,{\widehat{x}}^{\left(C-1\right)d}\right)$ according to Algorithm 2

3: For each c = 0,1,⋯, C − 1

4:   For each i = 0,1,⋯, d − 1

5:      Calculate ${\widehat{h}}_{cd+i}{\widehat{x}}^i$

6:   end for

7: end for

8: Let $\widehat{a}=\widehat{0}$

9: For each c = 0,1,⋯, C − 1

10:   ${\widehat{a}}^{\prime }=\widehat{0}.$

11:   For each i = 0,1,2,⋯, d − 1

12:      Update ${\widehat{a}}^{\prime }={\widehat{a}}^{\prime }+{ĥ}_{cd+i}{\widehat{x}}^i$

13:   end for

14:   Update $\widehat{a}=\widehat{a}+{\widehat{a}}^{\prime }{\widehat{x}}^{cd}$

15: end for

16: Outputs: $\widehat{a}=\widehat{f}\left(\widehat{x}\right).$.

Parallel computation using multiple slots

Since HME schemes with ciphertext space ${\mathbb{Z}}_q^{L_s}$ support single instruction multiple data (SIMD) with L _s slots, we can use parallel computation to reduce the number of homomorphic multiplications (HMs) and homomorphic additions (HAs). Denote $\widehat{a}=\left({\widehat{a}}_1,{\widehat{a}}_2,\dots ,{\widehat{a}}_{L_s}\right)$ and $\widehat{b}=\left({\widehat{b}}_1,{\widehat{b}}_2,\dots ,{\widehat{b}}_{L_s}\right)$ the two encrypted ciphertexts with L _s slots. SIMD is applicable to simultaneous computation of the addition $\widehat{a}+\widehat{b}=\left({\widehat{a}}_1+{\widehat{b}}_1,{\widehat{a}}_2+{\widehat{b}}_2,\dots ,{\widehat{a}}_{L_s}+{\widehat{b}}_{L_s}\right)$ and $\widehat{a}\cdot \widehat{b}=\left({\widehat{a}}_1\cdot {\widehat{b}}_1,{\widehat{a}}_2\cdot {\widehat{b}}_2,\dots ,{\widehat{a}}_{L_s}\cdot {\widehat{b}}_{L_s}\right)$ multiplication. In two ciphertexts, only two slots in the same position can operate with each other.

In Algorithm 1, multiple encrypted outputs can be calculated at the same time with parallel computation. When the result is returned back to the user, the user extracts the integer in each slot and search it in the lookup table. Noticeably, in the parallel computation, m^u and n^u should be selected as the upper bounds of all the dividends and divisors in the slots. Similarly, multiple slots can also be used in the secure approximation division protocol. The secure integer division developed in Algorithm 3 can be simultaneously conducted for L _s pairs of inputs $\left({\widehat{a}}_i,{\widehat{b}}_i\right),i=1,2,\dots ,L_s$ using multiple slots.

Results

In this section, we evaluate the proposed FORESEE framework, which was implemented with HElib [23], one of the most efficient open-source HME libraries based on the LWE theory [13, 14]. The evaluations were made on an Ubuntu 14.04 server with Intel Xeon CPU E5-2687W @ 3.10GHz and 256 GB memory. We present the performance in the terms of time and memory cost. First, we provide the results of secure errorless division on simulated data. Moreover, we provide the performance of chi-square statistics based on the secure approximation division.

Simulation study

Table 2 elaborates the experimental setups for the secure errorless division protocol. Given the ciphertext modulus p, the upper bound $\bar{m}$ (i.e., dividend) and $\bar{w}$ (i.e., divisor) is set to $⌊\sqrt{p}⌋$. A number of L _s slots are used for parallel computation. The lifting parameter for plaintext base is set to 1. The security level is 80. The number of columns in key switching is 2. Hamming distance is 64.

Table 2 Experimental setups for secure errorless divisio.

Using HElib, one ciphertext can contain multiple slots to have many integers encrypted into the ciphertext with the public key. Thus, the size of the public key is related to the number of multiple slots L _s in addition to the ciphertext modulus p and the number of levels in modulus chain L. Taking Table 2 for example, L _s for $\left(\bar{m},\bar{w}\right)=\left(70,70\right)$ is 3144, which is greater than most ones. Thus, its ciphertext sizes are much larger than the other configurations with close values of $\bar{m}$ and $\bar{w}$.

Using HElib, we are able to evaluate all the slots in the ciphertext in parallel. Table 3 shows the average execution time for the secure errorless division protocol. Based on the lookup table generated for various parameters $\left(\bar{m},\bar{w}\right)$, the proposed protocol is efficient for secure division operation over m ≤ 100 and w ≤ 100. However, its circuit depths increase rapidly with the growth of m and w, which limits its application for larger dividends and divisors.

Table 3 Time cost in seconds for key generation, encryption, and errorless division computation using various parameters.

Chi-square statistic computation

We employ the secure approximation division protocol in secure chi-square statistic computation. Two datasets from iDASH genome privacy protection challenge are used for evaluation, which contain 311 SNPs and 610 SNPs, respectively in the case-control groups, each consisted of with 200 individuals,

In homomorphic encryption, the ciphertext modulus p and the number of levels in modulus chain L are set to 25600000039 and 51, respectively. The public and private key sizes are both around 2.6 GB. The lifting parameter for plaintext base is set to 1. The security level is 80. The number of columns in key switching is 2 and the Hamming distance is 64. To reduce computational complexity, we use L _S = 864 slots in parallel computation. For secure integer division, ℳ is 25600000000. f and a are accordingly set to 28 and 15 for 200 individuals in each group.

Table 4 provides the time cost for homomorphic evaluation of both datasets in chi-square statistics computation, including key generation, encryption, total and average execution time. Using multiple slots, the secure approximation division protocol can achieve the chi-square statistics computation in less than one second in average. Table 5 evaluates the accuracy of computation in terms of the mean-squared error (MSE) and maximum error between the exact and the approximate chi-square statistics, where the MSE are less than 5 × 10⁻¹⁰ . The evaluation on maximum error also supports the conclusion.

Table 4 Time cost in seconds for key generation, encryption, and the computation of chiRsquare statistics using different parameters based on secure approximation division.

Table 5 Recall and precision with different p-value cutoffs in identifying significant SNPS on both datasets, where the mean-squared error (MSE) and maximum error between the exact and the approximate chi-square statistics

Remarkably, we also computed the p-value for each SNP based on the chi-square statistic and applied different p-value cutoffs as 0.05, 0.01, and 0.005. The secure approximation division protocol is demonstrated to find out all the significant SNPs for both datasets under different p -value cutoffs. As a result, the proposed protocol provides a good tradeoff between accuracy and complexity for secure chi-square statistic computation.

Furthermore, we compare the two proposed protocols in the chi-squared statistics computation. For secure errorless division protocol, we use the same parameters list above, except that L is set to 151 to guarantee the required circuit depth in implementation. Table 6 and 7 compare the computational complexity and storage cost for the two protocols, respectively. In Table 6 the secure errorless division protocol requires about 10, 20 and 5 times in complexity for key generation, encryption, and execution (computation), when compared with the secure approximation division protocol. Table 7 shows that the ciphertext key sizes for secure errorless division are about 8 times larger due to the greater L. These results imply that the secure approximation division protocol provides a good trade-off in terms of complexity and accuracy for chi-squared statistic computation.

Table 6 Time cost in seconds for key generation, encryption and the computation of chi8squared statistics using different parameters based on secure errorless division (ED) and secure approximation division (AD) protocols.

Table 7 L (i

Discussions

In this section, we analyze the computational complexity of the proposed FORESEEE protocol and discuss its potential extension and its limitation.

Complexity analysis

In this subsection, we make an analysis on the computational complexity of the secure errorless division and secure approximation division protocols. Cumulative circuit depth2 (CCD) and the numbers of homomorphic multiplications (HMs) are provided for both protocols in the FORESEE framework.

We begin with the complexity analysis for secure errorless division protocol (i.e., Algorithm 1). As shown in Table 8 the number of HMs to calculate ${{ŝ}_i}^{*}$ at each iteration in A1 line 5 is 1. The number of HMs to obtain s^* at each iteration in A1 lines 8 is also 1. Therefore, the CCD in calculating ${{ŝ}_i}^{*}$ in A1 lines 4-6 are $⌊\mathsf{\text{log}}\left(p-2\right)⌋$. The depths to obtain s^* in A1 lines 7-9 are $⌊\mathsf{\text{log}}\left(p-2\right)⌋+h+1$

Table 8 Complexity analysis in terms of cumulative circuit depth2 (CCD) and the number of homomorphic multiplications (HMs) for secure errorless division protocol (Algorithm 1).

Table 9 provides the CCD and number of HMs for secure approximate division (i.e., Algorithm 3). The number of HMs to obtain $\widehat{X}$ and ${\widehat{X}}^{\prime }$ are d − 1 and C − 2, respectively. To evaluate Equation (18), the total number of HMs are d + 2C − 3. By using binary tree product based optimization, the circuit depths required for computing $\widehat{X}$ and ${\widehat{X}}^{\prime }$ are $⌊\mathsf{\text{log}}\left(d-1\right)⌋+1$ and $⌊\mathsf{\text{log}}\left(c-2\right)⌋+1$,respectively. Finally, the total CCD for secure approximate division operation is $⌊\mathsf{\text{log}}\left(c-2\right)⌋+⌊\mathsf{\text{log}}\left(d-1\right)⌋+3$.

Table 9 Complexity analysis in terms of cumulative circuit depth2 (CCD) and the number of homomorphic multiplications (HMs) for secure approximation division (Algorithm 3).

Potential extension

In this paper, we proposed the FORESEE framework to address the problem of fully outsourcing chi-square statistic computation to a public cloud. However, the application scenarios for the FORESEE framework, especially the secure approximation division protocol, can be further extended to securely compute other statistics tests that involve division operations. One intuitive example is the Transmission disequilibrium test (TDT) developed to assess the genetic linkage between the genetic variants and disease status in family-based association studies. TDT is based on the binomial test with one degree of freedom, which is asymptotically equivalent to the chi-square hypothesis test.

Limitation

There are several limitations in the FORESEE framework. First, in secure approximation division protocol, the upper bound of approximation error depends on the ciphertext modulus G. Therefore, G should be large enough to guarantee the accuracy in computation, which degrades the efficiency of the FORESEE framework. Second, the computational and storage costs based on homomorphic encryption are still very high. For example, key generation and encryption is much more timeconsuming than computation. The ciphertext sizes are also a heavy burden for communication. Finally, it is still a challenging problem to generalize the secure errorless division protocol. In summary, there is still room to improve the proposed division protocols in the FORESEE framework through better algorithm design, efficient coding in the HElib and parallelization.

Conclusion

In this paper, we proposed a novel FORESEE framework for the secure outsourcing GWAS in the iDASH genome privacy protect challenge, especially for the chi-square statistic computation. The proposed framework consists of two protocols for secure division operation, namely secure errorless division and secure approximation division. The secure errorless protocol made a bijection between floating numbers and a set of encrypted positive integers. Thus, it could output the accurate study results based on a lookup table. On the other hand, the secure approximation division protocol adopted secure integer division to obtain approximate study results with a tunable accuracy. The protocol was able to balance the complexity and accuracy by using the group-based computation and binary tree product with improved efficiency. In comparison to existing HME-based schemes [15, 17, 22], both protocols enabled fully outsourced secure GWAS in an untrusted public cloud and could directly release study results to authorized users for decryption. Experimental results show that the secure approximation division protocol can capture all the significant SNPs in chi-square statistic computation with a moderate computational complexity.

Appendix I: Proof of Proposition 1

Since $m_1{w}_1^{p-2}\equiv m_2{w}_2^{p-2}\left(modp\right)$, we can obtain Equation (17) by multiplying w₁w₂ on the both sides.

$$w_2{m}_1{w}_1^{p-2}\equiv w_1{m}_2{w}_2^{p-2}\left(modp\right)$$

(17)

According to the Fermat's little theorem,

$$w_i^{p-1}\equiv 1\left(modp\right)i=1,2.$$

(18)

When w₁ and w₂ are coprime with p, we can find that

$$w_2{m}_1\equiv w_{1}m\left(modp\right)$$

(19)

Since $w_2{m}_1\le ⌊\sqrt{p}⌋*⌊\sqrt{p}⌋<p$ and $w_1{m}_2\le ⌊\sqrt{p}⌋*⌊\sqrt{p}⌋<p$, it holds for the prime p that

$$-p<w_2{m}_1-w_1{m}_2<p$$

(20)

From Equations (19) and (20), we obtain that $w_2{m}_1=w_1{m}_2$, which comes to Proposition 1.

Appendix II: Proof of Proposition 2

We recall the Fermat's little theorem that, given a prime p,

$$q^{p-1}\equiv 1\left(modp\right)$$

(21)

where p and q are coprime numbers. Since $u_i\in \left[1,p-1\right]$, the greatest common divisor for u _i and p is always 1. Given an integer u _i , we can derive ${v^{\prime }}_i=u_i^{p-2}$ to satisfy $u_i{v^{\prime }}_i\equiv 1$ (mod p) in Equation (12). Thus, by considering $v_i\equiv {v^{\prime }}_i$ (mod p), v _i . ∈ [1, p − 1] can be found for Equation (12). As a result, we draw the Proposition 2.