CN114512178A - 基于伊辛机量子退火的密码子优化方法 - Google Patents

Abstract

本发明涉及到基于伊辛机量子退火的密码子优化方法。首先需要对mRNA上的密码子进行编码。然后再构建用于密码子优化的伊辛哈密顿量。之后需要通过量子退火伊辛机将伊辛哈密顿量演化到其最小值。最后根据量子退火伊辛机的输出结果，得到每个氨基酸对应的密码子。选择最可能的密码子的组合的过程被称之为密码子优化，本申请可以解决密码子优化方面的多项式复杂程度的非确定性问题。

Description

基于伊辛机量子退火的密码子优化方法

技术领域

本发明主要涉及到氨基酸的密码子优化技术领域，更确切的说，涉及到一种基于伊辛机量子退火的密码子优化方法。

背景技术

同义密码子是指同一种氨基酸有两个或更多密码子，称为密码子的简并性。蛋白质序列通常可以由大量各种可能的核苷酸序列予以编码。氨基酸和同义密码子之间的简并映射需要潜在核苷酸序列的数量和多肽链的长度之间的指数关系。然而，编码同一蛋白的不同核苷酸序列在表达系统当中有可能会表现出显著不同的结果。此外，最近的研究表明密码子的选择可以影响下游过程例如蛋白质折叠和功能，这对重组蛋白治疗很重要。

密码子优化又是基因表达优化的关键步骤之一，涉及到mRNA二级结构、稀有密码子和核糖体结合位点等关键因素。简单的说，基因能否顺利地表达蛋白与稀有密码子含量及与mRNA结构是否阻碍翻译有很大关系。密码子优化涉及的优化点可以从基因合成与载体构建、基因转录、mRNA翻译、翻译后修饰等过程提取，目的只有一个，就是便于相关过程的高效达成，为蛋白表达提供坚实的基础。

密码子优化是一种基于启发式评分函数而设计出的可增加基因表达的程序，已有许多评分函数被提出。虽然经典遗传算法具很高的性能，但随着多肽链长度的增长，在固定次数的迭代中采样的解空间分数呈指数级下降。因此，对于与生物学相关的用例，解决方案空间的彻底抽样通常是难以处理的。本申请的出发点之一是：解决与生物学相关的传统用例所存在的弊端，提供新的解决手段和途径。

量子信息科学和技术的最新进展已经明确表明：量子器件在某些相对狭窄的应用范围内具有超越经典器件的潜力。譬如在解决一些组合优化问题方面，特定类型的组合优化问题在最具有短期效益的量子计算领域，其相关的硬件和算法等技术发展迅速。每一种类型的量子计算技术都有其独特的优势和劣势，而且应用程序往往是量身定制的，以配合相应量子设备的优势和避免劣势。在本申请当中研究了编程量子计算机的可行性，以有效的识别具有任意目标函数评分的高质量解决方案。

发明内容

本申请公开了一种基于伊辛机量子退火的密码子优化方法，包括以下步骤：

S1、对mRNA上的密码子进行编码；

S2、构建用于所述密码子优化的伊辛哈密顿量；

S3、通过量子退火伊辛机将所述伊辛哈密顿量演化到其最小值；

S4、根据所述量子退火伊辛机的输出结果，得到每个氨基酸对应的密码子。

上述的方法，其中：

在S1中，每一个氨基酸对应有数个可能的密码子，将诸多氨基酸各自对应的数个可能的密码子全部集合在一起进行排序，设可能的密码子总数为N，则N个可能的密码子映射成N位量子位；第i个密码子的量子位标记为q_i，i＝0,1,2…,N-1。

上述的方法，其中：

每一个氨基酸自身所对应的数个可能密码子的排序是连贯的；在每一个氨基酸自身所对应的数个密码子中：只允许有一个密码子被测量后的量子位的值为1状态、其余密码子的量子位被测量后的值为相对的0状态。

上述的方法，其中：

在S2中，所述伊辛哈密顿量包括：

表征密码子使用偏好的项、优化GC含量的项、重复的核苷酸的数目的惩罚项、额外的惩罚项。

上述的方法，其中：

所述的表征密码子使用偏好的项H_f表示为：

其中c_f是设定的可调常数，q_i∈{0,1}，

是包含对数反密码子使用频率值的向量；

以及C_i表示在第i个排序位置找到密码子C的频率，ε_f是一个分母设定值。

上述的方法，其中：

所述的优化GC含量的项H_GC表示为：

其中c_GC是设定的可调常数，s_i为第i个密码子中鸟嘌呤和胞嘧啶的个数的整数及s_j表示第j个密码子中鸟嘌呤和胞嘧啶的个数的整数，ρ_T代表目标GC浓度；

第j个密码子的量子比特标记为q_j，j＝0,1,2…,N-1。

上述的方法，其中：

所述重复的核苷酸的数目的惩罚项H_R表示为：

其中c_R是设定的可调常数，r(C_i,C_j)δ_ij表征着带有对重复的核苷酸惩罚条件的矩阵；

二次函数r(C_i,C_j)返回第i个排序位置的密码子C_i和第j个排序位置的密码子C_j之间的重复序列核苷酸的最大数目；δ_ij在C_i与C_j的顺序匹配的情况下返回1，否则返回0。

上述的方法，其中：

所述的额外的惩罚项H_P表示为：

其中ε是设定的可调常数，δ′_ij在序号i和j映射到同一氨基酸位置时无穷大，否则为0。

上述的方法，其中：

所述伊辛哈密顿量H表示为：H＝H_f+H_GC+H_R+H_P。

上述的方法，其中：

所述伊辛哈密顿量H转换成其展开的等式H_BQM：

其符合二元二次方程模型的形式。

上述的方法，其中：

在S4中得到每个氨基酸对应的密码子的方式包括：

在每个氨基酸所对应的数个可能的密码子中，对应量子位为1状态的密码子被保留而对应量子态为0状态的密码子被舍弃。

本申请还公开了一种基于伊辛机量子退火的密码子优化方法，与mRNA相对应的诸多氨基酸中的每一个氨基酸均存在着多个可能性密码子，包括：

将蛋白序列的诸多氨基酸各自所对应的多个可能密码子全部集合在一起进行排序；

单独的每一个氨基酸自身所对应的多个可能密码子的排序是连贯性的；

设可能的密码子的总数为N，将N个密码子映射到N位量子位；

构建用于密码子优化的伊辛哈密顿量并通过量子退火伊辛机将其演化到最小值，从而得到每一个氨基酸的各个可能的密码子分别对应的量子位的值；

在每一个氨基酸所对应的数个可能的密码子中，对应量子位为1状态的密码子被保留而对应量子态为相对的0状态的密码子被舍弃。

本申请还公开了了一种基于伊辛机量子退火的密码子优化方法，包括步骤：

S1、基于每一个氨基酸均存在着多个可能密码子，将蛋白序列的诸多氨基酸各自所对应的多个可能密码子全部集合在一起进行排序，将N个密码子映射到N位量子位；

S2、构建一个用于密码子优化的伊辛哈密顿量；

S3、通过量子退火伊辛机将所述伊辛哈密顿量演化到其最小值；

S4、根据所述量子退火伊辛机的输出结果，从每个氨基酸所对应的数个可能的密码子中遴选出唯一一个确定的密码子，以得到mRNA上被确定下来的一系列密码子。

针对密码子优化问题，本申请提出一种全新的基于伊辛机(Ising machine)的量子退火解决方法。将多肽序列反向翻译为可表达mRNA结构是一个NP难问题，多肽序列中的每个氨基酸最多可以用6种密码子来表示，选择最可能的密码子组合的过程被称之为密码子优化，此NP(NP-Hard)难问题即多项式复杂程度的非确定性问题。本申请利用量子退火技术可开辟一条全新的、高效的能解决密码子优化问题的解决方案。绝热量子计算在确定最优解决方案方面非常具有竞争力。

附图说明

为使上述目的和特征及优点能够更加明显易懂，下面结合附图对具体实施方式做详细的阐释，阅读以下详细说明并参照以下附图之后，本申请的特征和优势将显而易见。

图1是将诸多氨基酸各自所对应的多个密码子全部集合在一起进行排序。

图2是甘氨酸氨基酸的每个可能密码子映射到量子位元的示范性示意图。

图3是构造惩罚矩阵以对多个密码子量子位同时处于1的情况进行惩罚。

图4是在氨基酸对应的数个可能密码子中保留对应量子位为1的密码子。

图5是通过量子退火伊辛机将伊辛哈密顿量演化到其最小值而得到结果。

图6是根据量子退火伊辛机的输出得到每个氨基酸对应的密码子的过程。

图7是对应量子位为1的密码子被保留以及量子位为0的密码子被舍弃。

具体实施方式

下面将结合各实施例，对本发明的方案进行清楚完整的阐述，所描述的实施例仅是本发明用作叙述说明所用的实施例而非全部的实施例，基于该等实施例，本领域的技术人员在没有做出创造性劳动的前提下所获得的方案属于本发明的保护范围。

关于密码子优化方面：是一种通过增加靶基因的翻译效率来提高生物体内蛋白质表达水平的新技术。通常会通过避免稀有密码子，利用偏爱密码子(Preferred codons)、以及通过简化mRNA二级结构、优化重复的序列、消除限制酶切位点、调整GC含量等方法来重新设计基因，以提高翻译效率，进而提高蛋白表达水平。编码蛋白质的密码子可能会影响到翻译效率，进而影响蛋白质表达。实际上，研究表明：翻译效率与大肠杆菌和酿酒酵母中所有内源基因的密码子偏好相关。

关于密码子的偏好性：是指不同生物、甚至同种生物不同的蛋白质编码基因，对简并密码子使用频率并不相同。譬如举例而言，大肠杆菌中某个氨基酸，编码该氨基酸的密码子有两个如AGG和ATG，但是大肠杆菌却较偏爱利用AGG去指导氨基酸的合成。

密码子优化的必要性：细胞内对应于稀有密码子的tRNA较少，高频率使用这些密码子的基因翻译过程易受阻，大大降低蛋白表达水平。若是在基因的同义密码子使用频率与表达宿主相匹配的情况下，蛋白的表达水平则会显著提高。1980年代后期，科学家根据高表达基因的参考文献中的密码子使用频率创建了密码子适应指数(CAI)，通常情况下若出现CAI＜0.80的情况，普遍被认为需进行密码子优化。

关于GC含量：GC含量是影响表达的另一个重要因素，相比较AT间的2个氢键的相对键能而言，GC间存在着3个氢键，所以GC直接影响DNA的稳定性。由于GC配对能形成3个氢键而且比AU配对或AT配对都来的更加稳定，因此GC含量不仅仅是直接的影响PCR退火温度，也影响着mRNA二级结构的稳定性，还间接关系到mRNA翻译成蛋白质的效率。过低的GC含量可能会导致转录终止、过高的GC含量易使mRNA形成稳定的二级结构并阻碍蛋白表达。因此优化GC含量自然存在着必要性。

由于本文与量子相关，关于量子器件和量子数据的相关内容如下文所述。

本文所谓“量子器件”包含已知的量子计算设备和量子芯片等，也可以用量子硬件替代量子器件这类术语。典型的“量子器件”包括但不限制于：量子计算机、量子信息处理系统或量子密码系统、量子模拟器、处理量子数据的所有种类的装置、设备和机器。

本文所谓“量子数据”包含由量子系统携带、保存或存储的信息或数据，最小的非平凡系统是量子比特，即定义量子信息单位的系统。应当理解，术语“量子比特”包括在相应上下文中可以适当地近似为二能级系统的所有的量子系统。这种量子系统举例来说通常包括了典型的原子、电子、光子、离子或超导量子比特等等。

基于前文介绍的内容，本申请涉及到了密码子优化，密码子优化可允许在量子器件上进行的BQM(二元二次模型)运算，其哈密顿量的设计至少可以满足以下目的：

其一、优化GC含量。

其二、最小化顺序重复的核苷酸。

其三、优化密码子使用偏好。

参见图1，每一个氨基酸对应有数个密码子，换言之就是是蛋白序列中每个氨基酸可以用多种密码子来表示。首先需要将蛋白序列中的诸多氨基酸各自所对应的数个密码子全部集合在一起进行排序，设密码子总数为N，则N个密码子映射成N位量子比特并且可以定义第i个密码子的量子比特标记为q_i，其中i＝0,1,2…,N-1。

参见图1，暂以甘氨酸G、丝氨酸S、缬氨酸V、亮氨酸L和赖氨酸K等作为氨基酸代表来展示排序的过程。但值得注意的是，给出的有限氨基酸代表数仅只是用于对本申请进行阐释说明，蛋白序列上的氨基酸的实际数量通常比图示的代表数多得多，尤其是图示的氨基酸代表之间的排列组合关系也仅仅只是作为范例却不构成任何特定限制。

参见图4，甘氨酸G的第一至第四种可能密码子分别是GGU、GGC、GGA、GGG。

参见图4，丝氨酸S的第一至第六种可能密码子分别是UCU、UCC、UCA、UCG以及剩下的可能密码子分别是AGU、AGC。

参见图4，赖氨酸K的第一至第二种可能密码子分别是AAA、AAG。

参见图1，缬氨酸V的第一至第四种可能密码子分别是GUU、GUC、GUA、GUG。

参见图1，亮氨酸L以及更多的氨基酸的可能密码子不再一一赘述，属现有技术。

参见图1，需要将诸多氨基酸{G、S、V、L、K…}等各自所对应的可能密码子全部集合在一起进行排序。假设从甘氨酸G的四个可能密码子开始进行排序，譬如在所有密码子的全部集合之中，甘氨酸G的四种可能密码子排在第0～3位，丝氨酸S的六种可能密码子排在第4～9位，缬氨酸V的四种可能密码子排在第10～13位，…，赖氨酸K的两种可能密码子排在第x位和第x+1位。以此类推，密码子总数为N就有N位排序。注意这里出现的总数N是正整数，作为范例的x也是正整数。

参见图1，限于篇幅，只给出了少量密码子作为排序(SORT)的阐释示范。本文主旨是基于伊辛机量子退火，需要将N个密码子映射成N位量子比特(QUBIT)。

参见图1，甘氨酸G的四种可能密码子映射到量子比特q₀、q₁、q₂、q₃，第一个可能密码子映射到量子比特q₀、第二个可能密码子映射到量子比特q₂、第三个可能密码子映射到量子比特q₃、第四个可能密码子映射到量子比特q₄。

参见图1，丝氨酸S的六种密码子映射到量子比特q₄、q₅、q₆、q₇、q₈、q₉，第一个可能密码子至第六个可能密码子相应分别映射到量子比特q₄～q₉。

参见图1，缬氨酸V的四种密码子映射到量子比特q₁₀、q₁₁、q₁₂、q₁₃，第一个可能密码子至第四个可能密码子相应分别映射到量子比特q₁₀～q₁₃。

参见图1，按照密码子与量子比特之间的映射关系，合计有N位量子比特。第i个密码子的量子比特标记为q_i，其中自然数i满足i＝0,1,2…,N-1。

参见图1，每一个氨基酸自身所对应的数个密码子的排序是连贯的。甘氨酸G所对应的四个密码子的排序是连贯的，例如排序在N个密码子中的第0～3位。丝氨酸S所对应的六个密码子的排序是连贯的，例如排序在N个密码子中的第4～9位。缬氨酸V所对应的四个密码子的排序是连贯的，例如排序在N个密码子中的第10～13位。换句话说就是任何某一种氨基酸所对应的任意一个密码子不会排序穿插到其他氨基酸所对应的数个密码子的排序序号之间，例如甘氨酸G对应的任意一个密码子不会穿插到丝氨酸S所对应的数个密码子的排序序号4～9之间。

参见图1，鉴于将诸多氨基酸例如{G、S、V、L、K…}等它们各自所对应的数个密码子全部集合在一起进行了排序，密码子的排序规则也反映了与它们对应的一系列量子比特的排序关系。甘氨酸G的四个量子位按序连贯排列q₀、q₁、q₂、q₃之后是丝氨酸S的六个量子位按序连贯排列q₄、q₅、q₆、q₇、q₈、q₉，如图所示，再之后则是缬氨酸V的四个量子位也按序连贯排列q₁₀、q₁₁、q₁₂、q₁₃等等，诸如此类。因此，与每一个氨基酸相关的数个可能密码子所对应的数个量子比特或量子位的排序也是连贯的。

参见图1，注意量子设备使用量子位元来存储数据，其相关的数据在测量时可被数字解码为0和1，但在计算过程中也可能是0和1的叠加态。为了将经典的遗传数据编码到量子设备当中，值得强调的是，每一个可以映射到目标多肽序列的密码子都需要用物理量子位明确表示。测量后返回“1”的量子位代表多肽序列中每个位置选择的密码子。因此多肽序列中每个位置只能有1个量子位(密码子)处于“1”状态或二进制1，其余量子位测量后必须返回“0”状态或或二进制0，这里的位置是指氨基酸的位置。

参见图2，譬如甘氨酸G的四个量子位测量后，只允许有1个量子位(例如第三个可能密码子)处于“1”的状态，其余的三个量子位测量后必须返回“0”状态。在数据表现形式上可将其描述为q₀～q₃＝0010。

参见图4，譬如丝氨酸S的六个量子位测量后，只允许有1个量子位(例如第四个可能密码子)处于“1”的状态，其余的五个量子位测量后必须返回“0”状态。在数据表现形式上可将其描述为q₄～q₉＝000100。

参见图4，譬如赖氨酸K的两个量子位测量后，只允许有1个量子位(例如第二个可能密码子)处于“1”的状态，其余的一个量子位测量后必须返回“0”状态。在数据表现形式上可将其描述为q₁₀～q₁₁＝01。

参见图3，前文已告知每个氨基酸位置只能有1个量子位处于“1”状态，这种方案是通过构建一个二维惩罚矩阵来实现的，它会为映射到多肽序列中相同位置的密码子对增加无限的能量。例如通过构造惩罚矩阵：对多个密码子处于“1”状态的情况进行惩罚而增加无限能量。在例子中量子位q₂处于“1”状态，其余量子位为“1”状态，这种情况下将返回等于0的能量惩罚。但如果q₂处于“1”及例如q₁也处于“1”状态，这种情况下譬如可返回q₁等于∞的能量惩罚。能量惩罚的相关内容在下文会继续介绍。

参见图4，是将某段示意性的蛋白序列100的局部某些氨基酸{GSK…}的密码子映射到量子位元的例子。为蛋白序列中的每个氨基酸位置选择一个密码子：甘氨酸G的四种可能性密码子中选择了唯一一个量子位q₂所对应的密码子GGA、丝氨酸S的六种可能密码子中选择了唯一一个量子位q₇所对应的密码子UCG、赖氨酸K的两种可能性密码子中选择了唯一一个量子位q₁₁所对应的密码子AAG。依此类推。当为蛋白序列中的每个氨基酸位置选定好一个明确的密码子之后，即为每个氨基酸位置分配优化的密码子，籍此将多肽序列反向翻译为可表达mRNA结构之目的便可达成。为每个氨基酸位置分配唯一最佳的密码子的任务与伊辛机量子退火紧密相关。

参见图5，基于伊辛机量子退火的密码子优化方法，包括以下各步骤的内容。

参见图5，步骤S1、对mRNA上的密码子进行编码。

参见图5，步骤S2、构建用于密码子优化的伊辛哈密顿量H_Ising。

参见图5，步骤S3、量子退火伊辛机将伊辛哈密顿量演化到最小值min(H_Ising)。

参见图5，步骤S4、由量子退火伊辛机的输出结果，得到氨基酸对应的密码子。

参见图5，关于对mRNA上的密码子进行编码在前文已经阐释过：将氨基酸对应的可能密码子进行排序，例如对氨基酸对应的n个密码子进行排序，序号对应QUBIT而且规定这里的n个QUBIT有且仅有一个能取“1”状态。图4是以n＝4,6,2为例。

参见图5，构建伊辛哈密顿量将从“嵌入密码子使用的偏好、优化目标的GC含量及重复的核苷酸的数目的惩罚项、额外的惩罚项”等若干方面进行诠释(涉及S2)。

第一方面，表征密码子使用偏好的项H_f。

关于嵌入密码子使用偏好。密码子的使用频率因宿主系统而异因此，可对评分函数进行了定制，以匹配表达式系统。在本申请的相关研究中，作为研究工具，大肠杆菌的密码子使用频率从python-codon-tables0.1.10库中予以输入。可设定C_i表示在位置i找到密码子C(是codon的缩写，注意与GC碱基相关的字母标识进行区分)的频率。因此电位被设计为对稀有密码子返回较大的代价(此情况下C_i很小)，作为相反状况，而对系统容易获得的密码子产生可忽略不计的代价(此情况下C_i很大)。在可选范例中，其相关的函数是逆的对数乘以-1。这个函数产生了向稀有密码子添加较大惩罚和向可访问密码子添加较小惩罚的理想行为。然而，该函数在C_i＝0处没有定义。如果宿主系统确实不能访问给定的密码子，那么包含该密码子的任何序列都是不可表达的，表达概率为零。然而评分函数须限制为有限的十进制值，因此在函数的分母上可加上一个无限小的值ε_f，以避免未定义的值和避免分母为零的值。对于包含N个可能密码子的系统，哈密顿量H_f如下所述。

哈密顿量H_f关系式中，c_f是一个可调常数，q_i∈{0,1}是量子位的值，

是一个包含对数反密码子使用频率值的向量。给定二元q_i值，哈密顿量只惩罚被“选择”的密码子而且惩罚的对象通常是利用q_i＝1的量子位元表示。

函数涉及到密码子和反密码子等概念，密码子(codon)指信使RNA分子中每相邻的三个核苷酸编成一组，在蛋白质合成时，代表某一种氨基酸的规律。另外上述函数涉及到反密码子(anti-codon)：RNA链经折叠而看上去像三叶草的叶形，其一端是携带氨基酸的部位以及相对的另一端有3个碱基，每个tRNA(transferRNA)的这3个碱基的规律是可以与mRNA上的密码子互补配对，因而称之为反密码子。

第二方面，优化GC含量的项H_GC。

关于优化目标的GC含量。针对优化对象：核苷酸序列的GC浓度，必须引入一个与其浓度相关的式子ρ_GC以及引入一个成本函数Δ。若给定了目标GC浓度，那么最终的主要目的是需要最小化ρ_GC与目标GC浓度即ρ_t之间的差异。而满足这个约束条件的最简单的目标函数则是如下所述的一个二次函数Δ。

Δ＝c_GC(ρ_GC-ρ_T)²

注意二次函数Δ关系式中c_GC是一个可调常数。除了可调常数c_GC之外，本申请在后续内容中将会进一步给出包括c_f和c_R及ε和ρ_T等各类可调常数的可选值范例。所谓可调常数可认为是事先设计或设定的常数值而且数值的大小可调，如依经验而定。GC含量的计算方法是将长度为N的序列中G和C(这里的GC是碱基标识)的个数相加，将序列中核苷酸的个数归一化。具体计算的表达式如下。

以上是与GC浓度相关的浓度式子ρ_GC的表达式并引入了s_i的概念。

其中s_i为表示密码子i(排序为第i位的密码子)中G和C的个数的整数，或换而言之可认为s_i为第i个密码子中鸟嘌呤和胞嘧啶的个数的整数、s_j表示第j个密码子中鸟嘌呤和胞嘧啶的个数的整数。q_i为第i个的密码子所对应的量子位的值。

将式Δ函数的表达式展开可得到类似于二元二次模型(BQM)公式的形式。

双求和中矩阵必须要是一个上三角矩阵以满足QUBO(二次无约束二进制优化)问题的矩阵要求，将其分解成迹和非对角元素部分，二者的和被限制为上三角矩阵。迹部分应当符合一定的条件：迹需要

上的单个求和并且

Δ函数的局部内容如下所述。

由于矩阵是对称的，所有非对角项都以上三角形式乘以2。因此前述H_GC函数展开之后为下文所展现的形式。

用与BQM(如下文的伊辛哈密量H)一致的形式概括二次代价函数(如Δ函数)。

第三方面，重复的核苷酸的数目的惩罚项H_R。

关于重复的核苷酸的数目的惩罚项。为减少序列中重复核苷酸的数量，对氨基酸序列中所有序列位置(如i＝0,1,2…,N-1等合计N个序号处)的密码子进行比较。在本申请中构造二次函数r(C_i,C_j)，二次函数r(C_i,C_j)之目的是返回密码子C_i和C_j之间重复序列核苷酸的最大数目。第i个排序位置的密码子是C_i和第j个排序位置的密码子是C_j。该函数的另一个作用是，在空情况下通过减去1移到原点。

用数学语言描述：记录同时出现密码子C_i和C_j两者中的、重复的碱基的次数，同时还规定重复的碱基出现在的C_i的结尾点或C_j的起始点，则所述二次函数r(C_i,C_j)的值等于该次数的平方再减去1。若是假设重复次数为z那么二次函数返回z²-1。

例如有r(ATA,TCG)＝1²-1＝0。是以C_i＝ATA和C_j＝TCG为例。发现重复的碱基此时处于空情况下。那么通过减去1移到原点。即规定r(C_i,C_j)返回0。

例如有r(ATA,ACG)＝2²-1＝3。是以C_i＝ATA和C_j＝ACG为例。发现重复的碱基出现在的C_i的结尾点(即C_i的结尾A)和C_j的起始点(即C_j的起始A)。

例如有r(CGG,GGG)＝5²-1＝24。是以C_i＝CGG和C_j＝GGG为例。则重复的碱基出现在的C_i的结尾点(即C_i的两个G)和C_j的起始点(即C_j的三个G)。

二次函数r(C_i,C_j)返回第i个排序位置的密码子C_i和第j个排序位置的密码子C_j之间的重复序列核苷酸的最大数目，在上述举例部分进行了阐释。如果C_i和C_j之间的重复序列核苷酸的重复次数为z那么二次函数返回z²-1。特殊情况，如果C_i和C_j之间的重复序列核苷酸的重复次数为空那么二次函数返回0。

二次函数r(C_i,C_j)的设计是为实现带有惩罚功能的矩阵R_ij作铺垫。

重复的核苷酸惩罚存储在一个设计的矩阵R_ij中，矩阵R_ij与二次函数r(C_i,C_j)相关联的表达式如下所述，矩阵R_ij还与一个函数δ_ij相关联。

R_ij＝r(C_i,C_j)δ_ij

矩阵的Delta函数δ_ij在密码子i与j(即第i个密码子C_i和第j个密码子C_j)的顺序匹配的情况下返回1，否则返回0(若密码子C_i和C_j的顺序不匹配)。

给定核苷酸序列的总重复核苷酸惩罚项H_R。

由于矩阵R_ij为上三角矩阵，所以允许在不改变结果的情况之下，可以将成对和限制为上三角元素，使其与BQM模型兼容(如下文记载的伊辛哈密量H)。惩罚项H_R的原关系式中额外引入了一个可调常数c_R来衡量哈密顿量的贡献，得到如下的H_R形式。

第四方面，额外的惩罚项H_P。

关于额外的惩罚项。前述描述的能量术语是为了给特定的序列属性添加惩罚。这类目标函数的基态能量为零，因为引入密码子增加了分数。为了抵消这一趋势，在每个密码子的单体项中减去一个常数ε，从而增加了密码子引入系统的能量优势。该常数的绝对值必须超过单体相互作用h的最大值(参阅以下方程式)。

在伊辛哈密量H公式中，第i个密码子的量子比特标记为q_i，第j个密码子的量子比特标记为q_j，第k个密码子的量子比特标记为q_k，排序序号i,j,k＝0,1,2…,N-1。其中量子比特的值对于二元表示可以是{0,1}或对于自旋表示可以是{-1,1}。J_ik是双体项。对于伊辛模型来说其公式中的h代表了物理的自旋而J则代表自旋之间相互作用的能量。

为了否定为一个给定位置(例如在诸多氨基酸中任意一个给定的氨基酸处)分配多个密码子的可能性，引入了一个位点特定的delta函数δ′_ij，该函数δ′_ij对分配到同一给定位置的密码子对施加有效的无限惩罚。值得注意的是，这里所言的“否定为一个给定位置分配多个密码子的可能性”，是指需要从一个给定位置的多种可能密码子中选择唯一一个确定下来的密码子来决定这一个给定位置处的氨基酸，而不是指代“不允许某一个给定位置使用多个可能密码子来表示”的意思。

关于函数δ′_ij的第一种取值①情况：函数δ′_ij在序号i和j映射到诸多氨基酸中的同一氨基酸位置时无穷大。可阻止为任意某个氨基酸位置分配多个密码子(例如不允许第i个密码子C_i和第j个密码子C_j同时映射到同一氨基酸位置处)。

关于函数δ′_ij的第二种取值②情况：函数δ′_ij在序号i和j未映射到诸多氨基酸的同一氨基酸位置时等于零。第一种取值是∞的情况，第二种取值是0的情况。

对于一个含有N个可能密码子的系统而言，哈密顿量通过添加以下项H_P来修正。

哈密顿量添加的项H_P中的ε是事先设定的可调常数。

最后还需要构建目标函数。哈密顿量代表了一个核苷酸序列的总的能量，将之前的公式加起来得到以下改进的哈密顿量H表达式(与前文描述的H略有差异)。显而易见改进的哈密顿量综合的考虑了“表征密码子使用偏好的项、优化GC含量的项、重复的核苷酸的数目的惩罚项、额外的惩罚项”等内容。

H＝H_f+H_GC+H_R+H_P

哈密顿量H扩展并整理公式后得到H_BQM。

伊辛哈密顿量H转换成其展开的等式H_BQM的形式与BQM的形式是一致的，并且能直接在可解决BQM运算的框架中实现。

本申请上文涉及到各个方程式中的一些常量的可选取值范例如下。

可调常数c_f在一个可选但非限制性的示例中的取值等于0.1。

可调常数c_GC在一个可选但非限制性的示例中的取值等于1。

可调常数c_R在一个可选但非限制性的示例中的取值等于0.1。

可调常数ε在一个可选但非限制性的示例中的取值等于1。

可调常数ρ_T在一个可选但非限制性的示例中的取值等于0.5。

参见图6，哈密顿量H或H_BQM经过通过量子退火伊辛机(Isingmachine)的量子退火机制可将伊辛哈密顿量演化到其最小值min(H_Ising)，此时用H_Ising或(H_BQM)_Ising表达式作为所构建的用于密码子优化目的的伊辛哈密顿量。

参见图6，S1是编码过程。S2、S3是伊辛机进行量子退火演化的过程。S4是由量子退火伊辛机的输出结果来得到氨基酸对应的密码子的过程。

参见图6，关于量子退火演化，可借助于现有的D-Wave系统来解释：量子比特是超导回路的最低能量状态，此情况下的超导回路是D-WaveQPU的组成部分。这些能量状态有循环电流和相应磁场。同经典比特一样，量子比特的状态可以是0或1，因为量子比特是量子对象，所以它允许同时具有0态和1态的叠加态。量子退火过程结束时每个量子比特从叠加态坍缩成0或1(经典状态)，如S2、S3演化得到S4的结果。

参见图6，量子退火伊辛机的输出结果为{q₀,q₁,q₂,q₃…q_N-1}。全部输出结果中的量子位结果q[3:0]＝0010是甘氨酸G的四个可能性密码子的测量值，全部输出结果中的量子位结果q[9:4]＝000100是丝氨酸S的六个可能密码子的测量值，全部输出结果中的量子位结果q[13:10]＝0100是缬氨酸V的四个可能密码子测量值。全部输出结果中的量子位结果q[x+1:x]＝01是赖氨酸K的两个可能密码子的测量值。步骤S4的输出部分展现了有限的若干氨基酸的量子位退火结果，其他量子位不再一一赘述。

参见图6，甘氨酸G的四个可能密码子(对应的量子比特为q₀、q₁、q₂、q₃)中选择第三个可能密码子(对应的量子比特为q₂)作为最终被确定下来的密码子。换而言之只允许有一个密码子的量子态为1、其余密码子的量子态为0。在甘氨酸G所对应的四个可能性密码子中，对应量子态为1的密码子被保留而对应量子态为0的密码子被舍弃。

参见图6，丝氨酸S的六个可能密码子(量子位为q₄、q₅、q₆、q₇、q₈、q₉)中选择第四个可能密码子(对应的量子比特为q₇)作为最终被确定下来的密码子。换而言之只允许有一个密码子的量子态为1、其余密码子的量子态为0。在丝氨酸S所对应的六个可能性密码子中，对应量子态为1的密码子被保留而对应量子态为0的密码子被舍弃。

参见图6，缬氨酸V的四个可能密码子(对应量子位为q₁₀、q₁₁、q₁₂、q₁₃)中选择第二个可能密码子(对应的量子比特为q₁₁)作为最终确定下来的密码子。换而言之只允许有一个密码子的量子态为1、其余密码子的量子态为0。在缬氨酸V所对应的四个可能性密码子中，对应量子态为1的密码子被保留而对应量子态为0的密码子被舍弃。

参见图6，赖氨酸K的两个可能性密码子(对应的量子比特分别为q_x、q_x+1)中选择第二个可能密码子(对应的量子比特为q_x+1)作为最终确定下来的密码子。换而言之只允许有一个密码子的量子态为1、其余密码子的量子态为0。在赖氨酸K所对应的两个可能性密码子中，对应量子态为1的密码子被保留而对应量子态为0的密码子被舍弃。

参见图6，甘氨酸G的最终密码子(量子比特为q_G)选择第三个可能密码子。注意在退火程序中甘氨酸G可能会选择任意一个可能密码子，即q_G可能是q₀～q₃中的任意一者而不仅仅只限制于选取的q₂，选取q₂只是可能结果中的一种。根据量子退火伊辛机的输出结果可以得到甘氨酸G对应的密码子，即量子比特为q_G的可能密码子，视为mRNA上被确定下来的一个密码子(mRNA-G)。

参见图6，丝氨酸S的最终密码子(量子比特为q_S)选择第四个可能密码子。注意在退火程序中丝氨酸S可能会选择任意一个可能密码子，即q_S可能是q₄～q₉中的任意一者而不仅仅只限制于选取的q₇，选取q₇只是可能结果中的一种。根据量子退火伊辛机的输出结果可以得到丝氨酸S对应的密码子，即量子比特为q_S的可能密码子，视为mRNA上被确定下来的一个密码子(mRNA-S)。

参见图6，缬氨酸V的最终密码子(量子比特为q_V)选择第二个可能密码子。注意在退火程序中缬氨酸V可能会选择任意一个可能密码子，即q_V可能是q₁₀～q₁₃中任意一者而不仅仅只限制于选取的q₁₁，选取q₁₁只是可能结果中的一种。根据量子退火伊辛机的输出结果可得到缬氨酸V对应的密码子，即量子比特为q_V的可能密码子，视为mRNA上被确定下来的一个密码子(mRNA-V)。

参见图6，赖氨酸K的最终密码子(量子比特为q_K)选择第二个可能密码子。注意在退火程序中赖氨酸K可能会选择任意一个可能密码子，即q_K可能是q_x～q_x+1中任意一者而不仅仅只限制于选取的q_x+1，选取q_x+1只是可能结果中的一种。根据量子退火伊辛机的输出结果得到赖氨酸K对应的密码子，即量子比特为q_K的可能密码子，视为mRNA上被确定下来的一个密码子(mRNA-K)。

参见图6，展现了将多肽序列反向翻译为可表达mRNA结构的大致流程。

参见图7，最终mRNA上对应的所有密码子都被确定下来。

参见图7，例如mRNA上被确定下来的密码子是：与q_G、q_S、q_V、…q_K…等量子比特对应的密码子。例如由挑选出来的与q₂、q₇、q₁₁、…q_x+1…等量子比特相对应的密码子作为决定各个氨基酸的最终密码子，挑选出来的量子比特的量子位为1。

综上所述，本申请的核心是从每个氨基酸的可能密码子(possiblecodon)中筛选或曰优化出最佳的密码子，从优化GC浓度、最小化顺序重复的核苷酸、优化密码子使用偏好等方面着手，得到mRNA上对应的所有最终被确定下来的密码子。与现有技术相比本申请的优势是能够以较高正确率、较高处理速度的方式遴选出最优的密码子。

以上通过说明和附图的内容，给出了具体实施方式的特定结构的典型实施例，上述申请内容提出了现有的较佳实施例，但这些内容并不作为局限。对于本领域的技术人员而言在阅读上述说明后，各种变化和修正无疑将显而易见。因此，所附的权利要求书应当看作是涵盖本发明的真实意图和范围的全部变化和修正。在权利要求书范围之内的任何和所有等价的范围与内容，都应认为仍属本发明的意图和范围内。

Claims (12)

1.一种基于伊辛机量子退火的密码子优化方法，其特征在于，包括以下步骤：

S1、对mRNA上的密码子进行编码；

S2、构建用于密码子优化的伊辛哈密顿量；

S3、通过量子退火伊辛机将所述伊辛哈密顿量演化到其最小值；

S4、根据所述量子退火伊辛机的输出结果，得到每个氨基酸对应的密码子。

2.根据权利要求1所述的方法，其特征在于：

在S1中，每一个氨基酸对应有数个可能的密码子，将诸多氨基酸各自对应的数个可能的密码子全部集合在一起进行排序，设可能的密码子总数为N，则N个可能的密码子映射成N位量子位；第i个密码子的量子位标记为q_i，i＝0,1,2…,N-1。

3.根据权利要求2所述的方法，其特征在于：

每一个氨基酸自身所对应的数个可能密码子的排序是连贯的；在每一个氨基酸自身所对应的数个密码子中：只允许有一个密码子被测量后的量子位的值为1状态、其余密码子的量子位被测量后的值为相对的0状态。

4.根据权利要求2所述的方法，其特征在于：

在S2中，所述伊辛哈密顿量包括：

表征密码子使用偏好的项、优化GC含量的项、重复的核苷酸的数目的惩罚项、额外的惩罚项。

5.根据权利要求4所述的方法，其特征在于：

所述的表征密码子使用偏好的项H_f表示为：

其中c_f是设定的可调常数，q_i∈{0,1}，

是包含对数反密码子使用频率值的向量；

以及C_i表示在第i个排序位置找到密码子C的频率，ε_f是一个分母设定值。

6.根据权利要求5所述的方法，其特征在于：

所述的优化GC含量的项H_GC表示为：

其中c_GC是设定的可调常数，s_i为第i个密码子中鸟嘌呤和胞嘧啶的个数的整数及s_j表示第j个密码子中鸟嘌呤和胞嘧啶的个数的整数，ρ_T代表目标GC浓度；

第j个密码子的量子比特标记为q_j，j＝0,1,2…,N-1。

7.根据权利要求6所述的方法，其特征在于：

所述重复的核苷酸的数目的惩罚项H_R表示为：

其中c_R是设定的可调常数，r(C_i,C_j)δ_ij表征着带有对重复的核苷酸惩罚条件的矩阵；

二次函数r(C_i,C_j)返回第i个排序位置的密码子C_i和第j个排序位置的密码子C_j之间的重复序列核苷酸的最大数目；δ_ij在C_i与C_j的顺序匹配的情况下返回1，否则返回0。

8.根据权利要求7所述的方法，其特征在于：

所述的额外的惩罚项H_P表示为：

其中ε是设定的可调常数，δ′_ij在序号i和j映射到同一氨基酸位置时无穷大，否则为0。

9.根据权利要求8所述的方法，其特征在于：

所述伊辛哈密顿量H表示为：H＝H_f+H_GC+H_R+H_P。

10.根据权利要求9所述的方法，其特征在于：

所述伊辛哈密顿量H转换成其展开的等式H_BQM：

其符合二元二次方程模型的形式。

11.根据权利要求3所述的方法，其特征在于：

在S4中得到每个氨基酸对应的密码子的方式包括：

在每个氨基酸所对应的数个可能的密码子中，对应量子位为1状态的密码子被保留而对应量子态为0状态的密码子被舍弃。

12.一种基于伊辛机量子退火的密码子优化方法，其中每一个氨基酸均存在着多种可能的密码子，其特征在于，该方法包括：

将蛋白序列的诸多氨基酸各自所对应的多个可能密码子全部集合在一起进行排序；

单独的每一个氨基酸自身所对应的多个可能密码子的排序是连贯性的；

设可能的密码子的总数为N，将N个密码子映射到N位量子位；

构建用于密码子优化的伊辛哈密顿量并通过量子退火伊辛机将其演化到最小值，从而得到每一个氨基酸的各个可能的密码子分别对应的量子位的值；

在每一个氨基酸的数个可能的密码子中，相应量子位为1状态的密码子被保留而相应量子位为0状态的密码子被舍弃，以得到mRNA上被确定下来的一系列密码子。

Images