---

今天要分享的是一篇蛋白质同源建模(MODELLER同源建模)与虚拟筛选相关的文献——《Performance of virtual screening against GPCR homology models: Impact of template selection and treatment of binding site plasticity》发表在PLOS COMPUTATIONAL BIOLOGY期刊上. 文献链接:[https://doi.org/10.1371/journal.pcbi.1007680]。 以下介绍的内容本人近期在做的工作，若有兴趣可以一起学习讨论。期待...

·Abstract

G蛋白偶联受体(GPCRs)属于细胞表面受体最大的家族，在许多生理过程中起着重要作用。 由于GPCRs是重要的治疗靶点，因此在计算机筛选中应用基于结构的方法来加速药物的发现过程具有重要的意义。然而，GPCRs难以结晶，80%的家族蛋白缺乏结构。作者以先前确定的晶体结构为模板，利用同源建模方法来评估GPCR结构的预测。并利用分子动力学模拟进一步细化结构。计算结果表明，最接近目标的同源物不一定是最好的模板，并证明了可以获得具有识别配体出色能力的精确结合位点模型。这一结果突出了结构预测方法的优缺点，为虚拟筛选成功应用于未知结构的蛋白质提供了指导。
作者评估了使用同源建模来预测两个治疗相关的GPCRs的结构，以及改进针对模型结合位点的虚拟筛选性能的策略。基于16种不同的GPCRs晶体结构，建立了D₂R和5-HT_2AR的同源模型。建模结果表明，可以得到准确的预测，但不一定使用最相似的模板。另外评估模型、模板的虚拟筛选性能是基于配体与诱饵（decoy）的分子对接。结果表明，必须对多个模板和基于每个模板的多个模型进行评估，以确定最佳的结合位点结构。基于胺能GPCRs的模型表现出大量配体富集，随着结合位点准确性的提高，虚拟筛选性能有提高的趋势。最好的模型得到的配体富集程度甚至比得上或优于D₂R和5-HT_2AR晶体结构。另外探讨了考虑结合位点可塑性的方法以进一步提高预测。结构集成分子对接并没有超过最佳的个体结合位点模型，但可以增加虚拟筛选命中的多样性，并有利于GPCR靶标与少数已知配体。分子动力学（MD）的改进导致了结构精度的适度提高，快照的虚拟筛选性能与原始同源模型相当或更差。这些结果为成功应用GPCR同源模型发现基于结构的配体提供了指导（黑体部分是对应的计算步骤方法）。

·Introduction

由于GPCRs控制着许多生理过程，刺激或阻止其活性的化合物是了解受体功能和具有治疗应用价值的工具。开发针对GPCRs的药物的努力取得了显著的成功。目前，美国食品和药物管理局批准的药物中有34%通过GPCRs发挥作用，用于治疗心血管、神经精神和神经退行性疾病等多种疾病。64个GPCRs的结构已经被解析，基于GPCRs结构的大型化学文库虚拟筛选，然后对最高评价化合物进行实验测试，已经成功识别了许多治疗靶点的先导物，包括生物胺、嘌呤、肽和蛋白结合受体。分子对接可用于>1亿化合物的化学文库搜索，并可识别命中率高达73%的配体，这表明虚拟筛选可为药物发现做出重要贡献。

尽管通过试验确定的GPCR结构越来越多，但>80%的受体的晶体仍缺乏原子分辨率信息。解决这个问题的一种方法是使用同源建模预测蛋白质结构。模板与靶标相似性大于30%可认为足够使用同源建模，但最终结构准确性取决于应用程序。尽管相似性>35%模板可以获得，但只有少数研究组确定了配体结合模式接近于试验确定的复合物。如果只有远同源的模板，没有研究组可以精确建模配体结合模式。

建模GPCR技术将使基于结构的配体发现扩展到未探索的药物靶点成为可能。评估模型常用的策略是使用分子对接来评估结合位点是否能识别诱饵中的已知活性化合物。高富集已知配体的模型被认为受体结构的良好代表，适合于虚拟筛选。虽然已有一些针对同源模型的虚拟屏幕已经取得了成功，然而还需要注意几个事项。首先，目前还不清楚该方法是否仅限于具有紧密相关模板的目标。药物复合物的预测将对侧链构象和形成结合位点的环区很敏感，很难用遥远的模板进行建模。通常的经验法则是选择序列一致性最高的模板来构建同源模型。而同源建模的基准并没有发现模板的序列身份与虚拟筛选性能之间有明显的相关性，有趣的是基于远源模板的GPCR模型有时会比紧密相关的模型获得更好的配体富集。其次，针对同源模型的前瞻性对接筛选通常针对结合位点的单一刚性结构进行。由于GPCRs是柔性蛋白，静态侧链的模型将无法捕获诱导拟合效应，并可能无法识别不同构象所识别的配体。考虑到多种受体构象可能会提高配体的富集，但也会增加虚拟筛选的计算成本。最后，同源模型会包含由目标和模板之间的结构差异引起的误差，随着序列一致性的降低，这些误差可能会增加。用更严格的方法如分子动力学(MD)模拟改进同源模型可以得到改进的结构。然而，原始模型和MD改进的同源模型的准确性和虚拟筛选性能很少被比较。

作者在这项工作中，以D2多巴胺受体(D₂R)和5-HT2A血清素受体(5-HT_2AR)为靶点。以16个晶体结构模板为基础，生成了代表近缘和远缘受体的同源模型。根据模型与最近确定的D₂R和5-HT_2AR晶体结构的一致性进行评估，并将已知配体与结合位点对接以评估其虚拟筛选性能。还研究了是否可以通过使用同源模型的集合或MD模拟改进预测。最后，评估了细胞外环区域对虚拟筛选性能的影响，以及模板的选择可能对某些配体化学类型造成偏向性。

·Result

1.基于16种不同模板的同源建模

根据12个胺能受体和4个非胺能GPCRs的晶体结构，建立了D₂R和5-HT_2AR的同源模型（如下图）。非胺能受体是一组亲缘性较远的蛋白质结构，用来覆盖识别不同配体的GPCRs。2种肾上腺素能(β₁AR和β₂AR)， 2种多巴胺(D₃R和D₄R)， 1种组胺(H₁R)， 3种血清素(5-HT_1BR， 5-HT_2BR ， 5-HT_2CR)， 4种毒菌素(M₁R，M₂R，M₃R，M₄R )，4个非胺能受体：视紫红质(Rho)、CXCR4趋化因子(CXCR4)、A2A腺苷(A_2AAR)和大麻素1 (CB1R)受体。

为了比较同源模型与D₂R和5H_2AR精度, 大多数选择非活性状态下确定的晶体结构和模板。少数情况下使用中间构象结晶的模板（例5-HT_1BR和5-HT_2BR）。GPCRs高度保守的拓扑结构促进了序列比对，并为第二个细胞外环(ECL2)手动调整，因为该区域参与配体识别。16个模板覆盖了广泛的跨膜螺旋™和结合位点(BS)与D₂R和5-HT_2AR的序列同源，5-HT_2CR和D₃R的晶体结构与5-HT_2AR(75%)和D₂R(77%)的TM序列一致性最高。正如预期的那样，四个非胺能模板的TM和BS序列与目标受体的识别率较低(分别为21-37%和6-19%)。（附：TM区域与BS区域有各自的定义，在前篇博客专门介绍过，具体可以参考以下两篇文献
《Integrated Methods for the Construction of Three-Dimensional Models and Computational Probing of Structure-Function Relations in G Protein-Coupled Receptors》、《What Can Crystal Structures of Aminergic Receptors Tell Us about Designing Subtype-Selective Ligands?》。这两篇也可以在本文介绍的文献中找到分别是Ref60，Ref52。）

对于每个模板，使用MODELLER建模250个模型，并挑出其中DOPE分数最好的50个作为下一步分析。由于比对中Gap和Insertion，由同一模板导出的同源模型在侧链构象和主链loop各不相同。对于基于相同模板的模型，其结合位点上的侧链的平均RMSDs在0.9到2.0Å之间。非胺能模板的结合位点结构变化最大，而BS序列一致性最高的模板的结合位点结构变化最小。不同模板得到的模型由于helix和loop区域的相对方向不同，主链结构也会有差异。

比较同源模型与晶体结构

D2R(PDB：6CM4)和5-H2AR(PDB：6A94)的晶体结构可在PDB数据库中下载。模型与native结构的TM主干、BS主干和BS侧链的平均RMSDs各记为RMSD_TMBB、RMSD_BSBB和RMSD_BSSC。虽然晶体结构的分辨率不足以区分电子密度图中的所有侧链原子，但RMSD_BSSC被计算来评估是否抓取了结合口袋的整体形状。

下图是将低(<30%)和高(>50%)序列以及30%-50%序列相似性模板进行比较。

对于两个靶标，RMSD_TMBB与TM序列相似性有适度的相关性(D₂R和5-HT_2AR的皮尔森相关系数各为R = -0.80和-0.59)。（附皮尔森相关系数是一种线性相关系数，是最常用的一种相关系数。用来反映两个变量X和Y的线性相关程度，R值介于-1到1之间，绝对值越大表明相关性越强。0.5-1.0之间表明有较强相关性。）

下图将基于不同模板的50种型号的D₂R (A-C)和5-HT_2AR (D-F)的平均RMSD_TMBB (A和D)、RMSD_BSBB (B和E)和RMSD_BSSC(C和F)的晶体结构进行比较。实线表示线性回归，R为皮尔逊相关系数。

5-HT_2AR的模板序列一致性更高，趋向于采用更准确的BS模型(RMSD_BSBB和RMSS_BSSC的R = -0.74和-0.77)。相比之下，随着序列一致性的增加，D₂R的结合位点精度得到了较小的提高(RMSD_BSBB和RMSS_BSSC的R = -0.29和-0.60)。相关性较弱的主要原因是，D₃R和D₄R模板(71%和94%的BS序列一致性)的RMSD值并不优于5-HT_2CR和肾上腺素能受体模板(50-60%的BS序列一致性)。血清素和肾上腺素能受体模板比多巴胺受体结构能更好地预测结合位点区域的TM6。

对同源模型和晶体结构分子对接和配体富集

对已知配体进行分子对接筛选，进一步评价同源模型。对于16个模板中的每一个，都准备了50个模型与DOCK3.6对接。共有822个和650个已知的D₂R和5-HT_2AR配体与性能匹配的诱饵结合到同源模型中，每个化合物在结合位点的数千个方向上取样，这些结合位点在计算中被固定。成功连接的化合物使用基于物理的评分函数根据它们的预测结合能进行排序。通过对> 8000万复合物结构的预测，评价了模型的虚拟筛选性能。而后根据ROC曲线分析诱饵上配体的富集情况，并使用曲线下面积的调整对数(aLogAUC)进行量化。aLogAUC更倾向于早期富集，这在虚拟筛选中是可取的，如果是正值，说明对接评分函数的表现要优于随机选择。例如，aLogAUC值为10对应于从随机选择中识别的配体数量超过预期的两倍。（附：logAUC与alogAUC只差一个随机选择的值，即alogAUC = logAUC - 14.5(随机选择的logAUC值)，logAUC也只有一条公式，比较简单，可以参看此文献中的Ref57。）

为了对富集结果进行分类，定义aLogAUC值< 10,10 - 15，> 15-20， > 20-25， >25分别为poor（差）, fair（中等）, good（良好）, very good（优好）, excellent（极好）。ROC曲线和aLogAUC值的形状变化较大，反映了预测的结合位点结构的差异。即使基于相似性最近模板的模型，alogAUC最佳和最差模型差异为14个单位，配体富集程度从差到优。模板Rho产生的模型富集度接近或者低于随机选择的预期值。对于所有其他模板，至少有一个模型的配体富集优于随机选择。根据aLogAUC的中值和最大值对同源模型进行分析。中值富集作为模板质量的衡量，aLogAUC值最高的模型被认为是最能代表受体结构的模型。12个胺能受体有7个产生的模型具有良好到优的富集。有5个胺能受体模板产生的模型富集结果是一般。其中有四种是基于毒蕈碱受体结构，与靶标共享较低的序列相似性。D4R模板也获得了差(D₂R)到中等(5-HT_2AR)中位配体浓缩。在四种非胺能模板模型中，虚拟筛选性能最差。模板(CXCR4和CB1R)和模板(CB1R)分别使D₂R和5-HT_2AR的中位富集良好。其他远源模板导致配体富集差或中等。

为了进一步评估模型的质量，我们分析了所选模板的最佳丰富结构。对于胺能模板，最好的富集模型显示至少良好的富集，12个模型中有8个具有很好的到极好的富集。在筛选的样品中，也有良好的到很好的结构具有最大的丰富度基于CXCR4-、CB1R-和A_2AAR的模型，表明这些模型在虚拟筛查中是有用的。

接下来是对结合模式的分析。结合模式好与差的定义是：如果配体与保守的Asp^3.32形成盐桥，并带有芳香部分向TM5/6延伸，则结合模式为良好。下图以D₃R、5-HT_2CR、CXCR4、A_2AAR和CB1R为基础，对50个*已知配体的预测结合模式进行了检测，以获得最好富集的模式。其中Aminergic包括D₃R，5-HT_2CR。

需要注意的是，这些结合模式特征不仅存在于D₂R和5-HT_2AR晶体结构中，在所有实验确定的胺能GPCRs结构中都是守恒的。对于基于最近胺能模板的同源模型，98-100%的排前配体对接到合理的结合模式。相比之下，在非胺能模板中，只有基于CXCR4的D2R模型具有较高的配体结合模式(96%)。对于其他基于非胺能模板的模型，预测的结合模式被判定为较差，只有0-15%的好位姿。其实，虽然有富集性很好的A_2AAR为基础的D₂R模型，但配体结合位点并没有得到正确的识别。绝大多数的*配体与Glu95^2.65形成盐桥，而不是与保守的Asp^3.32形成盐桥。同样，在基于CB1R的富集最佳的5-HT_2AR模型中，到Asp^3.32的盐桥被捕获，但配体向由TM2和TM7，而不是TM5。

对D₂R和5-HT_2AR的晶体结构也进行了对接筛选，D₂R筛选产生了一个优秀的配体富集(aLogAUC = 26.6)，该富集分数超过了所有800个评估同源模型的富集分数。5-HT_2AR晶体结构也获得了优异的虚拟筛选性能(aLogAUC = 26.1)，但富集程度略低于基于5-HT_2CR和5-HT_2BR模板同源模型的最大富集程度。接近100%的排前配体在晶体结构中具有合理的结合模式。

另外作者还用16个模板结构做分子对接计算富集，在16套D2和5-HT2AR模型中，分别有14套和16套同源模型得到优于模板的同源模型。一些模板结构很好的富集配体，这一结果可以解释为这些受体也识别生物按，因此结合许多D₂R和5-HT_2AR配体。对于远源模板，配体富集得到很大改进.

·Discussion

从GPCR家族的两个药物靶点建模和针对预测受体结构的广泛分子对接筛选中得出了四个主要发现。首先，我们的结果与选择具有更高序列同一性的模板将导致更好的模型的概念基本一致。然而，结构准确性和序列一致性之间的相关性并不是特别强，因为最好的结合位点模型不一定是基于最接近的同源物。其次，分子对接可以用于识别在虚拟筛选和配体富集中表现良好的模型，配体富集随着结合位点的准确性而增加。基于遥远模板(<30%序列同一性)的同源模型，由于结合位点的建模错误，大多数情况下不适合虚拟筛选，这使得很难获得受体-配体配合物的准确预测。第三，通过筛选结构整体来考虑结合位点的塑性一般没有改善虚拟筛选性能，但在有利的情况下，结果是可比的最佳个体模型。最后，在某些情况下，MD改进导致了结构精度的适度提高，但仿真快照的虚拟筛选性能与同源模型相当或更差。

·Materials and methods

1.序列比对和同源建模

12个胺能受体多序列比对使用带有默认参数的MAFFT localpair算法生成。对于这个MSA，使用hmmbuild从HMMER套件中获得一个HMM配置文件。如果在TM区域和细胞外环区域存在间隙，则手动调整所得到的比对结果，D₂R和5-HT_2ARTM区域的定义是基于MSA中排列良好的区域和胺能受体的结构信息。模板和靶点的ECL2中的保守半胱氨酸两个氨基酸对齐。在进行同源建模之前，我们将N-和C-端、细胞内3环以及ECL2在TM4和ECL2中保守半胱氨酸之间的延伸从比对序列中排除。使用MODELLER建立同源模型。每个模板-目标对共生成250个模型，提取DOPE得分最好的50个模型进行进一步分析。使用PyMol对同源模型和晶体结构进行了RMSD计算，其中考虑了侧链对称。

2.已知配体和诱饵的制备

已知的D₂R和5-HT_2AR配体(活性< 1uM，分子量< 350 Da)从ChEMBL20中提取。选择分子量范围是为了将基准集中在与用于虚拟筛选的化学文库(如ZINC数据库)中存在的性质相似的化合物上。不进行基于配体功能活性的过滤。D₂R和5-HT_2AR属性匹配的诱饵(配体分别为822个和650个配体，诱饵分别为55146个和43777个)采用DUD-E方法得到为了。评估特定配体化学类型的浓缩，D₂R配体与三个模板的共结晶配体相似使用聚类和子结构搜索确定(H₁R/Doxepin: 45 ligands, D₃R/Eticlopride: 38 ligands, and 5-HT_2CR/Ritanserin: 48 ligands)附件表8。化合物使用锌数据库协议通过DOCK3.6对接。

3.分子对接筛选

使用DOCK3.6程序进行分子对接计算，受体是用AMBER力场的一个版本描述。所有的Asp、Glu、Arg和Lys残基均采用电离侧链的参数。通过对氢键网络的目测，选择了His残基的互变异构体，D₂R和5-HT_2AR所有模型的互变异构体都是相同的。D₂R(Nε protonated: His33, His106, His393, and His398; Nδ protonated: His442)、5-HT_2AR(Nε protonated: His70, His165 and His182; Nδ protonated: His183).。同源模型的结合位点是根据模板的共结晶配体确定的。化合物通过DOCK3.6柔性配体算法与45个匹配的球对接到一个刚性受体结构上，这些球被标记为化学匹配。对接化合物生成的方向数由bin size、bin size重叠和距离公差决定，分别设置为0.4 A、0.1 A和1.5 A。每个被对接的化合物的结合能被计算为范德华作用能、静电作用能和配体失溶剂作用能项的总和。通过ROC曲线分析了配体在诱饵上的富集情况。为了量化配体的富集，对ROC曲线进行半对数变换，然后对曲线下面积积分，得到LogAUC值。调整后的logAUC (aLogAUC)是通过减去随机选择得到的logAUC值来计算。

以上并非文章全部内容，还有模型集成对接以及MD模拟改进模型两部分内容未涉及，但以上内容是此工作的第一步，后面内容是在此基础上进行的。工作内容分为两步：同源蛋白质建模+分子对接。两部分都有各的难点。在准备建模过程中，如果按照文章所说的，一步步处理的话会发现很多问题要处理，多序列比对和HMMER的目的是为了确定TM区域，那么我们可以简单过一遍操作（理解为主），主要是利用文章和参考文献给的信息确定TM1-7区域了。如果死磕步骤，我觉得挺浪费时间的。建模前要对比对文件手动更改，移除loop，加断链键等（如下图）。

因为模型过多，就不考虑结构优化之类的。最后还要对TM1-7区域抓出来计算RMSD（如文中所说）。建模就差不多到这完成了。

下面是分子对接，我觉得最为难处理的步骤是从数据库中提取数据，文章并没有详细提及如何获取配体以及诱饵数据，也没有提供数据，因为第一次接触虚拟筛选，完全不知如何处理数据。文章说的从chembl20提取分子量<350Da以及活性<1um的活性分子。但是chembl20是chembl的一个子库，chembl是一个生物活性类药物小分子数据库，包含二维结构、计算性质和生物活性。它每年都在更新生物活性分子数据，这样也算是改善了数据质量。但按所说的分子量以及活性过滤的分子比文章的多。而下载chembl20 MySQL库在我服务器上无法下一步操作，提取不了数据，这就让人很是头疼（若是有建议方法，请留言，感谢！）。对接软件的选择也没有硬性要求，但最好是自己熟悉的。最后对接完可以按照产生的配体和诱饵分数排名，在进行富集计算，画出ROC曲线计算alogAUC值（ROC曲线有很多的理解方式，找自己能够理解的处理起来也会方便点）。

最后的富集结果也不好，估计就是数据处理不当。
接下来看能不能抓紧处理完这一点工作，想再做些蛋白质结构预测（contact-map）、蛋白质设计之类的工作。什么都不懂，但是！什么都想学。抛瓦~