引言

生物体内很多重要的生理过程，如细胞的周期调控、合成代谢、信号转导以及遗传信息的传递等，都有赖于蛋白质和配体分子的相互作用和识别^[1-3]。分子对接是研究蛋白质与配体分子间相互作用与识别的理论方法，对生命活动的分子机理研究、生物分子复合物结构预测和靶向药物筛选等都具有重要的意义^[4-6]。其中，创新药物开发是分子对接的重要应用领域。创新药物研究具有重大的社会效益和经济效益，而新药开发却是一件耗资巨大且效率较低的工作，一般平均要筛选十万种以上的化合物才能得到一种新药。目前，已知的化合物分子数量已达到数千万量级，使药物开发具备了大数据的基础。但依靠传统药理学和实验模型测试成千上万个化合物，开发的时间成本和资金成本巨大。分子对接则可以快速高效地对千万量级的化合物分子进行计算筛选，精准发现有潜力的备选新分子，因此可以有效推进新药研发的进程。

分子对接被用于研究分子间相互作用，经典的热力学认为, 蛋白质和配体分子间的相互作用是系统热力学平衡的过程, 相互作用所形成的复合物结构应该是结合自由能最低的构象^[7]。因此，一方面分子对接计算需要采用数学模型或函数尽量准确的计算结合自由能。另一方面，还需要发展有效的搜索算法快速找到自由能极低的构象。分子对接中的构象搜索是极其复杂的问题。假设受体结构固定不动，配体相对受体运动，在刚性对接的前提下，至少需要搜索平移和旋转6个自由度。如果进一步考虑受体或配体的分子柔性，则搜索空间将更复杂。

蛋白质-配体分子对接构象搜索方法主要包括快速穷举搜索和启发式搜索两大类。如果配体为生物大分子，如蛋白质、RNA或DNA，则相互作用的区域可以出现在分子表面任何位置(图 1(a))，因此往往需要进行全局搜索，可采用快速穷举搜索遍历各种位置，或采用启发式算法进行近似的全局搜索。如果配体为小分子化合物，分子对接则往往有一定的结合口袋(图 1(b))，可以限定构象空间的搜索范围，通常会采用启发式算法进行限定区域的搜索。

本文将结合国内外研究进展以及本课题组开展的一些研究工作, 就蛋白质和配体分子对接中构象搜索方法的进展情况进行分析和综述。

1 快速穷举搜索方法

穷举搜索是当处理的问题找不到明确解决方法和规律时，逐一枚举此问题的各种可能情况，通过一定的评判原则从中挑选出符合要求的情况作为问题的候选解。当配体为生物大分子时，蛋白质受体表面的任何区域都有可能成为相互作用界面，因此，尽可能穷举所有构象并进行评价成为蛋白质-配体分子对接中的常用的搜索方法。

1.1 快速傅里叶变换方法

蛋白质-配体分子对接穷举搜索中，应用最广泛的是Katchalski-Katzir等^[8]提出的快速傅里叶变换(Fast Fourier transform, FFT)方法。该方法首先将蛋白质分子和配体分子表示为三维格点数据，评判原则如分子几何或能量的匹配程度，采用相关函数来定量地描述，具体评判公式可表示为

$ {\mathit{\boldsymbol{c}}_{\alpha , \beta , \gamma }} = \sum\limits_{l = 1}^N {\sum\limits_{m = 1}^N {\sum\limits_{n = 1}^N {} } } {\mathit{\boldsymbol{a}}_{l, m, {\rm{ }}n}}\cdot{\mathit{\boldsymbol{b}}_{l + \alpha , m + \beta , n + \gamma }} $

(1)

式中：a和b分别是蛋白质a和配体b的三维格点数据，可以表示几何、静电或统计势等各种格点形式。以静电计算为例，假如a为蛋白质a附近的静电场格点，b为配体b的静电荷分布格点，则相关函数c表示两者之间的静电相互作用能。其中α，β，γ分别是笛卡尔坐标系下，配体b质心在3个方向上相对蛋白质a质心平移的格点数，即(α，β，γ)为配体b相对蛋白质a的平移向量。N是三维格点最大的格点值。直接计算相关函数的计算量较大，以式(1)为例，乘法计算量正比于N³，如果α，β，γ遍历所有的格点位置，则整个平移搜索的计算量为N⁶。相关函数的计算可以采用FFT方法来加速，具体的计算步骤如下：

(1) 对蛋白质三维格点a进行快速傅里叶变换并求复共轭操作

$ \mathit{\boldsymbol{A}}_{_{o, p, q}}^ * = {\{ {\rm{FFT}}({\mathit{\boldsymbol{a}}_{l, m, n}})\} ^*} $

(2)

式中:o, p, q={1, 2, …, N}，*表示复共轭，快速傅里叶变换的乘法计算量是正比于N³ln(N³)。

(2) 对配体格点b进行快速傅里叶变换

$ {\mathit{\boldsymbol{B}}_{o, p, q}} = {\rm{FFT}}({\mathit{\boldsymbol{b}}_{l, m, n}}) $

(3)

(3) 对A_{o, p, q}^*和B_{o, p, q}进行点乘计算

$ {\mathit{\boldsymbol{C}}_{o, p, q}} = \mathit{\boldsymbol{A}}_{o, p, q}^ * \cdot{\mathit{\boldsymbol{B}}_{o, p, q}} $

(4)

式中点乘的计算量正比于N³。

(4) 对C_{o, p, q}进行快速傅里叶逆变换(Inverse fast Fourier transform, IFFT)，从而得到相关函数c

$ {\mathit{\boldsymbol{c}}_{\alpha , \beta , {\rm{ }}\gamma }} = {\rm{IFFT}}({\mathit{\boldsymbol{C}}_{o, p, q}}) $

(5)

快速傅里叶逆变换的乘法计算量正比于N³ln(N³)，因此通过上述计算，将构象搜索原本o(N⁶)的总计算量降低到o(N³ln(N³))。然而，这里仅穷举和评价了配体分子b平移的所有情况，但还需要进一步考虑配体分子旋转的自由度。因此，可根据欧拉角旋转配体b，在每个角度下重复步骤(2)—(5)，搜索和评价配体的匹配程度，并最终挑出符合要求的近天然复合物构象。

FFT算法因其高效性而被广泛应用于快速穷举的计算几何匹配，静电相互作用和原子对偏好性。大多数蛋白质-配体分子对接程序都采用该方法，如3D-Dock ^[9]，ClusPro^[10]，ZDock^[11]，FTDock^[8]，DOT^[12]，GRAMM^[13]和MolFit^[14]等程序。

1.2 球面极坐标傅里叶变换方法

球面极坐标傅里叶变换(Spherical polar Fou rier, SPF)和快速傅里叶变换方法其本质一致，都是穷举采样构象空间。不同之处在于，球面极坐标采用正交球面极基函数的三维多项式展开来取代直接计算传统的三维笛卡尔坐标。故其具有的优势为：旋转配体操作可以通过变换其扩展系数来实现，极大降低了需要循环搜索旋转角度的计算量。甚至可以把旋转自由度也表示成系数，进行五维傅里叶变换。Hex^[15]和FMFT^[16]等程序主要采用球面极坐标傅里叶变换进行构象采样。

2 启发式搜索算法

启发式搜索算法是把对接体系中配体分子的平移和旋转操作先随机进行编码，并根据能量评分对操作后的配体构象进行优化和取舍，最终找到能量最低的配体分子构象。目前蛋白质-配体分子对接中常用的启发式搜索算法包括，蒙特卡洛(Monte Carlo，MC)算法^[17]、遗传算法(Genetic algorithm，GA)^[18]以及群智能(Swarm intelligence，SI)算法^[19]。

2.1 蒙特卡洛算法

MC算法有较好的物理学理论支持，因此在分子对接中应用也非常广泛。该方法首先在构象空间中随机地采样，并评估采样的函数值，通过函数值高低进行取舍，保留已得到的最优解作为最终解。MC算法不受研究问题搜索空间结构和分布的影响，在采样数趋近无穷时，理论上会收敛到全局最优解。该方法的物理基础与分子对接一致，即寻找能量更低的状态。由于不可能达到无穷次的采样，故采用重要性采样的方法，即Metropolis准则^[17]进行有偏好的采样。具体方法描述如下：

以配体分子的某个构象作为初始状态σ，通过随机平移和旋转产生新的状态σ′。计算两个状态的能量E(σ′)和E(σ)。如果新状态的能量更低，则有理由将其保留下来。反之，则需要进行判断，不能因为σ′状态的能量高就直接抛弃，这样会忽略了热运动的影响，并容易陷入局部极小。可通过计算Boltzmann因子的比值来判断，计算公式为

$ r = {\rm{exp}}\left( {\frac{{E\left( \sigma \right) - E\left( {\sigma '} \right)}}{{kT}}} \right) $

(6)

式中：T是绝对温度，k是Boltzmann常数。由于r在[0, 1]之间，且能量越接近，r值越接近1。故可以采用随机发生器产生[0, 1]之间的随机数ξ。通过比较r和ξ的大小，判断是否保留σ′。如果两个状态能量很接近，则r值接近1，此时有很大概率会保留σ′；反之，如果新状态σ′的能量过高，则被保留的概率将会很低。

MC算法在配体构象空间中随机采样，且不受构象空间结构和分布的影响，是一种通用的搜索方法。但该方法在实际使用时，可能需要较长的计算时间才能给出较好的解。RosettaDock程序^{[20, 21]}和GLIDE^{[22, 23]}程序就是采用MC算法进行全空间构象搜索，特别是RosettaDock程序，其生物大分子对接往往需要运行较长的时间。

2.2 遗传算法

GA算法将复杂优化问题视为生物的遗传和进化过程，先随机生成问题的解，然后进行交叉和突变操作，通过适应度选择使候选解不断优化，最终收敛于局部最优解甚至全局最优解。GA算法最早由Holland在1975年提出 ^[18]。该算法中，若干个体构成一个种群，个体采用二进制数据表示。先随机产生初始种群，即问题的初始解。然后通过交叉操作和突变操作来产生新的个体，通过适应度函数对新个体进行评价。根据轮盘赌的方式对新个体进行选择，适应度越高的则被选中的概率也越高，因此个体会逐步逼近适应度更高的情况，从而得到最优解。GA算法被广泛应用于分子对接构象搜索中，蛋白质分子一般固定不动，而配体相对于蛋白质分子的平移和旋转操作为个体的二进制数据，将蛋白质和配体相互作用的能量函数作为适应度函数。分子对接中，能量越低的构象对应的个体适应度越高，从而挑选出能量更低的复合物构象。

GA算法在分子对接搜索中的主要步骤是：(1)先产生初始种群，种群中的个体表示配体的平移和旋转操作；(2)平移和旋转配体构象，把对配体构象的操作，如平移和旋转进行二进制编码，对二进制数据进行交叉操作或者突变操作，即产生新的平移和旋转位置；(3)选择操作。用能量函数来评价配体分子的适应度，采用轮盘赌的方式对新个体进行选择，适应度好的个体有更大的概率被选择，最终保证种群的个数不变；(4)判断搜索是否已收敛，可以采用循环的步数或种群的能量函数得分情况判断收敛，如循环步数达到阈值则结束，或新种群与旧种群的能量函数分值差距在一定的阈值之内，也可认为收敛了。

GA算法的主要特点是不需要确定的规则，可以直接对个体对象进行操作，能自适应地调整搜索方向，不存在求导和函数连续性的限定，并具有较好的全局寻优能力。作为一种启发式寻优算法，GA算法以其简单通用、高效且易于并行处理等显著特点，在生物信息学各领域得到了广泛应用，已成为蛋白质-配体分子对接中最重要的构象搜索算法之一。GOLD^[24]和AutoDock^[25]等著名的分子对接程序就是采用GA算法进行构象搜索，并取得了较好的对接结果。

2.3 群智能算法

SI算法源于对自然界中的群体行为进行研究，包括蚁群、鱼群、鸟群、细菌群和其他微生物群等^[26]。SI系统由一群相互间交互和与周围环境交互的个体构成，并模拟群体自组织的行为。SI算法相对于GA算法来说搜索效率更高、全局寻优能力更强、收敛速度也更快，因而该方法很快被引入到分子对接的构象搜索应用中。

分子对接的构象搜索问题，SI算法中粒子群的应用较多，因此主要以粒子群为例介绍SI算法在分子对接构象搜索中的应用。在分子对接构象搜索中，目标搜索区域中每个粒子的位置代表问题的潜在可行解，即配体的旋转或平移操作，每个粒子根据周围粒子情况调整自己的位置和速度。能量函数作为适应度函数，评估每个粒子的当前位置。在整个过程中，每个粒子会根据其找到的最好位置，即个体最优解p_best，以及所有粒子的p_best中最好的位置，即种群最优个体g_best，来更新其速度和位置。粒子速度及位置更新公式为

$ \begin{array}{l} {\chi _i}\left( {t + 1} \right) = {\chi _i}\left( t \right) + {V_i}\left( {t + 1} \right)\\ {V_i}\left( {t + 1} \right) = \omega \cdot{V_i}\left( t \right) + {c_1}\cdot{r_1}\cdot({p_{{\rm{best}}}}\left( t \right) - {\chi _i}\left( t \right)) + {c_2}\cdot{r_2}\cdot({\rm{ }}{g_{{\rm{best}}}}\left( t \right) - {\chi _i}\left( t \right)) \end{array} $

(7)

式中：r₁，r₂为[0, 1]之间的随机数；c₁，c₂为常数。SwarmDock^[27]就是采用粒子群算法进行全空间构象搜索。本课题组对Autodock 3.0程序中全局搜索的GA算法采用蚁群算法进行了改写^[28]，提升了程序的搜索效率和寻优能力。

3 其他搜索方法

除了上述两大类构象搜索算法，还有其他一些课题组提出了不同于上述两类方法的搜索策略用于构象搜索。

3.1 基于图像识别的构象搜索算法

将蛋白质-配体分子识别看成三维图形的识别问题，一些课题组采用图像识别中的算法^[29]，如几何哈希方法。该算法先对配体分子进行预处理，对配体分子的表面关键点建立几何哈希表。然后进行分子识别，将蛋白质的表面特征与配体分子进行关键点匹配，评价配体不同方位的匹配程度，从而确定配体的构象。该方法主要针对分子表面特征的匹配，避免了分子对接中最耗时的平移和旋转操作，因而有效提升了搜索效率，采用该方法的代表程序主要有LZerD^[30]和PatchDock^[31]等。

3.2 基于分子动力学的构象搜索方法

分子动力学是模拟蛋白质和配体相互作用的有效计算方法，但由于其需要计算原子间的多种相互作用，并模拟动力学运动过程，因此极为耗时且容易陷入局部极小，在分子对接构象搜索中不太常用。HADDOCK^[32]程序将位点约束和分子动力学模拟相结合来进行分子对接。该方法利用预测的或实验提供的位点信息, 建立了模糊的相互作用约束(Ambiguous interaction restraints, AIRs), 并将其作为蛋白质与配体与关键位点相互作用的位置约束来加速分子动力学构象采样的效率，从而克服常规分子动力学的局限性。

3.3 基于片段组合的构象搜索方法

在考虑配体分子柔性构象产生的时候，一些课题组采用片段组合的方法进行构象搜索。该方法将配体分子看成多个刚性片段组成，刚性片段通过可旋转键连接。通过不同角度的组合来产生大量柔性配体构象，采用均方根偏差判断构象的多样性，并选取代表性构象保留。采用该方法的代表程序主要有Dock程序^[33]。

此外，在分子对接最终还需要进行构象的局部位置进行调整和优化。在局部构象调整和优化上，还有一些课题组采用了确定性优化搜索方法，如Autodock Vina采用拟牛顿法^[34]，ITS core-PP采用单纯型法^[35]等。

综上，蛋白质-配体分子对接的方法很多，表 1列出了目前应用较为广泛的蛋白质-配体分子对接软件, 并对其配体类型和构象搜索方法进行了简单介绍。

4 分子对接构象搜索的实验评测 4.1 国际评估实验

为了推动蛋白质分子对接技术的发展，2001年欧洲生物信息学研究所开始举办生物大分子相互作用预测评估实验(Critical assessment of predicted interaction，CAPRI)^[39-41]。截至目前，共进行了45轮实验。在该评估实验中，组委会首先选取尚未发表实验结构数据的蛋白质复合物作为预测内容，要求国际上参加实验的小组在规定时间内，从蛋白质单体结构出发，用理论方法对蛋白质复合物结构进行预测，并通过网络提交预测结果。随后，CAPRI评估小组会将这些结构与实验结构进行对比，来评估理论预测方法的准确性，并在CAPRI网上公布。大约每隔2~3年，召开一次CAPRI评估与研讨会。目前已成功举办了5届，它们分别是于2002，2004，2007，2009，2013和2016年在法国、意大利、加拿大、西班牙、荷兰和以色列举办的会议。

4.2 分子对接评估数据库

美国波士顿大学的研究小组建立的数据集是目前应用最广泛的分子对接评估数据集，已经发展到第5版了^[42]。数据集中有230个蛋白质复合物结构，并根据结合前后的构象变化分成151个刚性(Rigid-body)体系，45个中等难度(Medium difficult)体系，34个最难(Most difficult)体系。每个复合物结构有其相应的自由态单体(Unbound)结构。因此各小组可以测试自己的构象搜索方法，并比较构象变化对搜索算法的影响。

4.3 评估指标

配体均方根偏差(L_RMSD)是评价分子对接预测的复合物构象与实验解析的复合物结构相似程度的重要指标。L_RMSD为预测结构与实验结构中的受体叠落后，配体间主链原子(N，C，CA，O)位置的均方根偏差。CAP RI评估小组将预测结构分为高精度、中等精度、可接受和错误结果4种类型。高精度预测结果的要求是L_RMSD≤0.1 nm，中等精度预测结果是0.1 nm＜L_RMSD≤0.5 nm，可接受预测结果为0.5 nm＜L_RMSD≤1.0 nm，L_RMSD>1.0 nm的结构为错误结构。

5 结束语

蛋白质-配体分子对接一方面需要搜索能量低的构象，另一方面需要在尽量短的时间内搜索到各种可能的情况，因此快速有效的搜索算法是分子对接中重要的研究领域，迫切需要采用新的理论和计算方法提升现有程序的计算效率。为进一步提升构象采样的准确性和计算效率，以下几个方面可能将是未来的研究和发展方向：

(1) 蛋白质柔性采样方法。柔性依然是制约分子对接的一个技术瓶颈，蛋白质与配体相互作用的过程中会发生构象变化，考虑蛋白质主链的柔性会增加搜索的构象空间。可以考虑的方法包括精确确定蛋白质的柔性部位、缩小柔性搜索的范围，并发展快速有效的构象搜索算法。或者采用弹性网络以及图论等方法计算蛋白质的运动模式，作为新的自由度进行构象搜索和优化。另一方面，考虑蛋白质的柔性同时会增加构象数量，既有正确构象也引入了错误构象，因此还需要发展更为精确的能量函数，有效区分柔性对接产生的复合物构象。

(2) 基于模板的采样方法。对于蛋白质等生物大分子对接情况，直接采用快速穷举算法，尽管可以快速搜索整个构象空间，但同时也引入了较多错误的构象，为区分正确构象提供了难度。借鉴蛋白质折叠的思想，可以搜索蛋白质数据库中已有的结构模板，如果受体和配体的结构模板恰好有相应的复合物结构，则无需搜索构象空间，可直接根据搜索到的复合物结构进行对接构象模建或局部优化。即使没有完全一致的复合物结构，如果有类似的相互作用界面结构，也能为搜索提供有益的信息，减少构象空间的搜索范围。

(3) 面向新硬件平台上的构象采样方法。计算机处理器的发展日新月异, 特别是图形处理单元(GPU)的计算能力快速提高, 许多研究都根据新处理器改进原有算法。可以针对该处理器高效并行处理的特点, 发展基于GPU的构象采样加速方法来提升采样效率。另外，量子计算机等全新的计算技术进一步发展，将很可能颠覆复杂优化问题的计算框架和模型，因此也需要针对其特点建立相应的分子对接构象搜索方法。