该文主要解析基因组组装的核心算法原理，并探讨量子计算在该领域的创新应用。我们首先拆解De Bruijn图与OLC（Overlap-Layout-Consensus）两类主流组装策略的数学建模逻辑；随后针对OLC算法中重叠图构建引发的组合爆炸问题，提出基于量子计算的优化方案：通过将测序reads的拼接路径建模为旅行商问题（TSP），构建Ising模型的QUBO表达式，并借助相干伊辛机（CIM）实现高效求解。

什么是基因组组装？

基因组组装（Genome assembly）是使用测序方法将待测物种的基因组生成序列片段（即read），然后根据reads之间的重叠区域对片段进行拼接。首先，这些片段被拼接成较长的连续序列（contig），然后将contigs拼接成更长的允许包含空白序列（gap）的scaffolds。通过消除scaffolds的错误和gaps，我们可以将这些scaffolds定位到染色体上，从而得到高质量的全基因组序列。这个过程对于研究物种的基因组结构、功能基因挖掘、调控代谢网络构建以及物种进化分析等方面具有重要意义。

随着测序技术的进步和测序成本的降低，我们能够在短时间内获得大量的测序数据。然而基因组组装是一个复杂的过程，尤其是面临动植物等复杂且规模巨大的基因组的组装需要耗费大量的计算资源，因此，如何高效的分析和处理这些数据成为目前基因组学和生物信息学最亟待解决的问题。

图1 DNA 测序流程图(引用自：Next-Generation Sequencing (NGS)- Definition, Types (microbenotes.com))

基因组组装常见算法

基因组组装算法涉及多种方法，分为有参考基因组的组装和无参考的从头组装（de novo assembly），其中从头组装最常见的两个策略是基于德布鲁因图（De Bruijn graph）和Overlap-Layout-Consensus（OLC）的组装算法。

1. 德布鲁因图（De Bruijn graph）：

德布鲁因图是一种有向图，由节点和边构成。其要求是相邻两个节点的元素错开一个碱基。通过将测序得到的reads以k-mer方式排列，德布鲁因图可以帮助我们拼接成更长的Contig。例如，如果几条reads之间有overlap，那么就容易拼接成一个比较长的Contig。实际情况会更复杂，但我们通常保留主干路径，由主干路径组装成Contig1。

2. Overlap-Layout-Consensus（OLC）：

OLC方法首先构建重叠图，然后将重叠图收束成Contig，最后选择每个Contig中最有可能的核苷酸序列。重叠图表示对于一些序列，前后有若干个碱基是一样的。多个reads片段可以构成重叠图。寻找最短的Contig序列，使得所有的reads都是该序列的子集，通常采用贪心算法实现。

图2 OLC (左)和 DBG 算法示意图（右）(引用自：https://qinqianshan.com/bioinformatics/bin/olc-dbg/)

应用量子计算解决基因组组装问题

1. 问题转化

高质量的基因组装需要较高深度的测序才能满足要求，对于完整的人类基因组组装，通常建议的覆盖度在30x以上。以三代测序（如PacBio和Oxford Nanopore等）用于人类基因组组装为例，通常每个read可以达到几千到几万个碱基，覆盖30x的情况下，可能需要的reads数量在几百万到几千万不等。如此庞大的reads组装在计算重叠图(overlap graph)时通常面临巨大的计算复杂度和算力需求，因此经典计算通常引入重叠过滤、动态规划等方法来降低计算复杂度，如果组装高质量的基因组则需要耗费大量时间，这里我们引入量子算法来加速这一过程。

OLC算法中最关键的一步是通过构建重叠图来确定不同reads间的位置关系，这一步可以转换成组合优化问题的旅行商问题(Traveling Salesman Problem，TSP)。其中不同reads表示不同的城市结点，reads间的重叠长度作为结点间的权重，目标是将所有的reads按照它们之间的重叠长度进行排序，并尝试将它们组装成一个总体重叠长度最长的序列，如图三所示。

图3 测序reads 组装转换为有向图模型

在该模型中我们将所有reads设为节点，节点的边的权重为两个reads按顺序前后拼接的重叠长度，这样就构成一个带权有向图。 目标为寻找总体重叠长度最长的reads序列则转化为寻找权重最大的哈密顿回路。 将带权有向图用邻接矩阵表示：, 元素的下标表示在前、在后的拼接重叠长度。如果重叠长度不大于阈值，或者，则  。决策变量取值为,1表示排在第个位置，否则为0。

我们想得最大重叠长度，目标函数可写成：

 

2. 约束条件处理

一个reads必须而且只能出现在序列中的一个位置，即

序列中的一个位置必须而且只能被一个reads占用，即

将以上两个约束写成QUBO形式，有：

3. QUBO模型构建

如此，由以上的约束条件，我们可以构建一个遍历所有的点的路径的模型（又称哈密尔顿环）：

优化目标是求 H 的最小值。

4. 使用开物SDK进行Ising矩阵生成

加载依赖包

import numpy as np
import kaiwu as kw
from Bio import Entrez, SeqIO
import random
import matplotlib.pyplot as plt
import time
import copy
import re
np.set_printoptions(linewidth=1000)
seed = 0

准备数据，从NCBI上获取一条参考DNA序列

#　从NCBI获取乙肝病毒的X基因并返回序列
#　填写自己的email地址，可以帮助NCBI监控数据获取请求
Entrez.email = user_email
**gene_id ="NC_003977.2"**
# 从GenBank获取Accession=NC_003977.2的基因全序列,此处为乙肝病毒，长度为3182
handle = Entrez.efetch(db="nucleotide", id=gene_id,
                               rettype="gb", retmode="text")

# 使用SeqIO读取GenBank格式的记录
record = SeqIO.read(handle, "genbank")

# 解析获取的记录，并提取DNA序列
DNA_raw = str(record.seq)

# 找到病毒中X基因特征和起止位置，X基因长度为465
for feature in record.features:
if feature.type == "CDS" and "gene" in feature.qualifiers and feature.qualifiers["gene"][0] == gene_name:
            x_start = feature.location.start.position
            x_end = feature.location.end.position

# 释放handle
handle.close()
Gene_ref = DNA_raw[x_start:x_end]

内置一段下载好的序列

Gene_ref = 
'ATGGCTGCTAGGCTGTGCTGCCAACTGGATCCTGCGCGGGACGTCCTTTGTTTACGTCCCGTCGGCGCTGAATCCTGCGGACGACCCTTCTCGGGGTCGCTTGGGACTCTCTCGTCCCCTTCTCCGTCTGCCGTTCCGACCGACCACGGGGCGCACCTCTCTTTACGCGGACTCCCCGTCTGTGCCTTCTCATCTGCCGGACCGTGTGCACTTCGCTTCACCTCTGCACGTCGCATGGAGACCACCGTGAACGCCCACCAAATATTGCCCAAGGTCTTACATAAGAGGACTCTTGGACTCTCAGCAATGTCAACGACCGACCTTGAGGCATACTTCAAAGACTGTTTGTTTAAAGACTGGGAGGAGTTGGGGGAGGAGATTAGGTTAAAGGTCTTTGTACTAGGAGGCTGTAGGCATAAATTGGTCTGCGCACCAGCACCATGCAACTTTTTCACCTCTGCCTAA'

对参考序列模拟测序方法进行切割成`reads_number`份

reads = []  # 序列读取
cutpoint = []  # 切割位置
n = len(Gene_ref)

for i in range(0, reads_number):  # 随机将n个基因序列拷贝切割为两段，得到2n个reads
cut_1 = random.randint(0, len(Gene_ref)-1)
        cut_2 = (cut_1 + round(len(Gene_ref)/2) +
                         random.randint(-40, 40)) % n

# 在序列上任意点两个点，将基因序列分成大致相等的两半，长度±10%
Gene_cut_1 = min(cut_1, cut_2)
Gene_cut_2 = max(cut_1, cut_2)
cutpoint.append(Gene_cut_1)
cutpoint.append(Gene_cut_2)

reads.append(Gene_ref[Gene_cut_1:Gene_cut_2-1])  # 其中一半
reads.append(Gene_ref[Gene_cut_2-1:] + Gene_ref[0:Gene_cut_1])  # 另一半

定义两个序列比对函数，输入两段reads，返回重叠得碱基长度

def align(read1, read2, max_missmatch, cut_off):
        '''
Overlap模块，Reads1为头Reads2为尾进行比对，返回重叠长度。
      read1是前基因片段，read2是后基因片段，max_missmatch为容错位数,cut_off是重叠长度最低阈值。
        read1 = 'AAATTTCCTTC'
        read2 = 'TTCCAGGT'
        mm = 0
        align(read1,read2,mm,1)
        return 3
        align(read1,read2,mm,3)
        return 0
        '''

        l1 = len(read1)
        l2 = len(read2)

        for shift in range(0, l1-cut_off):
            miss_match = 0  # 不匹配位数
            r2i = 0
            for r1i in range(shift, l1):
                if read1[r1i] != read2[r2i]:
                    miss_match += 1  # 不对齐位数加一
                if r2i < len(read2)-1:
                    r2i += 1
                else:
                    break
                if miss_match > max_missmatch:  # 如果不匹配位数超过容错位，跳出
                    break

            if miss_match <= max_missmatch:  # 如果不匹配位数在容错范围内，返回重叠序列长度
                return min(l1-shift, l2)  # 返回重叠序列长度  min(l1-shift, l2)
        return 0

def aligncycle(self, read1, read2, max_missmatch, cut_off):
        l1 = len(read1)
        for shift in range(0, l1-cut_off):
            miss_match = 0  # 不匹配位数
            r2i = 0
            for r1i in range(shift, l1):
                if read1[r1i] != read2[r2i]:
                    miss_match += 1  # 不对齐位数加一
                if r2i < len(read2)-1:
                    r2i += 1
                else:
                    break

                if miss_match > max_missmatch:  # 如果不匹配位数超过容错位，跳出
                    break

            if miss_match <= max_missmatch:  # 如果不匹配位数在容错范围内，返回重叠序列长度
                return l1-shift  # 返回至初始位置

        return 0

把DNA测序reads转化为TSP邻接矩阵，边的权重为overlap，边的方向为行角标→列角标

#构建邻接矩阵
tspAdjMatrix = np.zeros((n, n))
        for r1 in range(0, n):
            for r2 in range(0, n):
                if r1 != r2:
                    tspAdjMatrix[r1][r2] = align(
                        reads[r1], reads[r2], max_mismatch, cut_off)

创建变量矩阵

# 创建qubo变量矩阵。QUBO变量矩阵Reads_Sort(i,j)意义为在第i个基因read排在第j位
#创建变量
Reads_Sort = kw.qubo.ndarray((n, n), "Reads_Sort", kw.qubo.binary)
w = tspAdjMatrix

添加约束项，构建目标函数

#read约束，每个read只属于一个位置
sequence_cons = kw.qubo.quicksum(
            (1-kw.qubo.quicksum(Reads_Sort[v, j] for j in range(n)))**2 for v in range(n))

#位置约束，每个位置只能有一个read
node_cons = kw.qubo.quicksum(
            (1-kw.qubo.quicksum(Reads_Sort[v, j] for v in range(n)))**2 for j in range(n))

 # 哈密尔顿环约束
ham_cycle = sequence_cons + node_cons

# 总COST为叠加位数；DNA片段选取总长度最小的，因此要COST最大，在目标函数里取负号
path_cost = - kw.qubo.quicksum(w[u, v] * (kw.qubo.quicksum(Reads_Sort[u, j] * Reads_Sort[v, j+1]
                                                              for j in range(n-1)) + Reads_Sort[u, n-1] * Reads_Sort[v, 0])
                                       for u, v in edges)

obj = constraint * ham_cycle + path_cost

模型编译，生成Ising矩阵

# 解析QUBO
obj = kw.qubo.make(obj)
# 转化为Ising模型
obj_ising = kw.qubo.cim_ising_model(obj)
# 提取Ising矩阵
matrix = obj_ising.get_ising()["ising"]

5. 使用光量子计算模拟器进行计算

worker = kw.cim.SimulatedCIMOptimizer(
    pump=1.3,
    noise=0.2,
    laps=5000,
    delta_time=0.05,
    normalization=0.3,
    iterations=50
)

output = worker.solve(matrix)
# Sort the results
opt = kw.sampler.optimal_sampler(matrix, output, bias=0, negtail_ff=False)
# Select the best solution
cim_best = opt[0][0]
# If the linear term variable is -1, perform a flip
cim_best = cim_best * cim_best[-1]
print(cim_best)

通过模拟器求解，得到Ising解向量

[ 1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1  1]

6. 还原Ising解向量的spin结果为QUBO变量值

以上我们得到了spin结果，需要将其转化为原始变量的值，才能进一步给出求得的路径。我们使用get_val模块进行转化。

# Get the list of variable names
vars = obj_ising.get_variables()

# Substitute the spin vector and obtain the result dictionary
sol_dict = kw.qubo.get_sol_dict(cim_best, vars)

# Check the hard constraints for validity and path length
seq_val = kw.qubo.get_val(sequence_cons, sol_dict)
node_val = kw.qubo.get_val(node_cons, sol_dict)
ham_val = kw.qubo.get_val(ham_cycle, sol_dict)
print('position cons: {}'.format(seq_val))
print('node_cons cons: {}'.format(node_val))
print('ham_cycle: {}'.format(ham_val))

# Calculate the path length using path_cost
path_val = kw.qubo.get_val(path_cost, sol_dict)
print('path_cost: {}'.format(path_val))

以上我们获得了路径的合法性和长度

position cons: 0.0
node_cons cons: 0.0
ham_cycle: 0.0
path_cost: 50.0

下面我们恢复变量的值并获取最后的DNA序列。

# 获得x的数值矩阵
x_val = kw.qubo.get_array_val(Reads_Sort, sol_dict)

# 找到其中的非0项的脚标
nonzero_index = np.array(np.nonzero(x_val)).T

# 对非零项的顺序脚标进行排序
nonzero_index = nonzero_index[nonzero_index[:, 1].argsort()]

# 得到满足限制条件解,将参数记录下来
validresults = validresults + 1
print(f'DNA sequence {validresults}')
DNA sequence 4 3 1 2

我们定义一个通过解获得得节点顺序转换为DNA序列的函数

def merge_seq(reads, DNAsort):

# consensus模块，按照解顺序组装DNA。给一个序列，返回按此序列拼接的reads
merge = reads[DNAsort[0]]  # 重建基因序列
        n = len(reads)
        sum_overlap = 0  # 拼接的总重叠长度
        r1i = 0
        r2i = 0

        # 按照解的顺序连接基因片段，先从第0个连接到第n-1个
        while r2i < (n-1):
            r2i += 1
            overlap = self.align(
                reads[DNAsort[r1i]], reads[DNAsort[r2i]], 0, 0)
            overlap1 = align(reads[DNAsort[r2i-1]], reads[DNAsort[r2i]], 0, 0)
            sum_overlap = overlap1 + sum_overlap

            # 如果重叠长度大于R2长度，则R1完全覆盖R2，直接跳过，不需要融合。R1不动，R2继续下一个
            if len(self.reads[DNAsort[r2i]]) > overlap:
                merge = merge + reads[DNAsort[r2i]][overlap:]
                r1i = r2i

        # 然后连接第此时最长的尾巴个和第0个
        overlap = aligncycle(
            reads[DNAsort[n-1-r2i+r1i]], reads[DNAsort[0]], 0, 0)
        if overlap != 0:
            merge = merge[:-overlap]
        if len(reads[DNAsort[0]]) > overlap:
            sum_overlap = overlap + sum_overlap
        else:
            sum_overlap = len(reads[DNAsort[0]]) + sum_overlap
        return [merge, sum_overlap]

DNAseq = merge_seq(nonzero_index[:, 0])

量子计算与经典方法得效果对比

我们比较了在不同问题规模下（reads数分别为7、8、9、10）CIM真机和传统求解器Simulated Annealing、Tabu Search的求解时间。对比发现，在同等规模问题上，CIM真机的求解时间要远少于Simulated Annealing、Tabu Search（左图）。同时，随着问题规模增大，CIM真机求解时间相较于Simulated Annealing、Tabu Search更加稳定（右）。

图4 不同问题规模下(reads数分别为7、8、9、10)“天工量子大脑”CIM真机和传统求解器Simulated Annealing、Tabu Search的求解时间对比

以reads 数为7的问题求解为例，我们可以发现“天工量子大脑”CIM真机相较于SA算法，其哈密顿量能快速下降以求得最优解。

图5 Reads数为7时的CIM(红线)、SA(蓝线)哈密顿量演化图

 

 

 

 

 

 

 

 

 

 

 

 

该文主要解析基因组组装的核心算法原理，并探讨量子计算在该领域的创新应用。我们首先拆解De Bruijn图与OLC（Overlap-Layout-Consensus）两类主流组装策略的数学建模逻辑；随后针对OLC算法中重叠图构建引发的组合爆炸问题，提出基于量子计算的优化方案：通过将测序reads的拼接路径建模为旅行商问题（TSP），构建Ising模型的QUBO表达式，并借助相干伊辛机（CIM）实现高效求解。

什么是基因组组装？

基因组组装（Genome assembly）是使用测序方法将待测物种的基因组生成序列片段（即read），然后根据reads之间的重叠区域对片段进行拼接。首先，这些片段被拼接成较长的连续序列（contig），然后将contigs拼接成更长的允许包含空白序列（gap）的scaffolds。通过消除scaffolds的错误和gaps，我们可以将这些scaffolds定位到染色体上，从而得到高质量的全基因组序列。这个过程对于研究物种的基因组结构、功能基因挖掘、调控代谢网络构建以及物种进化分析等方面具有重要意义。

随着测序技术的进步和测序成本的降低，我们能够在短时间内获得大量的测序数据。然而基因组组装是一个复杂的过程，尤其是面临动植物等复杂且规模巨大的基因组的组装需要耗费大量的计算资源，因此，如何高效的分析和处理这些数据成为目前基因组学和生物信息学最亟待解决的问题。

图1 DNA 测序流程图(引用自：Next-Generation Sequencing (NGS)- Definition, Types (microbenotes.com))

基因组组装常见算法

基因组组装算法涉及多种方法，分为有参考基因组的组装和无参考的从头组装（de novo assembly），其中从头组装最常见的两个策略是基于德布鲁因图（De Bruijn graph）和Overlap-Layout-Consensus（OLC）的组装算法。

1. 德布鲁因图（De Bruijn graph）：

德布鲁因图是一种有向图，由节点和边构成。其要求是相邻两个节点的元素错开一个碱基。通过将测序得到的reads以k-mer方式排列，德布鲁因图可以帮助我们拼接成更长的Contig。例如，如果几条reads之间有overlap，那么就容易拼接成一个比较长的Contig。实际情况会更复杂，但我们通常保留主干路径，由主干路径组装成Contig1。

2. Overlap-Layout-Consensus（OLC）：

OLC方法首先构建重叠图，然后将重叠图收束成Contig，最后选择每个Contig中最有可能的核苷酸序列。重叠图表示对于一些序列，前后有若干个碱基是一样的。多个reads片段可以构成重叠图。寻找最短的Contig序列，使得所有的reads都是该序列的子集，通常采用贪心算法实现。

图2 OLC (左)和 DBG 算法示意图（右）(引用自：https://qinqianshan.com/bioinformatics/bin/olc-dbg/)

应用量子计算解决基因组组装问题

1. 问题转化

高质量的基因组装需要较高深度的测序才能满足要求，对于完整的人类基因组组装，通常建议的覆盖度在30x以上。以三代测序（如PacBio和Oxford Nanopore等）用于人类基因组组装为例，通常每个read可以达到几千到几万个碱基，覆盖30x的情况下，可能需要的reads数量在几百万到几千万不等。如此庞大的reads组装在计算重叠图(overlap graph)时通常面临巨大的计算复杂度和算力需求，因此经典计算通常引入重叠过滤、动态规划等方法来降低计算复杂度，如果组装高质量的基因组则需要耗费大量时间，这里我们引入量子算法来加速这一过程。

OLC算法中最关键的一步是通过构建重叠图来确定不同reads间的位置关系，这一步可以转换成组合优化问题的旅行商问题(Traveling Salesman Problem，TSP)。其中不同reads表示不同的城市结点，reads间的重叠长度作为结点间的权重，目标是将所有的reads按照它们之间的重叠长度进行排序，并尝试将它们组装成一个总体重叠长度最长的序列，如图三所示。

图3 测序reads 组装转换为有向图模型

在该模型中我们将所有reads设为节点，节点的边的权重为两个reads按顺序前后拼接的重叠长度，这样就构成一个带权有向图。 目标为寻找总体重叠长度最长的reads序列则转化为寻找权重最大的哈密顿回路。 将带权有向图用邻接矩阵表示：, 元素的下标表示在前、在后的拼接重叠长度。如果重叠长度不大于阈值，或者，则  。决策变量取值为,1表示排在第个位置，否则为0。

我们想得最大重叠长度，目标函数可写成：

 

2. 约束条件处理

一个reads必须而且只能出现在序列中的一个位置，即

序列中的一个位置必须而且只能被一个reads占用，即

将以上两个约束写成QUBO形式，有：

3. QUBO模型构建

如此，由以上的约束条件，我们可以构建一个遍历所有的点的路径的模型（又称哈密尔顿环）：

优化目标是求 H 的最小值。

4. 使用开物SDK进行Ising矩阵生成

加载依赖包

import numpy as np
import kaiwu as kw
from Bio import Entrez, SeqIO
import random
import matplotlib.pyplot as plt
import time
import copy
import re
np.set_printoptions(linewidth=1000)
seed = 0

准备数据，从NCBI上获取一条参考DNA序列

#　从NCBI获取乙肝病毒的X基因并返回序列
#　填写自己的email地址，可以帮助NCBI监控数据获取请求
Entrez.email = user_email
**gene_id ="NC_003977.2"**
# 从GenBank获取Accession=NC_003977.2的基因全序列,此处为乙肝病毒，长度为3182
handle = Entrez.efetch(db="nucleotide", id=gene_id,
                               rettype="gb", retmode="text")

# 使用SeqIO读取GenBank格式的记录
record = SeqIO.read(handle, "genbank")

# 解析获取的记录，并提取DNA序列
DNA_raw = str(record.seq)

# 找到病毒中X基因特征和起止位置，X基因长度为465
for feature in record.features:
if feature.type == "CDS" and "gene" in feature.qualifiers and feature.qualifiers["gene"][0] == gene_name:
            x_start = feature.location.start.position
            x_end = feature.location.end.position

# 释放handle
handle.close()
Gene_ref = DNA_raw[x_start:x_end]

内置一段下载好的序列

Gene_ref = 
'ATGGCTGCTAGGCTGTGCTGCCAACTGGATCCTGCGCGGGACGTCCTTTGTTTACGTCCCGTCGGCGCTGAATCCTGCGGACGACCCTTCTCGGGGTCGCTTGGGACTCTCTCGTCCCCTTCTCCGTCTGCCGTTCCGACCGACCACGGGGCGCACCTCTCTTTACGCGGACTCCCCGTCTGTGCCTTCTCATCTGCCGGACCGTGTGCACTTCGCTTCACCTCTGCACGTCGCATGGAGACCACCGTGAACGCCCACCAAATATTGCCCAAGGTCTTACATAAGAGGACTCTTGGACTCTCAGCAATGTCAACGACCGACCTTGAGGCATACTTCAAAGACTGTTTGTTTAAAGACTGGGAGGAGTTGGGGGAGGAGATTAGGTTAAAGGTCTTTGTACTAGGAGGCTGTAGGCATAAATTGGTCTGCGCACCAGCACCATGCAACTTTTTCACCTCTGCCTAA'

对参考序列模拟测序方法进行切割成`reads_number`份

reads = []  # 序列读取
cutpoint = []  # 切割位置
n = len(Gene_ref)

for i in range(0, reads_number):  # 随机将n个基因序列拷贝切割为两段，得到2n个reads
cut_1 = random.randint(0, len(Gene_ref)-1)
        cut_2 = (cut_1 + round(len(Gene_ref)/2) +
                         random.randint(-40, 40)) % n

# 在序列上任意点两个点，将基因序列分成大致相等的两半，长度±10%
Gene_cut_1 = min(cut_1, cut_2)
Gene_cut_2 = max(cut_1, cut_2)
cutpoint.append(Gene_cut_1)
cutpoint.append(Gene_cut_2)

reads.append(Gene_ref[Gene_cut_1:Gene_cut_2-1])  # 其中一半
reads.append(Gene_ref[Gene_cut_2-1:] + Gene_ref[0:Gene_cut_1])  # 另一半

定义两个序列比对函数，输入两段reads，返回重叠得碱基长度

def align(read1, read2, max_missmatch, cut_off):
        '''
Overlap模块，Reads1为头Reads2为尾进行比对，返回重叠长度。
      read1是前基因片段，read2是后基因片段，max_missmatch为容错位数,cut_off是重叠长度最低阈值。
        read1 = 'AAATTTCCTTC'
        read2 = 'TTCCAGGT'
        mm = 0
        align(read1,read2,mm,1)
        return 3
        align(read1,read2,mm,3)
        return 0
        '''

        l1 = len(read1)
        l2 = len(read2)

        for shift in range(0, l1-cut_off):
            miss_match = 0  # 不匹配位数
            r2i = 0
            for r1i in range(shift, l1):
                if read1[r1i] != read2[r2i]:
                    miss_match += 1  # 不对齐位数加一
                if r2i < len(read2)-1:
                    r2i += 1
                else:
                    break
                if miss_match > max_missmatch:  # 如果不匹配位数超过容错位，跳出
                    break

            if miss_match <= max_missmatch:  # 如果不匹配位数在容错范围内，返回重叠序列长度
                return min(l1-shift, l2)  # 返回重叠序列长度  min(l1-shift, l2)
        return 0

def aligncycle(self, read1, read2, max_missmatch, cut_off):
        l1 = len(read1)
        for shift in range(0, l1-cut_off):
            miss_match = 0  # 不匹配位数
            r2i = 0
            for r1i in range(shift, l1):
                if read1[r1i] != read2[r2i]:
                    miss_match += 1  # 不对齐位数加一
                if r2i < len(read2)-1:
                    r2i += 1
                else:
                    break

                if miss_match > max_missmatch:  # 如果不匹配位数超过容错位，跳出
                    break

            if miss_match <= max_missmatch:  # 如果不匹配位数在容错范围内，返回重叠序列长度
                return l1-shift  # 返回至初始位置

        return 0

把DNA测序reads转化为TSP邻接矩阵，边的权重为overlap，边的方向为行角标→列角标

#构建邻接矩阵
tspAdjMatrix = np.zeros((n, n))
        for r1 in range(0, n):
            for r2 in range(0, n):
                if r1 != r2:
                    tspAdjMatrix[r1][r2] = align(
                        reads[r1], reads[r2], max_mismatch, cut_off)

创建变量矩阵

# 创建qubo变量矩阵。QUBO变量矩阵Reads_Sort(i,j)意义为在第i个基因read排在第j位
#创建变量
Reads_Sort = kw.qubo.ndarray((n, n), "Reads_Sort", kw.qubo.binary)
w = tspAdjMatrix

添加约束项，构建目标函数

#read约束，每个read只属于一个位置
sequence_cons = kw.qubo.quicksum(
            (1-kw.qubo.quicksum(Reads_Sort[v, j] for j in range(n)))**2 for v in range(n))

#位置约束，每个位置只能有一个read
node_cons = kw.qubo.quicksum(
            (1-kw.qubo.quicksum(Reads_Sort[v, j] for v in range(n)))**2 for j in range(n))

 # 哈密尔顿环约束
ham_cycle = sequence_cons + node_cons

# 总COST为叠加位数；DNA片段选取总长度最小的，因此要COST最大，在目标函数里取负号
path_cost = - kw.qubo.quicksum(w[u, v] * (kw.qubo.quicksum(Reads_Sort[u, j] * Reads_Sort[v, j+1]
                                                              for j in range(n-1)) + Reads_Sort[u, n-1] * Reads_Sort[v, 0])
                                       for u, v in edges)

obj = constraint * ham_cycle + path_cost

模型编译，生成Ising矩阵

# 解析QUBO
obj = kw.qubo.make(obj)
# 转化为Ising模型
obj_ising = kw.qubo.cim_ising_model(obj)
# 提取Ising矩阵
matrix = obj_ising.get_ising()["ising"]

5. 使用光量子计算模拟器进行计算

worker = kw.cim.SimulatedCIMOptimizer(
    pump=1.3,
    noise=0.2,
    laps=5000,
    delta_time=0.05,
    normalization=0.3,
    iterations=50
)

output = worker.solve(matrix)
# Sort the results
opt = kw.sampler.optimal_sampler(matrix, output, bias=0, negtail_ff=False)
# Select the best solution
cim_best = opt[0][0]
# If the linear term variable is -1, perform a flip
cim_best = cim_best * cim_best[-1]
print(cim_best)

通过模拟器求解，得到Ising解向量

[ 1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1  1]

6. 还原Ising解向量的spin结果为QUBO变量值

以上我们得到了spin结果，需要将其转化为原始变量的值，才能进一步给出求得的路径。我们使用get_val模块进行转化。

# Get the list of variable names
vars = obj_ising.get_variables()

# Substitute the spin vector and obtain the result dictionary
sol_dict = kw.qubo.get_sol_dict(cim_best, vars)

# Check the hard constraints for validity and path length
seq_val = kw.qubo.get_val(sequence_cons, sol_dict)
node_val = kw.qubo.get_val(node_cons, sol_dict)
ham_val = kw.qubo.get_val(ham_cycle, sol_dict)
print('position cons: {}'.format(seq_val))
print('node_cons cons: {}'.format(node_val))
print('ham_cycle: {}'.format(ham_val))

# Calculate the path length using path_cost
path_val = kw.qubo.get_val(path_cost, sol_dict)
print('path_cost: {}'.format(path_val))

以上我们获得了路径的合法性和长度

position cons: 0.0
node_cons cons: 0.0
ham_cycle: 0.0
path_cost: 50.0

下面我们恢复变量的值并获取最后的DNA序列。

# 获得x的数值矩阵
x_val = kw.qubo.get_array_val(Reads_Sort, sol_dict)

# 找到其中的非0项的脚标
nonzero_index = np.array(np.nonzero(x_val)).T

# 对非零项的顺序脚标进行排序
nonzero_index = nonzero_index[nonzero_index[:, 1].argsort()]

# 得到满足限制条件解,将参数记录下来
validresults = validresults + 1
print(f'DNA sequence {validresults}')
DNA sequence 4 3 1 2

我们定义一个通过解获得得节点顺序转换为DNA序列的函数

def merge_seq(reads, DNAsort):

# consensus模块，按照解顺序组装DNA。给一个序列，返回按此序列拼接的reads
merge = reads[DNAsort[0]]  # 重建基因序列
        n = len(reads)
        sum_overlap = 0  # 拼接的总重叠长度
        r1i = 0
        r2i = 0

        # 按照解的顺序连接基因片段，先从第0个连接到第n-1个
        while r2i < (n-1):
            r2i += 1
            overlap = self.align(
                reads[DNAsort[r1i]], reads[DNAsort[r2i]], 0, 0)
            overlap1 = align(reads[DNAsort[r2i-1]], reads[DNAsort[r2i]], 0, 0)
            sum_overlap = overlap1 + sum_overlap

            # 如果重叠长度大于R2长度，则R1完全覆盖R2，直接跳过，不需要融合。R1不动，R2继续下一个
            if len(self.reads[DNAsort[r2i]]) > overlap:
                merge = merge + reads[DNAsort[r2i]][overlap:]
                r1i = r2i

        # 然后连接第此时最长的尾巴个和第0个
        overlap = aligncycle(
            reads[DNAsort[n-1-r2i+r1i]], reads[DNAsort[0]], 0, 0)
        if overlap != 0:
            merge = merge[:-overlap]
        if len(reads[DNAsort[0]]) > overlap:
            sum_overlap = overlap + sum_overlap
        else:
            sum_overlap = len(reads[DNAsort[0]]) + sum_overlap
        return [merge, sum_overlap]

DNAseq = merge_seq(nonzero_index[:, 0])

量子计算与经典方法得效果对比

我们比较了在不同问题规模下（reads数分别为7、8、9、10）CIM真机和传统求解器Simulated Annealing、Tabu Search的求解时间。对比发现，在同等规模问题上，CIM真机的求解时间要远少于Simulated Annealing、Tabu Search（左图）。同时，随着问题规模增大，CIM真机求解时间相较于Simulated Annealing、Tabu Search更加稳定（右）。

图4 不同问题规模下(reads数分别为7、8、9、10)“天工量子大脑”CIM真机和传统求解器Simulated Annealing、Tabu Search的求解时间对比

以reads 数为7的问题求解为例，我们可以发现“天工量子大脑”CIM真机相较于SA算法，其哈密顿量能快速下降以求得最优解。

图5 Reads数为7时的CIM(红线)、SA(蓝线)哈密顿量演化图