QuantumPDB：一种基于蛋白质结构的高通量量子簇模型生成工作流程

《Journal of Chemical Information and Modeling》：QuantumPDB: A Workflow for High-Throughput Quantum Cluster Model Generation from Protein Structures

　　高分辨率图像下载 MS PowerPoint 幻灯片
酶的计算建模提供了对催化作用的分子级见解，但从实验结构开始进行量子力学（QM）计算的准备工作是高通量研究的一个重要瓶颈。为加速这一过程而开发的自动化工具可能无法适用于不同的活性位点化学性质和几何结构。为了克服这些限制，我们推出了 QuantumPDB，这是一个 Python 包，它可以自动化地围绕活性中心生成层次化的坐标/相互作用球体，从而直接从原始蛋白质结构创建 QM 簇模型。该工作流程集成了结构清洗、质子化状态分配和 QM 计算设置。它使用基于 Voronoi 划分的接触球体构建的具有化学意义的模型，能够准确表示复杂的活性位点几何结构。我们概述了我们的模块化代码，并描述了如何使用它来自动化高通量蛋白质筛选。为了展示其实用性，我们整理了一个包含 989 个全酶的数据集，并对这些酶中的 842 个进行了 QM 计算。计算出的性质分析表明，使用密度泛函理论模拟的酶环境能够一致地将底物电荷调节到中性，并减小底物的偶极矩。即使在活性位点主要由中性残基组成的情况下，这种现象也是普遍存在的。通过自动化和标准化多球体 QM 模型的构建，QuantumPDB 为大规模、数据驱动的蛋白质研究提供了一个强大的平台。

1. 引言
酶在生理条件下以高效率和选择性催化复杂的化学转化。(1?3) 它们的显著活性源于电子和静电性质，如极化、(4,5) 电荷转移、(6?8) 局部电场、(9?16) 以及精确的底物定位(17?19) 和构象动态(20?23)。为了模拟这些多样的性质，需要互补的计算策略。(24) 经典分子动力学（MD）用于模拟蛋白质动态(25)，但它无法描述电子变化，(26,27) 这些可以通过从头算方法（如量子力学（QM）簇模型(28,29)或多尺度量子力学/分子力学（QM/MM）方法(30,31)来提供。然而，实验晶体结构的准备工作存在重大挑战。蛋白质的晶体结构通常包含未解决的区域、晶体学伪影(32,33)、非蛋白质成分（如辅因子、配体、核酸、碳水化合物、离子和水分子)(30,34)，并且通常缺乏氢原子(35,36)。金属酶引入了额外的挑战，需要 QM 处理来准确描述活性位点的电子环境，这对其氧化状态、自旋状态和金属的精确 coordination 几何形状非常敏感。(37) 虽然单个系统研究存在许多困难，但这些准备工作挑战在试图从大量多样化的酶数据集中发现通用催化原理的高通量量子力学研究中被放大了。(38?41) 解决酶的高通量建模挑战需要一个强大的自动化工作流程，该流程标准化了从结构准备和 QM 簇模型构建到运行 QM 计算以及将计算出的 QM 性质提取到一致的数据集中的过程（图 1）。有几个软件包旨在促进高通量酶计算的结构准备和模型构建，例如 Protein3D(42,43) 用于 QM，QMzyme(44) 和 RINRUS(45?49) 用于 QM 和 QM/MM，以及 EnzyHTP，它是 Mutexa 生态系统的核心组件(50?53)，适用于 MD、QM 和 QM/MM 的所有层面。(39,54,55) molSimplify 中的 Protein3D 类(42) 自动化提取金属酶的主要坐标球体，但不提取第一层坐标球体之外的残基，而这些残基往往对反应性至关重要。(42,56) QMzyme 通过处理预平衡的分子动力学轨迹来生成 QM 和 QM/MM 输入文件，但不进行结构预处理。(44) RINRUS 根据残基相互作用网络选择 QM 区域，提供化学直观的 QM 区域，但缺乏自动化结构准备功能，例如建模缺失的环。(45) EnzyHTP 从 PDB 结构自动生成突变体，用于 MM、QM 和 QM/MM 模拟，但它依赖于基于距离的截止标准，可能会包含无关的残基或忽略 QM 区域中的关键相互作用。(39)

图 1
图 1. 酶的高通量 QM 研究的自动化工作流程步骤：1) 结构准备，2) 生成 QM 就绪的结构模型，3) 执行 QM 计算，4) 提取计算出的 QM 性质，5) 编译 QM 性质的数据集。

高分辨率图像下载 MS PowerPoint 幻灯片
尽管做出了这些努力，但从原始 PDB 条目到包含主要球体之外残基的 QM 就绪簇模型的完全自动化工具尚未出现。QMzyme 和 EnzyHTP 等工具中使用的基于距离的截止标准可能会导致无关残基的错误包含或关键相互作用残基的遗漏。这在非球形活性位点中尤为明显，例如乙酰胆碱酯酶的“峡谷”(57,58) 和碳酸酐酶的“锥体”(59,60)，其中球形 QM 簇可能会忽略理解特定蛋白质活性位点相互作用所必需的残基。

为了解决这些挑战，我们开发了 QuantumPDB，这是一个完全集成且可扩展的 Python 包，旨在直接从任何 PDB 格式的结构文件（包括直接来自 PDB 的文件）简化 QM 簇模型的生成。QuantumPDB 系统地纠正了起始结构中的结构问题，包括添加缺失的氢原子、建模未解决的环以及分配质子化状态。它使用通过 Voronoi 划分识别的基于接触的相互作用球体(61?63) 将残基分类为可解释和具有化学意义的主要、次要和三级接触球体，同时允许用户将活性位点修剪到目标系统大小。Voronoi 划分是一种生成由彼此接近的点组成的 Voronoi 单元的技术，在生物化学的其他应用中得到了广泛利用。(64?66) 使用 Voronoi 划分使我们能够识别蛋白质图并将其划分为紧密接触的区域。这些直接相互作用球体以及同心次级和三级球体的分类也有助于 QM/MM 或 QM/QM′（例如 ONIOM(67)）的计算，包括使用新兴的机器学习势能。(68) QuantumPDB 相对于其他替代方案的主要优势在于它是一个完全开源的代码，可以与多种 GPU 和 CPU 量子化学代码一起使用，并且可以扩展到其他代码。对于那些不熟悉自己蛋白质准备工作流程的新手用户来说，这也降低了入门门槛。此外，QuantumPDB 中的这些新功能也可以被其他支持这种模块化的代码（如 EnzyHTP）所利用。为了展示 QuantumPDB 的实用性，我们分析了不同酶系统中的底物-活性位点电荷转移和底物偶极矩，揭示了酶和底物之间电荷转移的一般趋势，突显了电荷转移的重要性。QuantumPDB 显著减少了劳动密集型的准备工作，为加速高通量酶分析工作流程提供了新的途径。

2. 代码结构概述
2.1 QuantumPDB 的设计架构
QuantumPDB 包是用 Python 以模块化方式实现的。该包主要通过命令行界面（CLI）与用户交互，该接口操作一个单一的 YAML 格式的配置文件。虽然底层功能可以独立调用，但典型的 QuantumPDB 工作流程分为三个不同的阶段，每个阶段由一个单独的命令启动：qp run、qp submit 和 qp analyze。qp run 命令执行初始的结构准备和 QM 簇模型的生成。随后，qp submit 处理 QM 计算的创建和提交，而 qp analyze 执行已完成 QM 计算的后处理分析。这种分离允许使用不同的参数重新运行分析，而无需重复初始设置和计算阶段。YAML 文件指定了工作流程的所有用户定义参数，包括结构准备、QM 簇生成、QM 计算设置和输出分析的选项（支持信息图 S1 和表 S1）。YAML 文件还包含布尔参数，允许用户启用或禁用管道的主要阶段——结构准备、簇生成和 QM 作业管理。用户提供结构输入，可以是用于自动下载的 PDB 访问代码或本地 PDB 文件。由于这种灵活性，该代码可以处理来自 X 射线晶体学、冷冻电镜 (cryo-EM) 和核磁共振 (NMR) 的结构，以及来自 MD 模拟和生成模型的计算结构。

为了提供软件框架的高层概述，该包被组织为五个顺序子包：qp_structure、qp.protonate、qp.cluster、qp.manager 和 qpanalysis，它们对应于 QM 准备流程的不同阶段（图 2）。在第一个 qp_structure 子包中，代码从本地或 PDB Web 服务器获取结构文件，并使用 Modeller (32) 软件包执行初始的结构修复，包括重建缺失的原子、残基和环。qp.protonate 子包然后使用 Protoss 服务器添加氢原子并分配质子化状态。(69?71) 如果输入的结构文件已经包含了完整的原子结构和其质子化状态（例如，来自全原子分子动力学模拟的快照），工作流程的模块化设计允许用户跳过准备步骤并立即开始簇模型的构建。qp.cluster 子包然后使用基于接触的相互作用球体围绕用户定义的核心残基或原子构建 QM 簇模型，从而灵活地构建残基球体。模型生成后，qp.manager 自动化创建 QM 输入文件、提交 QM 计算，并使用 TeraChem 和 ORCA 进行作业执行监控。最后，qp/fixtures 提供工具从 QM 输出中提取和计算性质，如部分电荷和偶极矩。每个子包的输出作为后续阶段的直接输入，并将输出写入标准目录层次结构中，确保从初始结构到最终结果的数据透明且可复制。

图 2
图 2. QuantumPDB 包的层次化工作流程。五个顺序模块及其主要功能。(1) qp_structure：获取 PDB 文件并修复缺失的原子和残基。(2) qp.protonate：分配质子化状态并评估原子占据情况。(3) qp.cluster：使用 Voronoi 划分生成相互作用球体。(4) qp.manager：创建 QM 输入文件并提交计算。(5) qp.analysis：对 QM 输出进行部分电荷和偶极矩分析。

高分辨率图像下载 MS PowerPoint 幻灯片
2.2 QuantumPDB 的子包
2.2.1 qp.structure
QuantumPDB 工作流程由 setup 模块启动，该模块从 YAML 配置文件中解析用户规格。通过执行 qp run -c config.yaml 命令后，get_pdbs 函数确定蛋白质结构的来源。用户可以提供 PDB 访问代码，然后使用 fetch_pdbs 来检索，或者用户可以指定本地 PDB 文件的路径。当提供 PDB 代码时，程序会验证是否可以通过标准 API 从蛋白质数据银行下载相关的 pdb 文件，如果失败则会显示错误消息提醒用户。对于高通量应用，可以提供 PDB 代码列表或 CSV 文件的路径，qp.structure 从 csv 的 pdb_id 列中读取 PDB 代码。PDB 通常包含同一蛋白质的多个副本，因此为了准确构建 QM 簇，建议用户应用 PDB 查询过滤器以确保高分辨率结构具有明确的活性位点（例如，通过 EDIA 分数判断(72)）。对于金属酶，由于金属的氧化态和自旋态通常会有所不同，用户还必须提供同一 CSV 文件中氧化态和自旋多重性的信息。所有电子状态在执行前都在 CSV 中指定，并在工作流程执行期间自动读取，从而无需额外用户输入即可批量处理多个结构。用户提供的这些信息是必不可少的，因为仅从结构数据中无法明确确定这些电子性质。missing module 然后使用 get_residues 函数通过解析 PDB 文件头中的 REMARK 465 和 REMARK 470 条记录来识别缺失的残基和重原子。clean_termini 函数移除 N-末端和 C-末端的未解决残基。然后使用 Modeller (32) 软件包构建结构的任何缺失部分，包括大的链内断裂。write_alignment 生成对齐文件，然后将其用作 Modeller 构建缺失的环、残基和重原子的输入。用户通过使用`optimize_select_residues`参数来控制添加残基后结构放松的程度，该参数指定是仅优化新建模的区域还是整个结构。如果同一结构文件中存在多个模型（例如NMR结构集合），Modeller库还会选择第一个模型。如果PDB文件中的原子条目中的`altLoc`字段指定了一个替代构象，Modeller会选择该替代构象。该模块的功能仅限于与Modeller库兼容的原子和残基，无法对底物中的缺失原子或非标准辅因子进行建模。如果Modeller无法建模缺失的原子，系统会发出警告。未来版本计划支持其他结构构建工具（例如AlphaFold）。虽然这段代码主要是为处理实验性的PDB文件而设计的，但用户可以选择预先修改文件，例如提供由分子动力学生成的文件或修改了活性位点的文件以包含催化活性中间体。类似的策略也可以用于生成突变体，可以通过将新残基与原始残基对齐或通过分子动力学放松来实现。QuantumPDB的未来发展将侧重于增加对催化中间体、突变体和分子动力学轨迹的支持。`qp.structure`子包的最终输出是一个完整的PDB文件（例如`pdb_modeller.pdb`），作为后续`qp.protonate`阶段的输入（图3）。

**图3**
**图3.** `qp.structure`和`qp.protonate`子包的架构概览，分别用绿色和蓝色表示。功能以橙色框显示，并标有相应的模块。输入和输出的结构文件用黑色圆圈表示，其他非结构输入和输出文件用黑色框表示。

**高分辨率图片**
**下载MS PowerPoint幻灯片**

**2.2.2. qp.protonate**
在使用Modeller进行结构精细化处理后，`qp.protonate`子包执行质子化状态分配和进一步的结构精细化。尽管Modeller库可以处理大多数重复结构的情况，但在不同的NMR构象模型中没有分离或用`altLoc`记录标记的一些模糊坐标仍可能在清理后的结构中引起冲突。因此，`clean_occupancy`函数通过选择一组单一的、自洽的坐标来解决剩余的替代构象和部分占位的原子（参见支持信息文本S1）。占位小于1的原子表示部分占用，并对应于结构替代方案，但在进一步的结构准备（即质子化）之前，我们需要选择一个单一的姿态。QuantumPDB基于一系列层次规则选择最可能的姿态，如果所有其他因素相同，则优先选择占据度更高和原子更明确的残基（参见支持信息文本S1）。该函数将所有具有替代构象的残基放入队列中并按顺序处理它们。对于每个残基，它评估与其邻居的重叠情况，仅保留优先级得分更高的构象。此层次优先级方案保留了用户定义的中心活性位点残基，其次是典型的氨基酸和其他残基类型（例如MSE、HYP）。对于同一优先级类别中的残基，保留平均占据度更高的构象（参见支持信息文本S1）。然后，代码将精炼后的结构提交给Protoss服务器（69-71），该服务器通过优化氢键网络和枚举异构体来添加氢原子并分配质子化状态（73-75）。这种质子添加适用于蛋白质残基和任何水分子。如果水分子或质子是通过分子动力学添加的，可以通过在输入YAML中设置`protoss: False`来跳过Protoss步骤。在添加氢原子的过程中，Protoss还会识别产生的空间冲突并自动移除冲突的残基。为了确保最终结构的物理真实性，QuantumPDB实现了自动化反馈循环来纠正这些空间冲突。`parse_output`模块解析Protoss日志文件以识别由于冲突而被删除的残基。如果检测到冲突，会导致冲突的残基会被报告给用户，工作流程会返回到缺失的模块，其中`delete_clashes`函数会在Modeller重新构建该区域以供再次提交之前删除这些问题残基（图3）。用户通过`max_clash_refinement_iter`参数控制精细化循环的次数。

一旦获得无冲突的结构，`adjust_activesites`函数应用一组基于启发式的结构校正，以提高金属酶活性位点的结构保真度。组氨酸咪唑环会被重新定向，以纠正潜在的晶体学错误赋值，即配位氮原子指向远离金属中心的情况。对于Protoss未识别的非标准残基（例如亚硝化水合酶结构中的S-羟基半胱氨酸（CSO），也会添加氢原子。最后，其他配位金属的残基（例如Cys、Tyr）以及辅因子和配体会被去质子化，以达到它们在配位金属时的正确电荷状态（参见支持信息文本S2）。`ligand_prop`模块确定配体、辅因子和简单离子的净电荷和自旋多重性。`compute_charge`函数从Protoss生成的SDF文件中读取形式电荷，根据结构变化进行校正，并为在生物条件下自旋和氧化态不变的阳离子添加电荷（参见支持信息文本S3）。`compute_spin`函数为预定义的自由基种类（如NO和O2）分配自旋多重性（参见支持信息文本S3）。具有可变自旋和氧化态的金属不通过`ligand_prop`模块处理，必须由用户在输入YAML中的CSV文件中提供。

**2.2.3. qp.cluster**
`qp.cluster`子包通过围绕用户定义的中心划分蛋白质环境来构建QM簇模型。这是对蛋白质的每个链分别进行的，因此如果PDB中存在多个链，则会产生多个QM簇。由于不同链的结构可能略有不同，用户需要评估每个链的相对质量和合理性。这种划分首先通过定位簇中心、使用Voronoi镶嵌建立原子接触网络并在用户定义的中心周围生成连续的相互作用球体来实现（图4和支持信息图S2）。用户通过`center_residues`参数或输入CSV中的`center`列来指定活性位点中心。有三种主要的方法可以越来越具体地选择中心。对于需要链和残基ID的高度具体选择，用户可以使用`[残基名称]_[链ID][残基编号`的表示法来指定任何类型的残基，例如链A中位置200的琥珀酸残基将被表示为SIN_A200。对于多个残基的目标选择，用户可以使用`center_residues`参数指定多个残基。例如，可以使用`[FE_SIN`语法将任何铁中心与其琥珀酸（SIN）辅因子合并，其中`merge_distance_cutoff`值用于指定合并距离。此外，还可以使用连字符将这些残基合并到一个中心（例如：His_A198-FE_A199-SIN_A200）。使用连字符而不是逗号，以便表示法也与输入CSV文件兼容。

**图4**
**图4.** 使用`qp.cluster`子包为TauD（PDB ID：1OS7）生成的基于接触的簇（76），其中用灰色表示初级球体，蓝色表示次级球体，紫色表示三级球体。原子的颜色如下：碳为灰色；氮为蓝色；氧为红色；硫为黄色；氢为白色；铁为深橙色。配位键用紫色虚线表示。

对于最通用的选择，用户可以将`center_residues`参数定义为一个单一值，例如[Fe]，这样输入结构中的每个Fe原子都将成为其自己的独立簇模型的中心。当不知道所需中心的确切位置和数量时，这种方法非常有用。当结构具有多个功能单元时，代码会找到所有匹配残基名称的中心并为所有中心创建簇。为了区分多核中心和来自不同功能单元的分离中心，用户可以指定一个`merge_distance_cutoff`值，将此距离内的中心残基合并为一个中心。然而，这种方法仅限于具有HETATM记录的残基，包括当它们作为底物时的氨基酸，因为指定蛋白质中的氨基酸将为该氨基酸的每个实例生成一个簇。

一旦用户选择了中心，`get_center_residues`函数会定位所有这些残基的实例。为了处理多金属活性位点或具有多个配体或辅因子的活性位点，可以提供一个`merge_cutoff`距离，将指定的中心残基合并为一个实体。然后`voronoi`函数使用SciPy库从所有原子坐标计算Voronoi镶嵌，以建立基于接触的网络。从相邻Voronoi单元之间的棱线构建原子级邻接列表，识别所有直接接触的原子。在稀疏区域（如配体结合口袋、蛋白质-蛋白质界面或溶剂暴露的表面），由于某些方向的邻近原子缺失，Voronoi单元的几何形状变得高度各向异性和拉长。这导致后续模型构建中的簇边界不规则。为了在低密度区域（如蛋白质表面和配体结合口袋）规范单元边界，`fill_dummy`函数在蛋白质周围的3D网格上放置虚拟原子，从而提高镶嵌的分辨率。这防止了残留物的不规则划分，并形成了密集、各向同性和几何形状规则的Voronoi单元（参见支持信息文本S4）。

在定义了原子接触网络之后，通过构建和细化相互作用球体来生成最终的QM簇模型。该过程包括迭代构建球体、根据用户指定的原子计数限制进行可选的剪枝、在簇边界处封顶断裂的键以及写入最终坐标文件（图4）。如果原始PDB文件中的任何水分子（例如晶体水）位于这些相互作用球体内，它们将会被包含在内，但代码不会生成额外的水分子来填补结构中的空缺。`get_next_neighbors`函数通过将残基基于与前一个球体的直接原子接触分配到连续的、不重叠的层中来迭代构建相互作用球体，其中接触被几何定义为Voronoi单元之间的共用边界。簇构建算法适用于整个结构级别，跨越多个肽链，这使得来自不同亚基的残基能够贡献于多聚酶寡聚体接口的簇。第一个球体由`radius_of_first_sphere`参数定义，默认值为4.0 ?的半径，以捕捉从短程强氢键相互作用到长程疏水相互作用的多种相互作用类型（参见支持信息文本S5）。与第一个球体不同，第二个及后续的球体不是直接由距离截止值定义的，而是基于Voronoi派生的接触网络构建的，其中每一新层由网络拓扑邻域内的残基组成。如果用户指定了`max_atom_count`，`prune_atoms`函数会系统地从整个簇中移除最远的残基，直到原子计数低于阈值，但不包括封顶原子。`cap_chains`函数解决簇边界处的开放价态。断裂的C–N肽键会根据`capping_method`关键字指定用氢原子（build_hydrogen）或N-甲基乙酰胺（NME）和乙酰（ACE）基团（build_heavy）进行封顶。对于其他断裂的非肽键，则用氢原子封顶（参见支持信息文本S6）。在未来的工作中，将重点研究仅保留侧链的替代方法，这在蛋白质主链上有许多切割的非常小的QM簇模型中可能更合适。最终簇模型及其组成球体由`struct_to_file`模块写入PDB和XYZ坐标文件。`compute_charge`函数计算簇的总电荷，`charge_check`函数执行一致性检查，以确保电子数、总电荷和自旋多重性兼容（参见支持信息文本S7）。每个球体的总电荷和残基组成保存在`charge.csv`和`count.csv`中，用于后续的QM计算。

**2.2.4. qp.manager**
在生成QM簇模型之后，`qp.manager`子包创建并管理QM计算的作业提交。用户通过执行 `qp submit -c config.yaml` 命令来启动此状态。该命令会为量子化学软件创建输入文件，处理高级计算的设置（如MM点电荷嵌入），并管理作业提交到高性能计算（HPC）调度器的过程（见图5）。`create` 模块中的 `create_jobs` 函数会根据用户提供的输入数据生成所需的文件，以适用于TeraChem（用于GPU加速）或ORCA（用于CPU计算）的QM计算。QM簇的总电荷和自旋多重性是通过将用户通过 `get_electronic` 函数提供的CSV数据与 `Protoss` 函数返回的净电荷和自旋值结合来确定的。利用 `job_scripts` 中的模板，该模块生成包含用户定义参数（如方法、基组、介电常数）的QM软件输入文件，以及支持SLURM和SGE调度器的提交脚本。这部分代码的目的是实现高通量计算的自动化（见第4节）。如果用户prefer使用自己的接口进行作业调度，可以通过关闭QM作业的创建或提交功能来绕过这部分代码（支持信息表S1）。

**图5**
图5. `qp.cluster` 和 `qp.manager` 子包的架构概览，分别用紫色和灰色表示。函数以橙色框显示，并标明了相应的模块。输入和输出的结构文件用黑色圆圈表示，其他非结构性的输入和输出文件用黑色框表示。

为了更准确地描述活性位点，软件会自动生成代表整个蛋白质环境静电影响的MM点电荷。如果在配置文件中启用了电荷嵌入功能，`charge_embedding` 模块会准备必要的点电荷文件。该模块首先使用 `rename_and_clean_resnames` 函数标准化组氨酸质子化状态和末端残基的名称，然后为整个蛋白质结构中的每个原子分配部分电荷。这种参数化不仅包括蛋白质残基，还包括任何水分子（使用标准TIP3P点电荷处理），以及结构中存在的常见单价和双价离子（如Na+、K+、Mg2+、Cl–）。所有电荷值都来源于内部字典 `ff14SB_dict`（与ff14SB力场一致），或者来自用户通过 `charge_embedding_charges` 配置选项指定的JSON文件中的部分电荷数据，从而支持使用其他力场。如果用户没有指定，将排除电荷字典中不存在的残基（如非标准氨基酸、糖类或辅因子），并会发出警告列表中遗漏的残基名称。为了生成点电荷文件，`remove_qm_atoms` 函数会识别属于QM簇的原子并将其从完整蛋白质结构中移除。最后，使用 `charge_embedding_cutoff` 参数指定在QM簇质心用户定义距离内的所有剩余MM原子会被写入 `ptchrges.xyz` 文件。选择是基于残基的；如果某个残基的任何原子都在这个距离范围内，则该残基的所有原子都会被包括在内。默认的截止距离为20.0 ?。在QM区域和MM区域边界处的断裂肽键会被氢原子或ACE/NME基团替换，从而在MM原子原本所在的位置引入QM原子。使用 `remove_qmatoms` 函数，距离任何QM原子0.5 ?以内的MM原子也会被从点电荷场中移除，以防止重复计算。

`submit` 模块负责作业的自动提交和监控。提交到队列中的作业数量由 `job_count` 参数控制。`manage_jobs` 函数通过检查是否存在隐藏的 `submit_record` 文件（即以 periods 开头的文件）来识别尚未提交的计算目录。该文件由 `create-marker` 函数在作业提交时创建，其中包含作业详细信息，防止后续运行中重复提交。

**2.2.5 qp.analysis**
QuantumPDB工作流的最后阶段是对QM计算结果的后处理和分析。`qp.analysis` 子包通过 `qp analyze -c config.yaml` 命令执行，提供用于诊断已提交作业状态和对完成的QM计算输出进行自动化电子性质分析的工具（支持信息图S3）。`checkup` 模块会自动评估通过 `qp.manager` 子包提交的所有作业的状态，该子包支持SLURM和SGE调度器。代码会遍历每个计算任务，解析QM输出文件和隐藏的 `submit_record` 文件，并将作业分类为已完成、正在运行、排队、积压或处于特定错误状态（如电荷或内存相关错误）。对于成功完成的作业，`multiwfn` 模块作为Multiwfn软件包的自动化包装器，用于从生成的波函数计算电子性质（如部分电荷方案和偶极矩）。然后，代码会找到所有包含 `.molden` 文件的已完成作业，并执行用户指定的分析任务，例如使用 `charge_scheme` 函数根据用户通过 `charge_scheme` YAML参数指定的一个或多个方案（如Hirshfeld、Mulliken、CM5）计算原子部分电荷。此外，`calc_dipole` 函数还会计算指定残基的分子偶极矩，将残基的质心置于原点。

**3. 适用蛋白质类型的回顾**
为了证明QuantumPDB适用于各种蛋白质结构，我们通过一些具有代表性和挑战性的案例来展示其实用性。虽然QuantumPDB尚未针对所有可能的用例进行测试，但我们认为通过关注这些应用，QuantumPDB在最常见的情况下应该表现良好。具体来说，我们展示了其在金属蛋白协调环境以及结合寡聚体和多残基配体或底物的蛋白质中的应用。

对于单核金属酶，用户应通过配置文件中的 `center_residues` 参数指定中心金属作为簇生成的起点，例如在非血红素铁酶TauD中的 [FE]（图4）。然而，许多金属酶的活性位点包含多个金属离子，这些离子通常彼此靠近并通过桥接配体连接。对于这些多金属活性位点，例如甲烷单加氧酶（MMO）中的双铁中心，将双金属核心视为一个单一的功能单元是必要的，以便构建一个包含协调残基和围绕双铁核心的第二协调层的化学意义上的QM簇模型。在MMO结构（PDB ID: 1FYZ）中，使用 `[FE2]` 作为 `center_residues` 并设置 `merge_distance_cutoff` 为4.0 ?，可以捕获PDB中彼此距离在4.0 ?以内的所有Fe(II)离子，并为每个合并的中心生成一个簇模型（图6）。或者，可以使用 `[FE2_A5001-FE2_A5002` 指定二聚体复合物中的单个铁中心，将链A中的Fe(II)残基5001和5002合并为一个中心。这种方法确保了第一个相互作用层围绕整个双铁核心构建，正确捕获了所有桥接和末端配体以及第一层内的协调残基。需要注意的是，合并金属中心的能力不仅限于距离相近且共享协调层的离子。例如，肽基甘氨酸α-羟化酶（PHM；PDB ID: 1PHM）中的两个铜位点相距超过11 ?，但可以通过指定两个铜离子（例如 [CU_A357-CU_A358] 并将 `merge_distance_cutoff` 设置为大于它们之间的距离（例如12.0 ?）来将其视为一个功能中心。这种统一的中心可以生成包含两个金属位点的第一个相互作用层（图6）。

**图6**
图6. QuantumPDB生成的代表性金属酶的QM簇模型。（左上）甲烷单加氧酶（MMO；PDB ID: 1FYZ）通过合并两个铁原子定义双铁中心。（右上）肽基甘氨酸α-羟化酶（PHM；PDB ID: 1PHM）通过合并两个铜原子定义长距离双铜中心。（左下）氧合肌红蛋白（PDB ID: 1MBO）将整个血红素辅因子（铁和卟啉）和结合的O2分子定义为中心。（右下）腈水解酶（NHase；PDB ID: 3A8O）定义铁中心，突出显示协调的骨架酰胺和非标准CSO和CSD残基。MMO、PHM、氧合肌红蛋白和NHase的初级和次级层原子数分别为102和613、46和421、453和1112、105和611。定义为中心的原子用黑色轮廓表示。初级、次级和三级层分别用灰色、浅蓝色和紫色显示。当金属协调残基位于第一层时，第一层用棒状图表示（MMO、PHM和NHase）；否则用表面表示（MbO2）。原子的颜色如下：氮蓝色；氧红色；硫黄色；氢白色；铁深橙色；铜芥末黄色；碳其他位置为灰色。协调键用紫色虚线表示。

对于不仅包含金属还包括协调配体和底物的活性位点，也可以直接定义簇中心。以血红素蛋白为例，活性位点由卟啉环、金属中心、配体（如O2）和远端协调残基组成。用户可以使用 `[FE_A155-HEM_A155-OXY_A555-His_A93` 表示法将金属、卟啉（HEM）、结合的O2和协调的组氨酸合并为一个中心（如氧合肌红蛋白的结构（PDB ID: 1MBO））。通过将铁-卟啉-O2和远端残基视为一个整体，在整个活性位点周围构建第一个相互作用层，从而保留了可能影响电子性质的重要非共价相互作用（图6）。由于许多含血红素的蛋白质还具有远端血红素协调残基，这些额外残基也可以包含在中心中。如果它们没有被包括在内，它们会自动被纳入第一个相互作用层中。由于代码至少生成两个相互作用层，这种区分并不重要。腈水解酶（NHases）的特点是含有基于半胱氨酸的硫代酸（CSO）和磺酸（CSO）酸以及脱质子化的骨架酰胺，这些情况下自动校正质子化状态是必要的（支持信息文本S2）。

QuantumPDB还通过其灵活的中心定义解决了围绕寡聚体（即多个共价连接的残基）底物构建模型的挑战。用户可以指定寡聚体的一个或多个亚单元，`get_center_residues` 函数会根据 `merge_distance_cutoff` 参数将它们合并为一个整体。这个统一的中心随后用于生成基于接触的相互作用层，适用于结合多残基配体或底物的蛋白质（图7）。例如，C型凝集素Langerin（CD207；PDB ID: 3P5F）的结合位点包含Ca2+离子、甘露糖寡糖，并且暴露在溶液中。金属离子、寡聚体配体和水化壳层的组合为自动QM簇生成带来了挑战，因为这三个组件都对准确描述系统至关重要。QuantumPDB将Ca2+–二糖复合物视为一个化学实体，并通过Voronoi镶嵌捕获第一溶剂层中的重要水分子。通过 `fill_dummy` 函数使用虚拟原子平滑接触层，比基于刚性距离的切割方法更能恰当地捕捉到关键的水化层（图7和支持信息文本S4）。如果没有虚拟原子，Voronoi划分在稀疏区域可能会导致远端残基被分配到簇中，从而在结果模型中产生不连续的区域（支持信息图S4）。QM簇模型用于由QuantumPDB生成的多残基中心系统。(左上角) C型凝集素Langerin (CD207, PDB ID: 3P5F) 与其钙离子及结合的甘露糖寡糖合并形成中心。(右上角) 环状肽药物Cyclosporine A与cyclophilin结合 (PDB ID: 1CWA)，整个11个残基的肽段被定义为中心。(左下角) 绿色荧光蛋白(GFP, PDB ID: 1EMA) 含有翻译后修饰的发色团(CRO)，该发色团由Ser65-Tyr66-Gly67三联体组成，被定义为中心。(右下角) Xylanase (XynII, PDB ID: 4HK8) 以多糖底物中的两个中心木糖单元为焦点，使模型集中在目标糖苷键上。作为中心的原子用黑色轮廓标出。一级、二级和三级结构球分别用灰色、浅蓝色和紫色表示。一级结构球以棒状表示突出特定配位残基(CD207)，或以表面表示展示大型多残基中心(CsA、GFP和XynII)周围的口袋形状。原子的颜色如下：氮蓝色；氧红色；硫黄色；氢白色；铁深橙色；中心残基的碳为橙色，其他位置的碳为灰色。配位键用紫色虚线表示。当中心被其他部分遮挡时，会显示切片视图。高分辨率图片下载 MS PowerPoint幻灯片

与碳水化合物链类似，用户也可以将氨基酸链定义为簇中心。例如，11个残基的环状肽药物Cyclosporine A (CsA) 与cyclophilin (PDB ID: 1CWA) 结合，可以将其整个定义为簇中心（图7）。此功能也适用于翻译后融合的氨基酸。一个关键的例子是绿色荧光蛋白(GFP, PDB ID: 1EMA) 的发色团，它由Ser65–Tyr66–Gly67三联体原位形成。这一过程生成了一个单一的共价嵌入残基(CRO，在其PDB文件中)，QuantumPDB可以利用[PHE_A64-CRO_A66-VAL_A68]的标记将其与相邻残基一起定义为簇中心（图7）。相反，QuantumPDB允许用户将QM模型集中在聚合物的催化相关亚基上。例如，在建模处理长碳水化合物链的酶（如xylanase (XynII, PDB ID: 4HK8) 时，用户可以使用QuantumPDB仅生成以催化相关部分为中心的簇。这样可以避免生成过大的QM簇模型，这些模型的大小与初始中心残基的数量成正比。例如，可以定义两个中间木糖单元为焦点，生成以目标糖苷键为中心的QM簇（图7）。

虽然我们构建了包含完整一级、二级和三级结构的模型作为示例，但资源有限的用户可能无法对这些最大的模型进行QM计算。在这种情况下，用户可以设置簇的总原子数限制，只有符合该限制的最接近的残基才会被包含在簇中。所有情况下生成的簇将具有不同的电荷，具体取决于每个球内包含的残基的性质。

4. 成功案例研究
4.1. 多样化酶数据集的整理
我们整理了一个酶及其天然底物结合的数据集，以展示QuantumPDB在大规模筛选中的实用性。为了系统地识别PDB中的酶结构，我们使用了酶的EC（EC number）分类系统。EC系统根据酶催化的化学反应提供层次化分类。我们使用PDB的RESTful应用程序编程接口（API）在2024年8月6日自动检索数据。该查询针对所有分类在七个主要EC类别（即氧化还原酶、转移酶、水解酶、裂解酶、异构酶、连接酶和转运酶）下的PDB结构。这得到了101,633个蛋白质结构的初步数据集（图8）。

为了获得所需的底物鉴定功能注释，我们使用了UniProt知识库成功识别了初始101,633个结构中的100,300个结构的蛋白质及其底物，并保留了这些结构（图8）。然后，我们精炼数据集，只保留适合建模具有高质量结构的催化活性复合物的结构。首先，为了丰富含有配体结合的结构，我们过滤掉了可能是无配体酶的条目。为了排除缺乏任何配体的蛋白质，我们扫描了HETATM记录，并排除了那些完全没有配体或其HETATM记录仅包含结晶缓冲液、离子、金属和常见辅因子的分子的结构。这一过滤步骤将数据集减少到了61,623个结构（图8）。其次，我们仅保留了通过X射线晶体学确定且分辨率小于3.0 ?且有相关出版物（即DOI）的结构，最终得到59,495个结构（图8）。第三，我们旨在去除任何异常大的蛋白质（即大于四分位数范围1.5倍或21,201个原子及以上），最终剩下57,580个结构（图8）。在结构质量过滤器之后，我们确保剩余结构不仅与异原子结合，更具体地与它们的天然底物结合。为了实现这种直接比较，我们使用了化学实体生物兴趣（ChEBI）标识符，这些标识符作为分子实体的词典。我们使用Rhea知识库（108-110）将晶体结构中存在的配体的化学身份与酶的注释生化反应中的已知反应物进行比较（图8和支持信息文本S8）。Rhea分析可以在PDB过滤之前进行，以缩小数据集范围，但这样需要在PDB服务器上下载结构两次，因此我们选择了更简洁的方法（支持信息图S5）。总体而言，我们的过程识别出了989个被认为与天然底物结合的完整酶结构，尽管这并不保证这些结构包含所有必要的辅因子/共底物或处于催化活性构象。从57,580个结构中大幅减少的原因是排除了包含抑制剂、产物或底物类似物的PDB条目，而不是Rhea反应参与者ID确定的天然催化活性底物。这一严格的化学身份匹配要求筛选出了一个高置信度的数据集，专门用于分析酶-天然底物的相互作用。

我们的过滤过程得到了一个高质量且多样化的989个酶的数据集，这些酶来自七个主要EC类别中的六个，非常适合随后的大规模量子力学分析，研究天然底物上的静电效应。在最初的七个EC类别中，转运酶（EC 7）是唯一未恢复的类别，这与它们在PDB中的低比例（1.1%）一致。虽然所有其他酶类别都有代表，但我们观察到它们的相对组成发生了变化：水解酶（EC 3），作为PDB中最大的酶类别，从38.2%减少到7.6%；而氧化还原酶（EC 1）从15.8%增加到36.2%（图8）。这种重新分布是我们对结构质量和明确底物分配进行严格筛选的直接结果。许多酶类别可以通过省略必需的辅因子或辅基（如O2或NAD+）与其天然底物一起结晶，从而避免转化。然而，水解酶使用水作为辅基，这无法从结晶环境中排除，因此很难捕获完整的酶-底物复合物。

4.2. 大规模QM分析
为了展示QuantumPDB在大规模研究中的实用性，我们研究了酶活性位点对其结合底物的影响。具体来说，我们使用QuantumPDB自动化QM簇构建，并使用GPU加速的ωPBEh-D3(BJ)/LACVP* DFT单点能量计算。我们的簇生成算法和用于蛋白质结构准备的外部库足以进行这种大规模分析。当对每个蛋白质进行串行处理时，簇生成的中位时间为2.79秒，Modeller和Protoss步骤的中位时间分别为7.09秒和9.05秒（支持信息图S6）。然后我们应用ωPBEh-D3(BJ)/LACVP* DFT单点能量计算来估计，在QM簇中包含更多QM原子如何通过电荷转移调节底物电荷并改变底物偶极矩（支持信息文本S9）。我们选择这些属性是因为它们是小型QM模型或具有最小QM区域的QM/MM模拟无法捕捉的因素。为了做出这一判断，我们评估了刚性底物在两种环境中的电子性质：(i) 在底物周围有132-812个原子的QM活性位点模型（包括一级和二级配位球），并进一步被MM点电荷包围；(ii) 仅在隐式水溶液中处理的底物（即ε = 80，支持信息文本S9）。在点电荷周围添加隐式溶剂或自由优化氢原子位置并不会显著改变计算出的性质（支持信息图S7）。为了量化电荷转移，我们使用qp.analysis接口与Multiwfn计算了在底物和QM活性位点上的实空间部分电荷，并对其总和进行了计算。底物的偶极矩作为完整QM系统中的一个片段来测量其被活性位点极化的程度（支持信息文本S10）。

在我们的1673个QM簇模型中（即每个酶至少有一个模型），我们观察到了广泛的电荷转移程度。我们将观察到的电荷转移量化为底物的形式电荷与其在酶活性位点上的部分电荷之和之间的差异（支持信息文本S10）。虽然对于相当一部分底物（即1673个中的381个，23.1%），活性位点计算出的净电荷与形式电荷的偏差小于0.1 e，但大多数底物的偏差更大。总体而言，蛋白质环境通常会减弱底物的电荷，无论其初始电荷是正还是负（图9）。尽管如此，任何初始电荷的底物的电荷都会发生变化。例如，中性底物在酶活性位点上可能变成正电荷或负电荷（图9）。电荷调制是显著的，底物的形式电荷与其在活性位点上的计算电荷之间的差异范围从+1.5到-3.2，平均值为-0.2（支持信息图S8）。尽管我们采用了范围分离的混合功能来减轻自相互作用误差并避免非物理电荷转移的影响（111-113），但电荷转移的定量性质可能对理论水平和基组敏感（支持信息文本S9-S10）。

图9显示了酶和底物之间的电荷转移。(左) 隐式溶液中底物的电荷与其在酶活性位点结合时的电荷的奇偶性图。实线黑色表示奇偶性，灰色虚线是最佳拟合线。展示了两个例子A和B。(中) 例子A显示了PDB ID: 5A60的活性位点，展示了底物的电荷转移。(右) 例子B显示了PDB ID: 6VI6的活性位点，也展示了底物的电荷转移。在A和B面板中，主要相互作用球表示为灰色表面，关键相互作用残基表示为棒状，次级球表示为蓝色表面。氢键表示为黄色虚线，配位键表示为紫色虚线。原子的颜色如下：碳（蛋白质），灰色；碳（底物），橙色；氮，蓝色；氧，红色；硫，黄色；磷，橙色；铁，深橙色；镁，绿色；氢，白色。高分辨率图像。下载MS PowerPoint幻灯片。

这种电荷稳定可以通过向底物提供或移除电子密度来实现。对于结合高度阴离子物种的酶，例如三磷酸隧道金属酶（TTM）ygiF（PDB ID：5A60），(114)，活性位点通过多达七个带正电荷的残基抽取电子密度来稳定高度带电的底物，这些残基直接与底物形成盐桥，使得底物的部分电荷为-1.8，这使得其电荷比形式电荷-5e更为中性（见图9和支持信息图S9）。在3-羟基 Anthranilate-3,4-二氧酶（HAO）（PDB ID：6VI6），(89）这种代表性的非血红素铁环裂解酶中，底物与一个路易斯酸Fe(II)中心配位，该中心从底物中抽取电子密度（见图9和支持信息图S9）。这导致活性位点的电荷从形式电荷-2变为-0.4，这对于降低O2活化的势垒及其后的芳香环裂解至关重要。(27,115)

为了更好地理解影响酶-底物电荷转移的因素，我们分析了活性位点组成与计算出的电荷差大小之间的关系，以及这与中性残基比例（FNR）的关系。(116) 如果活性位点的FNR大于0.8（对于第一和第二球体都适用），则将其标记为主要为中性（见支持信息文本S10）。进一步分析表明，即使在主要为中性的活性位点中也能观察到显著的电荷转移（见图10和支持信息文本S10）。这些主要为中性的活性位点通常也是疏水性的，这一趋势通过中性残基比例与平均Kyte–Doolittle疏水性得分之间的正相关关系（R2 = 0.56）得到证实（见支持信息图S10）。Kyte–Doolittle疏水性得分(117)是通过给疏水残基（例如亮氨酸）赋予高正值、给带电残基赋予负值得出的，因此疏水性活性位点的总得分更高。这种关系不仅限于完整的活性位点模型，在仅分析第一球体的FNR时也能观察到（见支持信息图S11）。两个有趣的例子是O-磷酸丝氨酸巯基化酶（PDB ID：3VSD）(118)和小鼠N1-乙酰多胺氧化酶（PDB ID：5MBX）(119)，它们结合的底物带有正电荷（形式电荷+3），尽管活性位点表面主要为中性，但仍然会引起电子密度的显著重新分布（见图10）。通过对酶活性位点中每个残基的Hirshfeld电荷进行总和分析，可以确认选定的模型区域主要为中性，只有少数区域有局部高电荷（即任何残基电荷>0.1 au，见图10）。因此，这种电荷转移源于许多中性极性和疏水残基的精确空间排列所产生的累积静电势，而不是强正电荷或负电荷。这些结果表明，仅基于残基组成或疏水性的酶-底物电荷转移假设是不充分的，需要显式的量子力学（QM）处理来全面考虑这些极化效应，即使基于残基类别的启发式方法预测这些效应很小。(120)

图10. 活性位点组成与底物电荷转移的关系。（左）所有球体的底物电荷差与FNR的散点图。点的颜色根据活性位点残基的平均Kyte–Doolittle疏水性来区分（蓝色表示更疏水；红色表示更亲水）。灰色虚线标记FNR = 0.8和电荷差 = 0.5作为一般的分界线。图中圈出了两个例子：A（PDB ID：3VSD）(118) 和 B（PDB ID：5MBX）(119)。（中）3VSD（A）的活性位点；（右）5MBX（B）的活性位点，底物以棒状表示，蛋白质表面根据每个残基的Hirshfeld部分电荷着色（颜色刻度：-1，红色；0，白色；+1，蓝色）。原子的颜色如下：碳，灰色；氮，蓝色；氧，红色；硫，黄色；磷，橙色；铁，深橙色；镁，绿色；氢，白色。高分辨率图像。下载MS PowerPoint幻灯片。

分子偶极矩所捕获的电荷分布也预计会影响底物的内在反应性。为了简化问题，我们计算了在没有酶环境和有酶环境情况下偶极矩大小的差异。我们发现，无论底物的电荷如何，酶活性位点都能一致性地减小底物的偶极矩大小（见支持信息图S12和支持信息文本S9）。这种偶极调制程度差异显著，并且不直接与底物的形式电荷相关，即通常情况下即使形式电荷不同，偶极矩大小也较大（见支持信息图S13–S15）。未来的分析可以通过显式计算偶极矢量来更全面地分析电荷重新分布的程度。此外，虽然我们计算出长底物（如两端带有电荷基团的(S)-3-羟基癸酰辅酶A（HSC）和连接的6-氨基己酸（ACA）分子的偶极矩非常大，但这些偶极矩也被酶活性位点一致性地减小了（见支持信息图S12）。我们还考虑了是否对底物电荷有显著影响的酶活性位点也会对偶极矩产生同等影响，但发现电荷差与偶极矩差之间没有相关性（皮尔逊相关系数r = 0.02，见支持信息图S16）。这表明不同的酶环境可以独立调节电荷和偶极矩，从而在增强反应性方面发挥不同作用。这可能是由于蛋白质不是均匀的介电环境，而是由能够独立调节这些电子特性的活性位点残基组成的高度特异性排列所导致的自然结果。

5. 结论

我们介绍了QuantumPDB，这是一个自动化的Python包，它将结构修复、质子化状态分配和QM聚类模型构建整合到一个可复现的工作流程中。该流程首先使用Modeller进行自动化结构修复，解决交替构象，然后使用Protoss在反馈循环中分配质子化状态，并识别和纠正空间冲突。该流程包括针对金属酶的基于启发式的校正，例如重新定位组氨酸咪唑环和脱质子化金属配位残基。虽然这个流程旨在具有通用性并包含了酶建模的最佳实践，但未来的工作可以通过将生成的模型与已推断或解析出氢位置的更高分辨率蛋白质结构进行对比来评估其准确性。QuantumPDB的一个关键特点是它使用基于Voronoi镶嵌的接触球体构建QM聚类模型，并通过虚拟原子进行规范化处理，以确保非球形活性位点的可靠划分，这些活性位点经常被基于距离的截断方法错误表示。QuantumPDB允许围绕灵活的中心定义构建聚类模型，该中心定义可以容纳单个原子、多残基复合物以及特定残基标识符的独特组合。通过识别同心接触球体，这种方法与多种方法兼容，包括QM/QM’和其他相关的ONIOM风格的方法。

为了展示其实用性，我们将QuantumPDB应用于一个包含989个全酶的数据集，对这些酶进行了DFT计算，得到了1,673个多球体QM聚类模型，这些模型采用了MM点电荷嵌入方法。对所得电子性质的分析显示，DFT模拟的酶环境会调节底物电荷，通常使其趋向中性，并一致性地减小底物的偶极矩大小相对于隐式水溶液的值。值得注意的是，即使在主要由中性残基组成的疏水性活性位点中，也会发生显著的电荷转移。这一结果强调了显式QM处理扩展酶活性位点的必要性，以捕捉这些效应。尽管如此，未来的工作还应考虑构象和动态的影响，特别是因为X射线晶体结构中解析出的天然配体可能并不处于正确的机械姿态。通过实现标准化、多球体的QM模型构建和分析，QuantumPDB为数据驱动的酶静电学和催化研究提供了一个强大的平台，从而弥合了详细单一系统QM研究与大规模生物信息学分析之间的差距。