Genome assembly using quantum and quantum-inspired annealing

Recent advances in DNA sequencing open prospects to make whole-genome analysis rapid and reliable, which is promising for various applications including personalized medicine. However, existing techniques for de novo genome assembly, which is used for the analysis of genomic rearrangements, chromosome phasing, and reconstructing genomes without a reference, require solving tasks of high computational complexity. Here we demonstrate a method for solving genome assembly tasks with the use of quantum and quantum-inspired optimization techniques. Within this method, we present experimental results on genome assembly using quantum annealers both for simulated data and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document}ϕX 174 bacteriophage. Our results pave a way for an increase in the efficiency of solving bioinformatics problems with the use of quantum computing and, in particular, quantum annealing. We expect that the new generation of quantum annealing devices would outperform existing techniques for de novo genome assembly. To the best of our knowledge, this is the first experimental study of de novo genome assembly problems both for real and synthetic data on quantum annealing devices and quantum-inspired techniques.

www.nature.com/scientificreports/ reconstructing the original genome. This approach is widely used in assemblers (e.g., see Ref. 13 ) and becomes suitable for the single-molecule sequencing technologies 14,16 . It is well known that finding a Hamiltonian cycle belongs to the class of NP-complete computational problems, for which finding an efficient solution is very hard. That is why the other graph representation is applied to the analysis of the sequencing data. This approach is related to the concept of De Bruijn graphs (DBG) 17 . The idea is to construct a graph based on the fragmentation of reads down to smaller sequences called k-mers, where k is the length of subsequence. These k-mers are aligned using (k − 1) sequence overlaps. Each node in the resulting k-mer graph represents a certain k-mer, while edges correspond to the overlaps between the k-mers. Thus, each read is represented by the sequence of the connected vertices (so-called overlapped k-mers). To obtain the original genome sequence one needs to find the Eulerian path, i.e., the path that visits each edge exactly once, but is allowed to revisit any vertex. The comparison of the OLC and DBG approaches are discussed in more detail in Ref. 15 . As was mentioned above, the time of de novo assembling may be crucial for many applications. Possible improvements and speed-ups of the DBG approach, employing the Eulerian path, have been considered, whereas methods based on finding the Hamiltonian path in the OLC approach are less studied.
Quantum computers are a new generation of devices that use quantum phenomena, such as superposition and entanglement, for solving computational tasks. It is believed that quantum computers 18 have a great potential to outperform existing technologies in various problems 19 , such as simulating complex systems 20 , machine learning 21 , and optimization 22 . Ongoing research are related to the question how quantum computers could be used for computational biology and bioinformatics 23 . One of the algorithms that can be realized using quantum computers is Grover's search 24 , which can be used as a subroutine for sub-sequence alignment with a quadratic speed-up 25 . There is increased activity at the interface of machine learning and quantum computing in the computational biology domain 23,26,27 . Being of extreme interest from the viewpoint of obtaining polynomial and exponential computational speedups, the suggested quantum algorithms 23-28 require both a significant number of qubits and quite low error rates, which are beyond the capabilities of existing noisy intermediate-scale quantum (NISQ) devices.
Although many different implementations and models of quantum computing are in development, one of the approaches, which is based on quantum annealing, deserves a special attention. This is due to the fact that the hardware for its implementation is available (it is produced by D-Wave System) [29][30][31] . However, the ability of realistic quantum annealing devices to demonstrate computation speedups is still a subject of ongoing research 31,[33][34][35][36][37] . An interesting outcome of these debates is the appearance of a new generation of quantum-inspired (digital annealing) algorithms, which are essentially classical but appear as a result of analysing the role of quantum phenomena in solving computational tasks [38][39][40] . Results on the comparison between available quantum annealers and quantum-inspired algorithms on realistic optimization problems have been obtained 40 . We note that quantum annealing has been applied to various real-world tasks, including computational biology 41 , exploration of the conformational landscape of peptides and proteins 42 , and genome sequence alignment 28 .
Here we investigate the de novo genome assembly problem within the framework of quantum annealing. On the one hand, the OLC approach is less sensitive to sequencing errors and repeats than the DBG approach which leads to the higher quality of the assembly. On the other hand, the main disadvantage of the OLC-based algorithms is their computational inefficiency 15 . Particularly, as we need to solve the NP-hard problem of finding a Hamiltonian path in the graph. Quantum computations have already shown their potential suitability for solving such problems. The above mentioned arguments motivate us to focus on the OLC formulation of the de novo genome assembly problem. The main step in our study is to map the genome assembly problem in the framework of OLC graphs to a quadratic unconstrained binary optimization (QUBO) problem, which can be then efficiently embedded in the quantum annealing architecture (see Fig. 1). We also show that the genome assembly problem can be solved with the use of quantum-inspired optimization algorithms. We note that our idea is to use quantum optimization for de novo sequencing, while other problems related to the analysis of genetic data are beyond the scope of the present work.

QUBO reformulation of the genome assembly problem
The growing interest in the use of quantum computing devices (and in particular to quantum annealers) is related to their potential in solving combinatorial optimization problems. It is widely discussed that the potential of quantum annealing is rooted in the quantum effects that allow us to efficiently explore the cost-function landscape in ways unavailable to classical methods. Therefore, the important stage is to map the problem of interest to a Hamiltonian, which maps the binary representation of a graph path into a corresponding energy value. The existing physical implementation of quantum annealing is the D-Wave processor, which can be described as Ising spin Hamiltonian. The Ising Hamiltonian can be transformed into a QUBO problem. Thus, we have to find the mapping to a problem that we would like to solve on the D-Wave quantum processor to the QUBO form. However, establishing correspondence between a problem of interest and the QUBO form may require additional overhead. In particular, in our case the transformation of the OLC graph to a QUBO problem requires the use of additional variables (see below).
In addition, the D-Wave quantum processor has its native structure (the chimera structure). That is why after the formulation of the problem of interest in the QUBO form an additional stage of embedding problem in the native structure of the quantum device is required. So additional overhead in the number of physical qubits, which is related to the representation of logical variables by physical qubits of the processor (that takes into account the native structure), is required (see below).
Formulation of the Hamiltonian path problem. Along the lines of Ref. 43 N} is a set of vertices, and E is set of edges consisting of pairs (u, v) with u, v ∈ V . The solution of the Hamiltonian path problem is represented in the form of N × N permutation matrix X = (x v,i ) , whose unit elements x v,i represent the path going through the vth node at the ith step. Then, we assign each element x v,i of the matrix X a separate logical variable (spin) within an optimization problem. Note, that this representation results in polynomial overhead in the number of logical variables of the QUBO problem: The solution for N-vertex graph requires N 2 logical variables. The resulting Hamiltonian of the corresponding QUBO problem takes the following form: where A > 0 is a penalty coefficient. The first two terms in Eq. (1) ensure the fact that each vertex appears only once in the path, and there is a single vertex at each step of the path. The third term provides a penalty for connections within the path that are beyond the allowed ones. With this QUBO formulation, we are able to run the genome assembly task using quantum annealers and quantum-inspired algorithms. We note that the applicability of the method requires the existence of the Hamiltonian path in the corresponding graph, which is not universally the case for arbitrary genetic data given by an OLC-graph.
Formulation of the Hamiltonian path problem for acyclic graphs. In general, Hamiltonian path mapping is suitable both for cyclic and acyclic directed graphs. However, it is often the case that the OLC graph contains no cycles. It is then possible to further simplify transformation and reduce the qubit overhead. Here, we demonstrate more compact mapping that requires only M logical variables, where M = |E| < N 2 is the number of edges. For the acyclic OLC graph G = (V , E) , let us define a set of binary variables {x u,v } (u,v)∈E that indicate whether an edge (u, v) is included in the Hamiltonian path. Then the corresponding Hamiltonian should include the following two components: www.nature.com/scientificreports/ The first and the second terms in Eq.
(2) assure that each vertex is incident (if possible) with a single incoming and outgoing path edges correspondingly. Although this realization is helpful and can be used for solving genome assembly problems on quantum annealers without polynomial qubit overhead (the encoding requires M variables, where M is the number of edges in the corresponding OLC-graph), the asymptotic computational speed-up versus classical algorithm is not exponential (since there exist efficient classical algorithms for solving this problem).
Embedding to the processor native structure. As it is mentioned before, the logical variables in Hamiltonians (1) and (2) are not necessarily equivalent to physical qubits of a quantum processor. This is due to the fact that quantum processors have its native structure, i.e. topology of physical qubits and couplers between them. For example, each D-Wave 2000Q QPU is fabricated with 2048 qubits and 6016 couplers in a C16 Chimera topology [29][30][31] . Within a C16 Chimera graph, physical qubits are logically mapped into a 16 × 16 matrix of unit cells of 8 qubits, where each qubit is connected with at most 6 other qubits. In order to realize Hamiltonians (1) and (2), whose connections structure may differ from the one of the processor, we employ an additional embedding stage, described in 32 . It allows obtaining a desired effective Hamiltonian by assigning several physical qubits of the processor to a single logical variable of the original Hamiltonian (see Fig. 1c,d). The embedding to the native structure introduces considerable overhead in qubit number relative to the fully-connected model, yet allows solving problems with existing quantum annealers.

Results
Here we apply our method for the experimental realization of de novo genome assembly using quantum and quantum-inspired annealers. As a figure of merit we use the time-to-solution (TTS), which is the total time required by the solver to find the optimal solution (ground state) at least once with a probability of 0.99. We first define R 99 as the number of runs required by the solver to find the ground state at least once with a probability of 0.99. Using binomial distribution one can calculate R 99 as follows: where θ is an estimated success probability of each run.
Then we define TTS it in the following way: where t a for D-Wave is 20µs (default value). Quantum-inspired optimization algorithms can be also used for solving QUBO problems. In our experiments, we employ SimCIM quantum-inspired optimization algorithm 38 , which is based on the differential approach to simulating specific quantum processors called Coherent Ising Machine (CIM; see "Methods"). SimCIM runs on conventional hardware and is easily parallelizable on graphical processing units (GPU). This is the time for simulating a single annealing run using our implementation of SimCIM, measured on Intel core i7-6700 Quad-Core, 64GB DDR4, GeForce GTX 1080.
φX 174 bacteriophage genome. We start with the paradigmatic example of the φX 174 bacteriophage genome 3 . In order to realize de novo genome assembly, we construct the adjacency matrix for OLC graphs and use pre-processing for packing this graph into D-Wave processor. We then transform each adjacency graph into the QUBO matrix according to Eq. (2). For the case of using the D-Wave system we embed the QUBO problem in the native structure of the annealing device (which naturally adds overhead in the number of qubits; see "Methods"). In order to embed φ X 174 graph into the D-Wave processor, we have preliminary conducted manual graph partitioning using classical algorithms implemented in the METIS tool 44 . Then the problem can be solved with the use of quantum annealing hardware by D-Wave and quantum-inspired optimization algorithm.
For each instance, a total of 10 3 anneals (runs) were collected from the processor, with each run having an annealing time of 20 µ s. The total number of instances is 1000 (the process of their generation is described in "Methods"). The results are presented in Table 1. Up to our best knowledge, this is the first realistic-size de novo genome assembly employing the use of quantum computing devices and quantum-inspired algorithms. The presented time is required for finding the optimal solution since only one solution has a right interpretation. We note that the time required for the data pre-processing is not included in Table 1. Details of graph size for each part after manual graph partitioning is presented in Table 2 in "Methods".
Statistics presented for 1000 simulated Phi-X 174 bacteriophage OLC graphs. We use D-Wave hybrid computing mode, which employs further graph decompositions with parallel computing on both classical and quantum backends. D-Wave gives access to the following timing information in the information system: T run (run time), Benchmarking quantum-assisted de novo genome assembly using the synthetic dataset. In order to perform a complete analysis of the suggested approach, we realize the quantum-assisted de novo genome assembly for the synthetic dataset. We generate a synthetic dataset, which consists of 60 random reads of length from 5 to 10 (for details, see "Methods"), 10 problems are generated for every sequence length. We then split each read into k-mers of length 3 and compute adjacency matrix for the corresponding OLC graph using Eq. (2). Finally, we transform each adjacency graph into the QUBO matrix according to our algorithm and minimize it using quantum annealing hardware by D-Wave and quantum-inspired optimization algorithms.
Our goal is to check the applicability of existing quantum annealers to the task of genome assembly, evaluate the upper bound on the input problem size (particularly, the length of the original genome), compare the performance of the D-Wave quantum annealer with a software annealing simulator SimCIM. The choice of tools is motivated by their maturity in terms of quantum dimensionality and compatibility with the original formulation in terms of the optimization problem. The similar routine is realized with the use of quantum-inspired annealing. We test the suggested approach with the simulated data first with the D-Wave quantum annealer (see Fig. 2) and compare our results with quantum-inspired optimization algorithm SimCIM.   Figure 2. Comparison of the performance of quantum and quantum-inspired methods for de novo genome assembly based on synthetic data (10 problems were generated for every sequence length): we compare TTS for quantum device D-Wave and quantum-inspired optimization algorithm SimCIM. www.nature.com/scientificreports/ We note that the D-Wave annealer shows an advantage in genome assembly for short-length sequences, while it cannot be applied for sequences of length 6 and more due to the fact that the decoherence time becomes comparable with the annealing time. What we observe is that the performance of the annealing system is dependent on the properties of input data. While the exhaustive investigation of the nature of D-Wave performance with respect to input data goes beyond the scope of our research, we consider the connectivity of the input graph as one of the critical factors. The D-Wave 2000Q processor is based on Chimera Topology with native physical connectivity of 6 basic qubit cells.
During the experiments on the synthetic dataset, we were able to embed the problems with a maximum sequence of up to 10 nucleotides. This length corresponds to the fully connected graph with 8 nodes (K8,8) and this is the maximum possible graph size, which can be embedded to the chimera lattice using clique embedding tool from DWave Ocean SDK. While Ocean SDK allows using other types of embedding (e.g., minor miner) we observed that clique embedding demonstrates more stable results due to the deterministic nature of embedding in comparison to the minor-miner tool, which is intrinsically based on randomized heuristics contributing to larger deviations across experiments. However, no viable solution that could reconstruct the original sequence was found for sequences longer than 6.

Discussions
In our work, we have demonstrated the possibility of solving the simplified bioinformatics problem of reconstructing genome sequences using quantum annealing hardware and quantum-inspired algorithms. We have implemented the experimental quantum-assisted assembly of φ X 174 bacteriophage genome. On the basis of synthetic data, we have shown that the existing quantum hardware allows reconstructing short sequences of up to 7 nucleotides. In order to use quantum optimization for realistic tasks, the ratio of the decoherence time to the annealing time should be considerably improved. We note that while the decoherence time is not a fundamental limitation of the technology, the realization of quantum annealers with sufficient decoherence time remains a challenge. While D-Wave machines use superconducting quantum circuits, setups based on ultracold Rydberg atom arrays [45][46][47] and trapped ions 48,49 can be also used for the efficient implementation of quantum annealing and other quantum optimization algorithms. Specifically, the system of Rydberg atom arrays has been studied in the context of solving the maximum independent set problem 46,47 , which is NP-hard. For longer sequences, as we have demonstrated, it is possible to use quantum-inspired algorithms that are capable of solving more complex problems using classical hardware.
We note that our work is a proof-of-principle demonstration of the possibility to use existing quantum devices for solving the genome assembly problem. The problem scale considered in this paper is still far from real sequences ( ∼130 kilo-base pairs for primitive bacterias) and is lacking numerous complications, such as errors in sequence reads and handling repeating sequences. However, the proposed method demonstrates that newly evolving computing techniques based on quantum computers and quantum-inspired algorithms are quickly developing and can be soon applied in new areas of science.
Limitations of existing quantum hardware do not allow universally outperform existing solutions for de novo genome assembling. At the same time, one of the most interesting practical questions is when one can expect computational advantages from the use of quantum computing in genome assembling tasks.
We note that in real-life conditions a number of additional challenges arise. Examples include errors (random insertions and deletions, repeats, etc.), genome contaminants (pieces of the genome not related to the subject of interest), polymer chain reaction artifacts, and others require additional post-processing steps. These problems are beyond the scope of our proof-of-principle demonstration and they should be considered in the future. Another complication comes from the fact that temperature and other noise effects play a significant role in the case of the use of realistic quantum devices. Thermal excitation and relaxation processes affect performance. Our further directions include optimization of the QUBO model for more compact spin representation and integration of error model into our algorithm. Solving these two issues can enable reconstruction of real sequences using the quantum approach.

Methods
Quantum annealing protocol. The beginning Hamiltonian of the D-Wave processor is a transverse magnetic field of the following form: where σ x i is the Pauli x-matrix, which acts on ith qubit. The problem Hamiltonian can be encoded to the following Ising Hamiltonian: where h i describe local fields, J ij stands for couplings, σ z i are the Pauli z-matrices, and E is the set of edges. One can see that H P is of diagonal form, so σ z i can be treated as spin values {σ z i = ±1} . For a given spin configuration σ z i the total energy of the system is given by H P , so by measuring the energy one can find a solution to the problem of interest.
Quantum annealing can be applied to any optimization problem that can be expressed in the QUBO form. The idea is then to reduce the problem of interest to the QUBO form. where {σ z i = ±1} . For solving the problem on the D-Wave quantum processor, all h i and J ij values are scaled to lie between −1 and 1. As a result, the processor outputs a set of spin values {σ z i = ±1} that attempts to minimize the energy, and the lower energy indicates better solution of the optimization problem. We note that Ref. 43 provides a method for QUBO/Ising formulations of many NP problems.
Quantum-inspired annealing using SimCIM. SimCIM is an example of a quantum-inspired annealing algorithm, which works in an iterative manner. It can be used for sampling low-energy spin configurations in the classical Ising model. The algorithm treats each spin value s i as a continuous variable, which lie in the range [−1, 1] . Each iteration of the SimCIM algorithm starts with calculating the mean field which act on each spin by all other spins ( b i is an element of the bias vector). Then the gradients for the spin values are calculated according to �s i = p t s i + ζ � i + N(0, σ ) , where p t , ζ are the annealing control parameters and N(0, σ ) is the noise of the Gaussian form. Then the spin values are updated according to s i ← φ(s i + �s i ) , where φ(x) is the activation function After multiple updates, the spins will tend to either −1 or +1 and the final discrete spin configuration is obtained by taking the sign of each s i . Bacteriophage simulations. We use Grinder 50 to simulate raw reads from φ X 174 bacteriophage complete genome (NCBI Reference Sequence: NC 001422.1). To simplify the task and make it feasible for quantum computing we generate 50 reads in each run of simulations. In our proof-of-concept research, we are focused on finding the Hamiltonian path in OLC graph using quantum and quantum-inspired annealing.
We generate the raw reads with no sequencing errors and the length of each read is equal to 600 base pairs. We build the OLC graph using the pairwise alignment of the raw reads implemented in minimap2 package 51 . We run minimap2 with the predefined set of parameters ava-ont and k = 10 . We apply miniasm 14 to the same data as the benchmark assembler, which uses OLC graphs.
For experiments with quantum annealing, we use public access to D-Wave 2000Q via Leap SDK. We evaluate the impact of tunable parameters (particularly, annealing time) on the final solution quality; however, no significant improvement was discovered against default values, so annealing time was set to 20 µ s (default value). The number of annealing runs is set to 10 3 (maximum possible value). During our experiments we use mostly the standard configuration of the D-Wave processor, so we do not have any specific requirements on the weights/ couplers in the model.
x for |x| ≤ 1; x/|x| otherwise www.nature.com/scientificreports/ Synthetic dataset graphs, which consist of reads no longer than 7 nucleotides (25 graph nodes), are small enough to fit into quantum annealer, so we can use DW 2000Q 5 backend (pure quantum mode of operation; see the following section). However, the size of the φ X 174 bacteriophage graph (248 vertices) is too large. In order to embed φ X 174 graph into the D-Wave processor, we have preliminary conducted manual graph partitioning using classical algorithms implemented in the METIS tool 44 . It allows splitting the original graph into 3 sub-parts. They are carefully selected so that only a single edge remains between them. The longest path is then calculated separately for each part and concatenated into a single path of the original graph. Each part is still large enough to be computed using the purely quantum mode, so we use the D-Wave hybrid computing mode -hybrid v1 backend. D-Wave hybrid computing mode employs further graph decompositions with parallel computing on both classical and quantum backends. Specifics of such decomposition are not publicly available and physical qubit count is also not shown to the end-user. Details of graph size for each part after manual graph partitioning is presented in Table 2. According to D-Wave Leap specification, hybrid v1 backend automatically combines the power of classical and quantum computation.
Simulations with synthetic dataset. In order to evaluate the performance of the algorithm in a controlled setup, we generated several hundreds of random nucleotide sequences with variable length and performed corresponding transformations as shown in Fig. 3. Further, we eliminated graph duplicates or other trivial cases, where graph structure contained no auxiliary edges. Synthetic dataset graphs up to the length of 7 nucleotides (25 graph nodes) are small enough to fit into quantum annealer, so we can use DW 2000Q 5 backend (pure quantum mode of operation). Finally, we selected 60 sequences that produce unique OLC graphs with comparable complexity.
Note Added. Recently, we became aware of the work reporting studies of the quantum acceleration using gate-based and annealing-based quantum computing 52 .