solving snp problem by Qaoa

QAOA, or Quantum Approximate Optimization Algorithm, is a quantum algorithm designed to solve optimization problems by finding the minimum value of a given objective function. The algorithm involves a time evolution process of a quantum system governed by a Hamiltonian, which is a mathematical representation of the system’s total energy.

                                                                     𝐻 = 𝐻1 + 𝐻2 + 𝐻3 + ⋯ + 𝐻𝑁

The Hamiltonian is decomposed into a sum of simpler Hamiltonians, which can be implemented using quantum gates. The system’s time evolution is then approximated using the Trotter-Suzuki formula, which allows for implementing the Hamiltonian using a sequence of simpler operations.

                                                                     𝑒𝐴 + 𝐵 ≈ (𝑒𝐴/𝑛 𝑒𝐵/𝑛)𝑛

By repeating this process multiple times, the algorithm is able to converge to a solution that approximates the minimum value of the objective function.

          𝑈(𝐻,𝑡,𝑛) = ∏ ∏𝑒−𝑖𝐻𝑘𝑡/𝑛  where 𝑗 = 1 to 𝑛 for k times for the Hamiltonian of the system given by  𝐻 = ∑𝐻𝑘 for k times

Solving Single Nucleotide Polymorphism by QAOA

So basically, in biochemistry, sometimes we have to figure out which sequences in a sample are causing problems and get rid of them. We call this the SNP Assembly Problem. We make a graph where each sequence is a point, and if two sequences are causing problems together, we connect them with a line. We aim to eliminate as few sequences as possible to solve all the problems. This is called finding a minimum vertex cover in the graph, but it’s a tough nut to crack because it’s an NP-complete problem. For example, one common source of genetic differences in humans is called Single Nucleotide Polymorphism (SNP), which is when there’s a change in one base of the DNA. Actually, SNPs make up about 90% of all the genetic differences in human DNA.

The SNP Assembly Problem is about finding the fewest number of SNPs to remove in order to eliminate all conflicts between them. An SNP assembly consists of a set of SNPs, a set of fragments, and a relation between the two indicating if an SNP appears in a fragment or not and, if it does, what value it has. Two SNPs are considered in conflict if there are two fragments where exactly three of their corresponding values match, and the fourth value is the opposite. The goal is to find which SNPs to remove to resolve all conflicts. An example of the problem is shown in a table, and the relation between SNPs and fragments is only known for a subset of them based on experimental data.

Assuming we have obtained the conflict graph:

Now, we can use it as an input to the Quantum Approximate Optimization Algorithm (QAOA) to determine which elements in the graph sequence are duplicated and need to be removed.

				
					import pennylane as qml
from pennylane import qaoa
from pennylane import numpy as np
from matplotlib import pyplot as plt
import networkx as nx



# this is conflict matrix that we are going to have QAOA algorithm to solve it for us
edges = [(0, 4), (4, 3),(3, 1), (3, 5), (2,5)]
graph = nx.Graph(edges)

nx.draw(graph, with_labels=True)
plt.show()


# Conflict Graph As Of Input To Minimum Vertex Cover
cost_h, mixer_h = qaoa.min_vertex_cover(graph, constrained=False)

print("Cost Hamiltonian", cost_h)
print("Mixer Hamiltonian", mixer_h)


def qaoa_layer(gamma, alpha):
    qaoa.cost_layer(gamma, cost_h)
    qaoa.mixer_layer(alpha, mixer_h)


wires = range(6)
depth = 3

def circuit(params, **kwargs):
    for w in wires:
        qml.Hadamard(wires=w)
    qml.layer(qaoa_layer, depth, params[0], params[1])

dev = qml.device("default.qubit", wires=wires)


@qml.qnode(dev)
def cost_function(params):
    circuit(params)
    return qml.expval(cost_h)


optimizer = qml.GradientDescentOptimizer()
steps = 1000
params = np.array([[0.5, 0.5, 0.5], [0.5, 0.5, 0.5]], requires_grad=True)


for i in range(steps):
    params = optimizer.step(cost_function, params)

print("Optimal Parameters")
print(params)



@qml.qnode(dev)
def probability_circuit(gamma, alpha):
    circuit([gamma, alpha])
    return qml.probs(wires=wires)


probs = probability_circuit(params[0], params[1])

plt.style.use("seaborn")
plt.bar(range(2 ** len(wires)), probs)
plt.show()

This code uses the PennyLane library to implement the Quantum Approximate Optimization Algorithm (QAOA) to solve the Minimum Vertex Cover problem for a given conflict matrix represented by a graph.

First, the code defines a graph using the NetworkX library and visualizes it using Matplotlib. Then, it calculates the cost and mixer Hamiltonians for the QAOA algorithm using the qaoa.min_vertex_cover function.

The QAOA circuit is defined using the qaoa_layer function and the circuit function, which takes the QAOA parameters as input and applies a series of Hadamard gates and QAOA layers to the qubits.

The QAOA optimization loop is implemented using PennyLane’s automatic differentiation capabilities and the GradientDescentOptimizer optimizer. It uses the cost_function defined earlier to calculate the cost of the current parameters and updates the parameters with each optimization step.

Finally, the code calculates the probabilities of measuring the qubits in each computational basis state using the optimized QAOA parameters and visualizes the results using Matplotlib.

Result:

The states |19⟩ = |010011⟩ have the highest probabilities of being measured, just as expected. Then I should remove sequences 1, 4, and 5 from the conflict graph to avoid conflicts with the samples. Easy, peasy!

Back