Efficient privacy-preserving string search and an application in genomics

Motivation: Personal genomes carry inherent privacy risks and protecting privacy poses major social and technological challenges. We consider the case where a user searches for genetic information (e.g. an allele) on a server that stores a large genomic database and aims to receive allele-associated information. The user would like to keep the query and result private and the server the database. Approach: We propose a novel approach that combines efficient string data structures such as the Burrows–Wheeler transform with cryptographic techniques based on additive homomorphic encryption. We assume that the sequence data is searchable in efficient iterative query operations over a large indexed dictionary, for instance, from large genome collections and employing the (positional) Burrows–Wheeler transform. We use a technique called oblivious transfer that is based on additive homomorphic encryption to conceal the sequence query and the genomic region of interest in positional queries. Results: We designed and implemented an efficient algorithm for searching sequences of SNPs in large genome databases. During search, the user can only identify the longest match while the server does not learn which sequence of SNPs the user queried. In an experiment based on 2184 aligned haploid genomes from the 1000 Genomes Project, our algorithm was able to perform typical queries within ≈ 4.6 s and ≈ 10.8 s for client and server side, respectively, on laptop computers. The presented algorithm is at least one order of magnitude faster than an exhaustive baseline algorithm. Availability and implementation: https://github.com/iskana/PBWT-sec and https://github.com/ratschlab/PBWT-sec. Contacts: shimizu-kana@aist.go.jp or Gunnar.Ratsch@ratschlab.org Supplementary information: Supplementary data are available at Bioinformatics online.


Introduction
String search is a fundamental task in the field of genome informatics, for which a large variety of techniques have been developed (see, for instance, Altschul et al., 1990;Kent, 2002;Li and Homer, 2010).
Traditionally, those techniques have been optimized for accuracy and computational efficiency, however a recent boom of personal genome sequencing and analyses has spotlighted a new criteria, namely, privacy protection. As reported in many studies, a genome is considered V C The Author 2016. Published by Oxford University Press.
to be one of the most critical pieces of information for an individual's privacy. In fact, it is largely different from any other personal information because it works as an identifier of an individual while it possesses the information that has strong correlation with the phenotype of the individual (Erlich and Narayanan, 2014;Roche and Annas, 2001). Therefore, in principle, privacy protection is an inevitable problem when handling personal genomes. As a practice, the most popular approach is protecting genomes physically; genomic sequences have been kept at few collaborator sites, and only a limited number of researchers are allowed to access them. This conservative approach severely limits the great potential of existing genomic resources. In order to mitigate the stagnation caused by privacy issues, it appears crucial to develop practical methods that enable searching and mining genomic databases in a privacy-preserving manner.
So far, several groups have tackled related problems. Jha et al. (2008) developed secure multi-party computation protocols for computing edit distance. Blanton and Aliasgari (2010) proposed a protocol to search DNA string against a DNA profile represented by finite automata. Bruekers et al. (2008) proposed a protocol to detect a match between two short DNA sequences for the purpose of genetic test. Baldi et al. (2011) also aimed for genetic test to develop a method for computing set intersection cardinality. Freedman et al. (2005) proposed a protocol for searching pre-defined keywords from databases. Naganuma et al. (2012) proposed a substring search protocol for public databases while keeping user's query private. Ayday et al. (2013) developed a system by using several cryptographic techniques to find a subset of short reads which includes a fixed-length query string at specific position. He et al. (2014) proposed an algorithm for finding relatives by secure identity-by-descent matches.
We propose a general approach which utilizes an efficient iteratively queriable data structure together with cryptographic techniques. Among many variations of such data structures, the Burrows-Wheeler Transform ( BWT Langmead et al., 2009;Li and Durbin, 2009;Li et al., 2009) and related techniques such as the positional BWT (PBWT; Durbin, 2014) have dramatically improved the speed of genomic database analyses. Those data structures commonly have an indexed dictionary called a rank dictionary. By referring to the rank dictionary in iterative operations, one can efficiently search the database. For the case of BWT, a match between query and database is reported as a left-open, right-closed interval ðf ; g, and the interval is computed by the look-up of the rank dictionary. In our approach, we access the rank dictionary in privacy-preserving manner by using additive homomorphic encryption and oblivious transfer (OT).
Cryptographic approaches often require significant computational resources. The goal of this work is to illustrate that privacypreserving queries are within reach when using current cryptographic techniques and standard computing hardware. We demonstrate that a typical query would only take about 4.6 s on the user side using a single thread and % 10.8 s on the server having four cores, while preserving privacy of the query string and the database.
The rest of the paper is organized as follows. In Approach, we describe the main ideas of our approach without going into technical details. In Methods, the detailed algorithm of recursive oblivious transfer is given followed by the description of a practical algorithm, named PBWT-sec, for privacy-preserving search in large-scale genotype databases. We also describe complexity and security properties of the proposed algorithm. We provide the more intricate details of a more efficient version of the algorithm in Supplementary Sections S1-S2. In Experiments, we evaluate the performance of PBWT-sec on datasets created from data of the 1000 Genomes Project (The 1000 Genomes Project Consortium, 2012) and compare it to an alternative method for fixed-length k-mer search. Finally, we conclude our study in Section 5.

Problem setup
We consider the setting in which a user would like to search a genomic sequence in a database with the aim to either determine whether this sequence exists in the queried database and/or to obtain additional information associated with the genomic sequence. An example is the use in a so-called genomic beacon (for instance, those created within the Beacon Project of the Global Alliance for Genome & Health (GA4GH).) Another application is the search of a specific combination of variants, for instance, in the BRCA1 or BRCA2 genes, with the aim to determine whether that combination of variants is known or predicted to be deleterious (see, for instance, GA4GH's BRCA Challenge). For privacy reasons, the user would like to conceal the queried sequence, which would be particularly relevant for the second example. For both examples it would be important that the server's database is protected.

Information flow of searches on recursive search data structures
Let us describe the information flow between a user and a server for such problems. In this work, we perform searches on the (positional) Burrows-Wheeler transform of a genomic database of length N. (P)BWT stores string information very efficiently and still allows computations (this is a property of Succinct Data Structures, see Jacobson, 1988).
To search for a query string q over the alphabet R, one iteratively operates on intervals that can later be used to identify the matching genomic regions based on the (P)BWT. A substring match is represented by an interval ðf ; g. The number of matches is given by the length of the interval g-f. It is known that the ðk þ 1Þth interval ðf kþ1 ; g kþ1 corresponding to a ðk þ 1Þ-mer match can be updated from the kth interval ðf k ; g k and the ðk þ 1Þth letter of the query q.
We will provide more details on how to update f and g in Section 3.3. To understand the key ideas, it is sufficient to understand that the updates can be written in the form of where c ¼ q½k þ 1 and v c 2 N N is a large, static lookup table. Hence, the iterative algorithm of updating ðf k ; g k by using the query q, can be written as a recursive algorithm: This can be done analogously for g kþ1 . In this work we will refer to data structures that can be queried in the recursive way described above as recursive search data structures. Figure 1 illustrates the information flow of a search on the recursive search data structure.

Oblivious transfer for privacy-preserving search
In a search on the recursive search data structures, the user needs to conceal not only a query string q but also f and g because f i is f iÀ1 th element of v q½i , and q½i is inferred from those two values. Analogously, q½i is also inferred from g i and g iÀ1 . The server needs to minimize output because the user reconstructs a part of v from the server's output. In this study, we achieve such security requirements by a cryptographic technique called oblivious transfer.

Oblivious transfer
Oblivious transfer (OT) is a cryptographic technique for two parties: the user and the server, and enables the user to specify 0 t < N and obtain only tth element of the server's vector v without leaking any information about t to the server (Rabin, 1981). Figure 2 illustrates an outline of the oblivious transfer. Among several efficient algorithms (Lipmaa, 2005(Lipmaa, , 2008Zhang et al., 2013), we used those which are based on additive homomorphic encryption. The detailed algorithm will be given in Section 3.2.

Concealing the query
The user's query consists of ðf i ; g i , and q½i þ 1 for ith iteration. A key idea of our approach is to look-up elements of v c by OT and obtain the next interval ðf iþ1 ¼ v c ½f i ; g iþ1 ¼ v c ½g i without revealing ðf i ; g i to the server. In our approach, we also use a masking technique such that the user tries v c ½f i for all c 2 R, and the server only returns v c ½f i where c ¼ q½i þ 1 without knowing the value of q½i þ 1. Technical details will be given in Section 3.2.

Concealing the database
While this approach protects a user's privacy, the server leaks information of v c which may be sufficient to reconstruct parts of the genotypes in the database. In order to rigorously protect the server's privacy, we propose a technique that allows for recursive oblivious transfer where the user does not learn intermediate results but only if a unique match was found. It is based on a bit-rotation technique which enables the server to returnf k :¼ Rðf k Þ andĝ k :¼ R 0 ðg k Þ which are random values to the user. Only the server can recover f k and g k in encrypted form (i.e. the server does not see f k and g k when recovering them), and thus the user can recursively access v c ½f k and v c ½g k correctly. The details of this approach are given in Section 3.2.
In this work, we designed an algorithm based on these techniques that can be used for privacy-preserving search in large genotype databases.
Note that there are still privacy risks for the server, though returning only a unique match minimizes the information leakage from the server. For example, assume there is a database storing a genomic study of drug addicts that implements the PBWT-sec, and a person (Bob) participated in the study. If someone obtains a DNA sample from Bob and queries the databases, he/she will reveal that Bob is a drug addict. As described in the above case, there is always a limitation for protecting the server's privacy as long as the server returns the search results, and there is associated information such as phenotypes (Shringarpure and Bustamante, 2015). We emphasize that this issue is common for any database search application and is not specific to our proposed method.

Additively homomorphic encryption
Our main cryptographic tool in this paper is an additive-homomorphic public-key encryption scheme ðKeyGen; Enc; DecÞ, which enables us to perform additive operations on encrypted values. Here, the algorithm KeyGen generates a public key pk and a secret key sk; EncðmÞ denotes a ciphertext obtained by encrypting message m under the given pk; and DecðcÞ denotes the decryption result of ciphertext c under the given sk. The scheme also has the following additive-homomorphic properties: • Given two ciphertexts Encðm 1 Þ and Encðm 2 Þ of integer messages m 1 and m 2 , Encðm 1 þ m 2 Þ can be computed without knowing m 1 , m 2 and the secret key (denoted by Encðm 1 ÞÈEncðm 2 Þ). • Given a ciphertext EncðmÞ of a message m and an integer e, Encð e Á mÞ can be computed without knowing m and the secret key (denoted by e EncðmÞ). In particular, EncðÀmÞ can be computed in this manner.
This scheme should have semantic security; that is, a cipher text leaks no information about the original message (Goldwasser and Micali, 1984). For example, we can use either the Paillier cryptosystem (Paillier, 1999) or the 'lifted' version of the ElGamal cryptosystem (ElGamal, 1985); now the second operation can be realized by repeating the first operation È. Figure 3 illustrates an outline of performing an additive operation on a user's value m 1 and a server's value m 2 by the additively homomorphic encryption. In the first step, the user generates two keys: a secret key and a public key, and the user sends the public key to the server. In the second step, the user encrypts m 1 by the public key and sends a ciphertext Encðm 1 Þ to the server. In the third step, the server Fig. 1. Information flow of a search on a recursive search data structure such as (P)BWT. For ith iteration, the user sends q½i; f iÀ1 , and g iÀ1 , the server returns v q½i ½f iÀ1 and v q½i ½g iÀ1 , and the user updates The server sends a ciphertext c to the user. In the fourth step, the user obtains m 1 þ m 2 by decrypting c.
It goes beyond the scope of this paper to review the details of these cryptographic techniques and the reader is referred to a book (Yi et al., 2014) on homomorphic encryption. A typical addition operation in the ElGamal cryptosystem takes about 2 Â 10 À7 s on a single CPU based on AIST's ElGamal encryption library (https://github. com/aistcrypt/Lifted-ElGamal).

Recursive oblivious transfer by random rotations
To protect the privacy of the database, we propose a technique for recursively querying a data structure without obtaining information about intermediate results. Let us define the recursive oblivious transfer problem as follows: MODEL 1. A user has a private value 0 x 1 < N and a server has a private vector v of length N. Let us denote x kþ1 ¼ v½x k and the user is allowed to access the server ' À 1 times. After the calculation, the user learns only x ' and the server learns nothing about x 1 ; . . . ; x ' .
Here we explain our idea by extending a simple linear communication size OT where the user aims to know the tth element of the server's vector v. Figure 4 illustrates the oblivious transfer algorithm based on additive homomorphic encryption. In the initialization step, the user generates a public key and a secret key and sends the public key to the server. The user creates a bit vector: and sends the following encrypted vector to the server.
The server computes and sends c to the user.
The user computes DecðcÞ and obtains v½t by using the secret key because Now we consider the case that the server does not leak v½t, but allows the user to access v½v½t. Our idea is that the server generates a random value r 2 f0; 1; . . . ; N À 1g and returns the cipher text where ða þ bÞ mod N denotes addition in a number field modulo N. The user decryptsĉ to know a randomized result ðv½t þ rÞ mod N , and performs the next query: Note thatq is the r-rotated permutation of the 'true' query: Therefore, denote Permðq; rÞ as the permutation of a vector q such that ith element moves to ðði À rÞ modN ÞÞth position, the server  . . . ; x5g illustrated in (a) is sorted by the algorithm described in Durbin (2014) to obtain the positional prefix arrays A illustrated in (b). Each element P i;j of PBWT matrix illustrated in (c) is ðj þ 1Þth letter of sequence A i;j . By computing rank operations with regard to kth query letter on P Á;kÀ1 , one can update an interval corresponding to k-mer match between the query and X. In this figure, the search starts from fourth allele. The first interval ðf 1 ; g 1 is initialized by rank operations on P Á;3 with regard to first query letter '1'. The second interval ðf 2 ; g 2 is obtained by rank operations on PÁ;4 with regard to the second query letter '0' and ðf1; g1. Similarly, the third interval ðf3; g3 is obtained by rank operations on PÁ;5 with regard to the third query letter '0' and ðf2; g2. See Sections 2.2 and 3.3 for more details can correctly recover 'true' query q 0 in its encrypted form by the following permutation:Ẽncðq 0 Þ ¼ PermðẼncðqÞ; rÞ. In this way, the server correctly computes an encrypted v½tth element by without learning any information about the user's query. By recursively applying these calculations, the user can obtain x kþ1 according to Model 1. The complete algorithm implementing this idea is given in Algorithm 1. It uses a function ROT for rotating the server's results to conceal intermediate query results in order to protect the database.

PBWT-sec: Privacy-preserving search on genotype databases
In this section, we introduce a practical genotype database search based on recursive oblivious transfer and PBWT. We only introduce the algorithm to search for the longest match starting from tth column, however, variations are possible and would allow for a variety of different search types (see also Durbin, 2014).
To formulate the problem, let us consider a set X of M haplotype sequences x i , i ¼ 1; . . . ; M over N genomic positions indexed by k ¼ 1; . . . ; N, and a query q which is a user's haplotype sequence over the same N genomic positions. We denote kth allele of a sequence x i by x i ½k. Given two indices k 1 and k 2 , we say that there is a match between q and x i from k 1 to k 2 , if q½k 1 . . . q½k 2 À 1 ¼ x i ½k 1 . . . x i ½k 2 À 1. We say that the match is set-longest at k 1 if there is no match between q and any sequence x j (possibly with j ¼ i) from k 1 to k 2 þ 1.
The goal is to find a set-longest match at a given position t between q and X in a privacy-preserving manner. Here, we consider the case that the user's private information is the query string and the position t is not the user's private information. We later introduce the case that the both the query string and t are user's private information. The formal description of the model is described as follows: MODEL 2. The user is a private haplotype sequence holder, and the server is a holder of a set of private haplotype sequences. The user learns nothing but a set-longest match at a given position t between the query and the database while the server learns nothing about the user's query. t is not a user's private information and the server knows it.
Let us remember how to search the set-longest match in nonprivacy-preserving manner. PBWT involves a matrix P 2 N MÂN that stores well-compressible information in an efficiently searchable form. It is created from the genotype matrix X by algorithms described in Durbin (2014) such that ith column is ði þ 1Þth letters of sequences sorted by i reverse prefix (i.e. sorted from ith letter to first letter). In order to compute the match starting from the first allele, P has 0th column P Á;0 ¼ ðx 1 ½1; . . . ; x M ½1Þ T . By using rank dictionary operations on P (see below), one can search a match between a query and X. When operating on P one computes updates of intervals using the following two quantities (see Durbin, 2014 for more details): (i) The rank dictionary for sequence S for letter c 2 R at position t: where R is the alphabet of S. (ii) The table CF counting occurrences of letters that are lexicographically smaller than c in S by CF c ðSÞ ¼ X r<c Rank r ðS; NÞ : Based on these two quantities, we can compute the updates ðf kþ1 ; g kþ1 using two simple operations where we denoted the kth column vector by P Á;k . Let us define a look-up vector v c for the column k where for c 2 R. Then, updating an interval is equivalent to two look-ups in the vector v c : Given a position t and a PBWT P of the database sequences, the first match is obtained as an interval ðf 1 ¼ v c ½0; g 1 ¼ v c ½M where c ¼ q½1 and v c is a look-up vector for ðt À 1Þth column of P (see the definition of v c in Equation (1)). The match is extended by one letter by an update of the interval. The update from the kth interval to ðk þ 1Þth interval is conducted by specifying c ¼ q½k þ 1, re-computing v c for ðk þ 1Þth column of P and referring v c ½f k and v c ½g k as f kþ1 and g kþ1 (see Equation (2)). The set-longest-match is found Algorithm 1. Recursive oblivious transfer 1: function PrepQuery(t, N) 2: Server computes:ĉ ROT(ẼncðqÞ; v, r, r 0 , N) 28: Server sets: r 0 r 29: Server sendsĉ to user 30: User computes: t DecðĉÞ 31: end while 32: User obtains x ' ¼ t.
when f ¼ g. An outline of the search using PBWT is illustrated in Figure 5.
In order to achieve the security described in the model 2, for each update, the user has to specify c without leaking c to the server, and obtain only v c ½f and v c ½g without leaking f and g. To satisfy the second requirement, the user accesses the server's v c through the function ROT, which allows the user to obtain a specific element in the specified vector. To achieve the first requirement, the server computes all possible intervals (i.e. computing (f, g] for the all case of c ¼ 0; . . . ; jRj À 1). This allows the user to obtain the correct interval, however, the sever leaks extra information (i.e. intervals for c 6 ¼ q½k). To avoid this, the user sends Encðq½kÞ, and the server adds a conditional randomization factor r Â ðq½k À cÞ to f and g with different random value r for all c 2 R. Note that this factor becomes equivalent to 0 iff. q½k ¼ c, and user only obtains the interval for c ¼ q½k.
In order to identify the set-longest match, the user has to know if f ¼ g. The user cannot compute the identity of f and g directly from the server's return, because ROT returns a value which is a random value to the user (but the 'true' value is recovered in encrypted form only at the server side). Therefore, the server also sends an encrypted flag d which shows whether or not f ¼ g. Since f and g are represented as indices of q 0 f ¼ Permðq f ; r 0ðf Þ Þ and q 0 g ¼ Permðq g ; r 0ðgÞ Þ (see the functions PrepQuery and ROT), the server computes d by following: where r i is a random value. DecðdÞ is equal to 0 iff. q f ¼ q g . See Supplementary Algorithm S5 which defines a function isLongest. In addition to finding a set-longest match at t, it is convenient to find a longest substring which matches to at least e sequences. This operation enables to avoid detecting unique haplotype and provides eanonymity result and is implemented by replacing the function: isLongest by another function: isLongestGT e which computes flags each of which shows if the interval matches to 0; . . . ; e À 1 respectively and returns shuffled flags, and the user knows the result by checking if there is a flag which is equal to zero. See Supplementary Algorithm S5 for more details. The detailed algorithm of PBWT-sec is shown in Algorithm 2.

Concealing the search position
By the algorithm introduced above, the match position t needs to be provided to the server. Let us consider the case that t needs to be concealed (e.g. the used would not like to reveal which gene is analyzed). In practical genotype database search, it is often sufficient for the user to hide t in a set of multiple columns. Therefore, here we assume the following security model. MODEL 3. The user is a private haplotype sequence holder, and the server is a holder of a set of private haplotype sequences. The user has a vector of D positions T ¼ ðt 1 ; . . . ; t D Þ. The user learns nothing but a set-longest match at a given position t 2 ft 1 ; . . . ; t D g between the query and the database while the server learns nothing about the user's query string. The server knows T but cannot identify which element the user queries.
Conceptually, the user could query multiple positions at the same time to conceal the search position. In the extreme case the user would query all search positions to avoid leaking any information about t. However, every answered query would leak more information from the database and querying would become computationally prohibitive. We therefore propose joint processing using OT that simultaneously uses multiple search positions. Let us define V c as another look-up vector for a letter c as follows: Public input: Problem size M & N; alphabet R ¼ f0; 1; ::; jRj À 1g, start position t 2 f1; . . . ; Ng Private input of user: A query sequence S of length ' Private input of server: PBWT matrix P 2 N MÂN 0. (Key setup of cryptosystem) User generates key pair ðpk; skÞ by key generation algorithm KeyGen for additive-homomorphic cryptosystem and sends public key pk to server. 1. (User initialization) Set initial interval ðf ¼ 0; g ¼ M.

(Recursive search)
Initializes query and position index: a. (Query entry) The user performs the following steps: Prepares next query: search. Note that (V c ½o j ; . . . ; V c ½o j þ M) corresponds to v c for t j th column. The algorithm for the Model 3 is designed by replacing the lookup tables v c by V c (see Step 2a, item 1 in Algorithm 2) and initializing f and g by o x and o x þ M, respectively, where t ¼ t x (see Step 1 in Algorithm 2). As a result the tables get D times larger which has an impact on computing requirements and data transfer size (see Section 3.7). We therefore suggest using this algorithm for small D.

Reducing communication size
As we will describe in the Complexity analysis in the following section, the PBWT-sec algorithm using standard OT requires OðMjRjÞ in communication size per iteration in the best case, which makes the core algorithm less practical. We propose to use an algorithm for sublinear-communication OT (SC-OT) proposed in Zhang et al. (2013). Using this approach we can reduce the communication size of PBWTsec to Oð ffiffiffiffiffiffiffiffiffiffi ffi MjRj p Þ (best case). Here, we only outline the key ideas of SC-OT and its adaptation of PBWT-sec. In the SC-OT, the one encodes the position t by a two dimensional representation: where Á d e denotes the ceil of the argument. The user sends Encðt 0 Þ andẼncðqÞ to the server, wherẽ The server obtains random valuesr k ; k ¼ 0; . . . ; ffiffiffiffi ffi N p AE Ç À 1, and computes and sends c ¼ ðc 0 ; . . . ; c ffiffiffi N p d eÀ1 Þ to the user. The user knows the result by the decryption: Decðc t0 Þ. Note that Encðt 0 À kÞ ¼ Encð0Þ iff. t 0 ¼ k, therefore the decryption of c i becomes a random value when i 6 ¼ t 0 .
In order to apply bit-rotation technique naturally to SC-OT, the server needs to return v½t in the same two dimensional representation. The key idea here is that the server creates v 0 and v 1 where v 0 ½i ¼ v½i= ffiffiffiffi ffi N p AE Ç and v 1 ½i ¼ ðv½iÞ mod ffiffiffi N p d e ; i ¼ 0; . . . ; N À 1, and searches on both v 0 and v 1 . Similar to the linear communication size function ROT, the removable random factors are added to server's returns. More details on SC-OT is given in Section S1. The complete algorithm for privacy-preserving search based on SC-OT is given in Supplementary and Algorithm S2.

An exhaustive baseline algorithm
There are a few related works in regard to finding a DNA substring match (Baldi et al., 2011;Cristofaro et al., 2013), however, the goal of PBWT-sec is to find the set-longest prefix match from a set of aligned sequences while those works aim to find a fixed-length approximate substring match between two sequences. Therefore, we will compare our algorithm with a baseline algorithm which can find the set-longest prefix match on the basis of exhaustive enumeration of k-mers. This baseline algorithm serves as a proxy for the other conceptually similar algorithms.
In order to identify the match, the user queries the server about the presence of a k-mer. Here, the server stores all k-mers, there are OðjRj k Þ of them, and we use SC-OT. Such a strategy is efficient for short queries as jRj k is not too large. However, the resource requirements will be dominated by queries for large k and the algorithm quickly gets intractable.

Complexity
Most of the computing and transfer on server and user side is related to the encryption/decryption and the computational cost of the search is negligible. While PBWT requires essentially Oð1Þ to update the intervals per iteration, PBWT-sec needs to conceal the query and requires MjRj operations on the server, where M is the number of sequences in the database and jRj is the size of the alphabet. When multiple queries are performed at the same time, i.e. D > 1, the effort increases linearly in D, i.e. the server sides compute effort is Oð MDjRjÞ per iteration. When using SC-OT, the communication size and effort for the user is Oð ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi MDjRj p Þ (see Section 3.5 and Supplementary Section S1 for details). Table 1 summarizes the time, data transfer overhead and space complexities of the PBWT-sec, when the server's PBWT is M Â N matrix consisting of a set of alphabet letters R and the user's query length is ' and the number of queries positions is D (including D -1 decoy positions; see Section 3.4 for details). For the purpose of comparison, we consider the method outlined in Section 3.6 that achieves the same security and utility as PBWT-sec. Since the complexity of the exhaustive approach is exponential to the query length, its performance deteriorates quickly for long matches. On the other hand, the time and data transfer overhead complexity of the PBWT-sec are linear and sub-linear to the query length, which enables the user to find a long match efficiently.

Security notion
In this paper, we assume the security model called Semi-honest model where both parties follow the protocol, but an adversarial one attempts to infer additional information about the other party's secret input from the legally obtained information. The semantic security of the encryption scheme used in the protocol (see Section 3.1) implies immediately that the server cannot infer any information about the user's query q during the protocol. Also, the user cannot infer any information about server's return except for the result.
Another security model is called Malicious model where an adversarial party cheats even in the protocol (e.g. by inputting maliciously chosen invalid values) in order to illegally obtain additional information about the secret. Here we briefly describe one example of an illegal access based on the Malicious model. In our protocol, the user needs to create a bit vector q of N that includes a bit that is 1 and the rest of the N -1 bits are 0. If the malicious user creates a non-bit vector: where x is a large integer, the server returns c ¼ Encðv½i þ x Á v½jÞ. When x is larger than any element of v, the user can infer v½i by Both algorithms use SC-OT. M is the number of haplotype sequences (server side), D is the number of queried positions (including D -1 decoy position to conceal the query position), ' is the length of query and jRj is the alphabet size.
ðDecðcÞÞ mod x and v½j by DecðcÞ=x. (For example, if x ¼ 100 and DecðcÞ ¼ 821, the user can detect v½i ¼ 21 and v½j ¼ 8.) Thus the server leaks two elements of v by a single query.
In this study, we do not discuss such cases in detail; however, we would like to mention that it is possible to design an algorithm for the Malicious model with a small modification. In order to avoid such attacks, the server needs to verify if the user sends a bit vector which includes only one bit that is 1 and rest of the bits are 0. To achieve this, we suggest using a cryptographic technique called Non-Interactive Zero Knowledge Proofs which enables the user to convince a server that each ciphertext EncðmÞ corresponds to a value m 2 f0; 1g, but does not leak any information about which of 0 and 1 is m. Among several algorithms, Sakai's algorithm (Sakai et al., 2013) has such a property. By using the algorithm, the server knows whether or not q½i 2 f0; 1g. To return a correct result only if q includes only a single 1, it is sufficient for the server to add w ¼ r ð where r is a random value. Note that w ¼ Encð0Þ iff. q i ¼ 1 and q j ¼ 0 for 0 j < N and i 6 ¼ j.

Experiments
In this section, we evaluate the performance of the proposed method on the datasets created from the chromosome 1 data from the 1000 Genomes Project phase 1 data release which consists of 2184 haploid genomes (The 1000 Genomes Project Consortium, 2012). In our experiments and as in Durbin (2014), we used alleles having SNPs, but we did consider indel variants. We used all 2184 genomes of original data for all the experiments.
We implemented the proposed algorithm in C þþ based on an open source C þþ library of elliptic curve ElGamal encryption provided by AIST. Our implementation supports communication over the network. We used the standard parameters called secp192k1 (SECG curve over a 192-bit prime field), according to the recommendation by The Standards for Efficient Cryptography Group. For comparison, we also implemented an exhaustive baseline method (see Section 3.6) that achieves the same security and utility as PBWT-sec. In order to perform a fair comparison, both PBWT-sec and the exhaustive method used the same SC-OT module where computation of c k (see Algorithm 1) is simply parallelized by OpenMP.
In the first experiment, the user selected a true start position together with 49 decoys (see Section 3.4 for details), and both PBWTsec and the exhaustive method were run with the same computational setting: the user used a single thread of a laptop computer equipped with an Intel Core(TM) i7 3.00 GHz CPU and 16 GB memory, and the server used more than eight threads of another laptop equipped with an Intel Core(TM) i7 2.60 GHz CPU (four cores with hyper-threading) and 8 GB memory. Those two computers communicated over the network. Figures 6 and 7 show run time and data transfer overhead of PBWT-sec and of the exhaustive method. The observed run time and data transfer size of PBWT-sec is linear in the query length, while that of the exhaustive approach is exponential. For query lengths larger than 30 bit, the computation of the exhaustive method did not finish within 24 h. These results fit the theoretical complexity described in Section 3.7. We also evaluated performance of the runtime of PBWT-sec when the user selected 0, 4, 9, 14, and 49 additional decoy positions. The search with a typical query of length 25 SNP positions and no decoy required no more than 15.5 s including communication overhead (Table 2).
The user's run time of PBWT-sec is relatively small, making it suitable for a practical case where computation power in a server side is generally stronger than that of user side. Since the memory usage of PBWT-sec does not depend on query length, it uses less PBWT-sec (user to server) PBWT-sec (server to user) Exhaustive (user to server) Exhaustive (server to user)  The server used all the four cores with hyper-threading while the user used a single thread. All included communication overhead. D is the number of positions queried simultaneously to conceal the query position (if required). than 60 MB while that of the exhaustive method exponentially increases according to the query length and required 6 GB when the query length is 25 bit.
Although the exhaustive method is efficient for short queries, we consider that PBWT-sec is more practical when taking into account that the bit length of a unique substring for a human genome is greater than 31 bits. Moreover, since there are large linkage blocks, even queries with more than 100 bits would not always lead to unique matches in the 1000 genomes data. Hence, the exhaustive search strategy would either not always be able to return a unique match or would be very inefficient. The proposed iterative privacypreserving technique is efficient also for long queries.
In the second experiment, we evaluated the performance of the run time of PBWT-sec on a compute node equipped with four CPU sockets (Intel Xeon 2.40 GHz CPU; total of 32 cores with hyperthreading). In this experiment, the user also selected 0, 4, 9, 14 and 49 additional decoy positions. For environmental reasons, we did not perform communication over the network and the data was transferred by file I/O which is also included in run time.
Although the current implementation is a prototype and there is room for improvement in terms of parallelization, the server's run time was at an acceptable level in practical configurations (Table 3). We note, that with improvements in parallelization, the server run time may be reduced to 3-4 s.

Conclusion
In this paper, we have proposed a novel approach for searching genomic sequences in a privacy-preserving manner. Our approach combines an efficient data structure that can be used for recursive search and a novel approach for recursive oblivious transfer. It achieves high utility and has strong security features and requires acceptable compute and communication resources.
The developed novel algorithm can find the longest match between a query and a large set of aligned genomic sequences indexed by PBWT. We implemented our algorithm and tested on the dataset created from the 1000 Genomes Project data (The 1000 Genomes Project Consortium, 2012). Compared to an exhaustive baseline approach, our algorithm, named PBWT-sec, was orders of magnitude more efficient both in run time and data transfer overhead for practical query sizes. When the prototype program was run on laptop machines, the total run time including communication time over the network was 15.5 s for searching on 2184 genomes without concealing the query position. Searches with a concealed query position using a compute node took between 18.6 and 133 s depending on the level of privacy.
As the original data structure supports many useful search options such as wild card search and set maximal search, PBWT-sec could also support those options by using the same techniques used in the original structures in combination with cryptographic techniques, including OT. Moreover, the approach could be easily applied for BWT and has a potential to be applied for other recursively searchable data structures.
To the best of our knowledge, the proposed algorithm is the first that allows set-maximal search of genomic sequences in a privacypreserving manner for user and database. We note that the implementation can still be improved and the overall run time can likely be reduced to not more than a few seconds per query. This would make it practical to use our approach in a genomic Beacon (see GA4GH's Beacon Project) that would allow the privacy-preserving search for combinations of variants. It also appears practical to use our approach to enable search by a user that has access to his/her genomic sequence and would like to query the database, for instance, for information related to disease risk without sharing this information with anybody. Finally, the algorithm can also be used to facilitate sharing of genetic information across institutions and countries in order to identify large enough cohorts with a similar genetic backgrounds. This is in spirit of the mission of the Global Alliance for Genome and Health. Wall time includes server (% 90%) and user time (% 10%). D is the number of positions queried simultaneously to conceal the query position (if required).