MODELLER
基本示例：
阴道毛滴虫乳酸脱氢酶的单模板建模
一个新型的关于乳酸脱氢酶的基因从阴道滴虫（TVLDH）中被确定下来。与其它的乳酸脱氢酶蛋白相比，这个类似蛋白与这个物种（阴道滴虫）的苹果酸脱氢酶（TvMDH）有着更高的相似性。我们假定通过相对近的时期的趋同进化，TVMDH从TvLDH中衍生了出来。为了研究这个序列的结构背景，并对一个通过定点突变来阐明在酶特异性的趋同进化的表观假设中的特异性改变的实验提出建议，我们对TvLDH与TvMDH建立了比较模型。天然的与突变后的酶表达后，对它们的活性进行了比较。
https://www.wendangku.net/doc/c711398076.html,LDH相关结构的搜索
首先，需要将目标TvLDH序列转化为MODELLER能读取的PIR格式（文件"TvLDH.ali").
第一行包含序列码，形式为">P1;code".用冒号分隔拥有十个域的第二行通常包含结构文件的信息，如果可用的话。只有这些区域中的两个被用于序列，"sequence"(表明这个文件包含一个结构未知的序列)和"TvLDH"（模式文件名）。该文件的其余部分包含TvLDH的序列，使用*标记它的结尾。使用标准的单字母氨基酸码。（注意必须使用大写形式；一些小写字母用于非标准残基。见Modeller中的
modlib/restyp.lib文件以获取更多的信息）
搜索已知结构的潜在相关序列通过MODELLER的profile.build()命令来执行。以下脚本，逐行执行以下操作(见文件"build_profile.py")：1.为这一建模运行初始化'environment'，通过创建一个新'environ'对象。
几乎所有的MODELLER脚本都需要这一步，因为新的对象（在这里我们称之为'env'，但你可以将它命名为任何你喜欢的）对构建大多数其他有用的对象是必须的。
2.创建一个新的'sequence_db'对象，命名为'sdb'。'sequence_db'被用于
容纳蛋白质序列的大型数据库。
3.读取一个包含序列相似性在95%的非冗余PDB序列的text格式文
件到sdb数据库中。这一序列在"pdb_95.pir"文件中，和之前创建的序列一样，这个文件也是PIR格式。
4.写一个二进制机器特异性文件，包含上一步骤中所读的所有序列。
5.读取二进制格式文件。注意，如果你打算使用一个数据库多次，
你应该只使用前两个步骤一次，以生成二进制数据库。在随后的运行中，你可以忽略这两个步骤并且直接使用二进制文件，因为读取二进制文件比读取PIR文件快很多。
6.产生一个新的'alignment'对象，命名为'aln'，从文件"TvLDH.ali"中
读取我们的二进制序列"TvLDH"，并且将其转化到一个配置文件'prf'中。配置文件包含着与序列比对相似的信息，但是更加的简洁并且有利于数据库搜索。
7.为我们的二进制配置文件'prf'搜索序列数据库'sdb'。序列数据库中
所匹配的被添加到配置文件中。
8.写一个二进制序列和它的同源物的配置文件（见文件
"build_profile.prf").等效信息也被以标准序列格式写出。
注意，因为这个脚本是通过Python编程语言所写的，它使用Modeller-特异性命令，所以在你的命令行中应该通过使用一个和下面相似的命令来运行该脚本
mod9v1 build_profile.py
注意，"mod9v1"脚本运行Modeller版本9v1，使用不同版本的Modeller就将"mod9v1"改为相应的版本。
profile.build()命令拥有许多选择项。在这个例子中rr_file被设置使用BLOSUM62相似矩阵（Modeller发布中所提供的文件"blosum62.sim.mat"）。因此，参数matrix_offset和gap_penalties_1d 被设置成对BLOSUM62 矩阵是适当的。例如，通过设置参数n_prof_iterations等于1我们将仅运行一个搜索重复。这样，不需要检查配置文件的偏差(check_profile set to False)。最后，参数max_aln_evalue被设置成0.01，说明e值小于或等于0.01的序列将被包含在最终的配置文件中。
2.筛选一个模板
"build_profile.py"脚本的输出被写入到"build_profile.log"文件。MODELLER总是产生一个log文件。log文件中的错误和警告能够通过分别搜索"_E>" and "_W>" 附加条件而发现。MODELLER也将配置文件写入到text格式的"build_profile.prf"文件。输出文件的一个摘录（忽略对齐序列）接下来能被看到。前6个注释行表明MODELLER 中用来构建配置文件的输入参数。后面的注释行相当于通过profile.build()所检测的相似性。
profile.build()输出中最重要的列是第二，第十，第十一以及第十二列。第二列报告了和目标序列相比较的PDB序列的代码。每一行的PDB 代码是一组PDB序列的代表，这些序列彼此拥有95%或更高的序列相似性并且拥有少于30个残基或30%序列长度差异。第十一列报告了TvLDH和一个通过序列长度标准化了的PDB序列（第十列所显示）之间的序列一致性百分比。通常，一个序列的一致性值在大约在25%以上表示一个潜在模板，除非序列太短（例如，少于100个残基）。第十二行中通过序列的e值来更好地检测显著性序列。在这个例子中，六个PDB序列对e值为0二进制序列显示出极显著相似性。正如所预期的，所有的采样点和苹果酸脱氢酶相一致(1bdm:A, 5mdh:A, 1b8p:A, 1civ:A, 7mdh:A, and 1smk:A)。在六个相似结构中为了筛选出最适合我们二进制序列的模板，我们将使用https://www.wendangku.net/doc/c711398076.html,pare_structures()命令来评估这些可能模板间的结构和序列相似性。(文件"compare.py").
在这种情况下，我们创建了一个（最初是空的）对齐对象'aln'然后使用一个Python 'for'循环来引导MODELLER来读取每一个PDB文件。（注意，为了使其能够实现，和这个脚本相同的目录下你必须有所有的PDB文件，从PDB网站上或者从该页顶部链接归档下载）。我们用model_segment参数要求仅对一个单链从每个PDB文件中读取（默认情况下，所有链都会从文件中读取）。当每一个结构都被读入，我们使用append_model法将结构添加到序列比对中。
在循环的结尾，所有的结构都在序列对比中，但是它们不是理想的彼此对齐（append_model生成了简单的1：1连续序列）。所以，我们通过使用malign来计算一个多序列比对的方法来改进这个比对。malign3d命令于是执行六个三维结构的迭代最小二乘叠加，使用多序列比对作为它的起点。该compare_structures命令根据malign3d构建的对比来比较这些结构。它并没有执行一个比对，而是计算原子位置和距离之间的RMS和DRMS偏差，主链与侧链二面角之间的不同，序列一致性的百分比，以及多个其它检测。最后，id_table命令写一个拥有成对序列线距的文件，它能够直接作为dendrogram命令的输入文件（或PHYLIP软件包中的聚类程序）。dendrogram命令计算一
个来自从成对距离的输入矩阵的聚类树，它帮助可视化这些候选模板间的差异。log文件中的节选在下面显示（文件"compare.log"）。
以上的对比显示1civ:A 和7mdh:A在序列和结构上最相似。然而，7mdh:A拥有更好的结晶分辨率（2.4? ：2.8?），排除1civ:A。第二组结构(5mdh:A, 1bdm:A, and 1b8p:A)拥有相同的相似性。在这个组中，5mdh:A有最低的分辨率，仅考虑1bdm:A, and 1b8p:A。1smk:A是所有可能的模板中最不同的结构。然而，和所查询序列相比它是序列相似性最低的一个。我们最终挑选了超过1b8p:A and 7mdh:A的1bdm:A，因为它有更好的晶体结构R-factor（16.9%）和与所查序列之间更高的全序列一致性（45%）。
3.调整TvLDF和模板一致
MODELLER中的align2d()命令是调整TvLDF的序列与1bdm:A的结构一致的好方法。尽管align2d()基于动态规划算法，它和标准的序列-序列比对方法是不同的，因为当构造一个比对时它考虑了来自模板的结构信息。这一任务是通过一个可变间隙罚函数来实现的，这一函
数趋向于在溶剂暴露和弯曲区域，二级结构片段之外，在空间上相互靠近的两个位置中放置间隙。结果，序列对比错误和使用标准序列对比技术所产生的错误相比减少了将近三分之一。当序列间的相似度下降及缺口数增多时这一提升变得尤为重要。在这个例子中，模板—目标之间的相似度太高以至于几乎所有的使用合理参数的对比方法将产生相同的对比结果。以下的MODELLER脚本使得文件"TvLDH.ali"中的TvLDH序列与文件"1bdm.pdb"中的1bdm:A结构相一致。(文件"align2d.py").
在这个脚本中，我们又创建了一个'environ'对象来作为后面命令的输入。我们创建一个空白比对'aln'，然后一个新的蛋白质模型'mdl'，到我们读取的1bdm PDB结构文件中的A链片段中。append_model()命令将这一模型的PDB序列传输到对比中，将其命名为"1bdmA" (align_codes)。然后使用append()命令将"TvLDH.seq"文件中的"TvLDH"序列添加到对比中。align2d()命令接着被用于执行对比这两个序列。最后，对比结果被写出为两种格式，PIR ("TvLDH-1bdmA.ali") and PAP ("TvLDH-1bdmA.pap")。PIR格式再随后的建模阶段中能够被MODELLER使用，而PAP对比格式更易于可见分析。由于目标-模板间的高相似性，对比结果中仅存在很少的缺口。在PAP格式中，
所有的相同位置都用“*”标记(文件"TvLDH-1bdmA.pap")。
4.建模
一旦一个目标序列与模板间的对比分析被构建，MODELLER使用automodel类来完全自动地计算目标序列的3D模型。以下的脚本将会生成基于1bdm:A模板结构和文件"TvLDH-1bdmA.ali"中的序列对比的五个相似TvLDH模型。(文件"model-single.py").
第一行装载了automodel class并将它准备使用。然后我们创建了一个automodel对象，称之为‘a’并且设置参数来引导模型建立步骤。alnfile 命名包含PIR格式的目标-模板序列对比文件。knowns定义了alnfile 中的已知模板结构("TvLDH-1bdmA.ali")。sequence定义了alnfile中的目标序列的名称。assess_methods要求一个或者更多的评估分数（在
接下来的章节中详细讨论）。starting_model和ending_model定义了被计算的模板数(它们的索引将从1运行到5)。文件中的最后一行叫做make法，它实际上计算了该模型。
最重要的输出文件是"model-single.log",它报告了警告，错误，以及其他有用的信息，包括在最终模型中保留违犯的用于建模的输入限制。这一log文件中的最后几行显示如下。
如你所见，log文件总结了所有建立的模型。对于每一个模型，它列出了文件名，这些文件名包含了PDB中的模型坐标。这些模型能够被任何读取PDB格式的程序所查看，比如Chimera。log还显示了每一个模型的分数，我们将在下面进一步讨论（注意，实际数量在你的机器上也许稍有不同-这点毋庸担忧）
5.模型评估
如果多个模型被用于计算相同的目标，最优模型能够通过多种方法被筛选出来。比如，你可以挑选拥有最低MODELLER目标函数或DOPE 或SOAP评估分数，或者拥有最高GA341评估分数的模型，这些报告在上面的log文件的末尾。（目标函数，molpdf也被计算了，并且同样报告在每个生成的PDB文件中。DOPE，SOAP,以及GA341数，或者其它任何评估分数，只有当你将他们列在assess_methods中时才被计算。为了计算SOAP分数，你首先需要在SOAP网站上下载SOAP-Protein潜在文件，然后通过删除‘#’来在model-single.py中
取消备注SOAP相关行。）molpdf,DOPE,及SOAP分数不是绝对的评价标准，在这个意义上它们仅能被用于排列计算自相同序列对比的模型。其它分数是可转移的。比如GA341分数总是从0.0（最低）到1.0（本体型）；所以从‘bad’模型中区分‘good’模型时GA341没有DOPE或SOAP好。
一旦最终模型被筛选出来，它能够被很多的方式进一步评估。用于模型评估的程序能够在这页的MODEL EV ALUATION中被发现。
在进行该模型的外部评价之前，应该为运行错误("model-single.log")和限制违犯检查来自建模运行的log文件。文件"evaluate_model.py"评价了一个拥有DOPE势的输入模型。（注意，这里我们任意挑选了第二个生成模型-你可以尝试其他模型）
在这个脚本中我们使用complete_pdb脚本来读入一个PDB文件并准备它做能量计算（该自动化允许了这些可能性，PDB 文件拥有无标准顺序的原子，或者具有不同的原子子集，比如包括氢在内的所有原子，但是MODELLER仅仅使用重原子，反之亦然）。我们创建了一个所有原子的筛选，因为多数MODELLER能量函数能够在一个模型原子的子集上运行。DOPE能量这时候通过assess_dope命令来计算，并且我们又请求了一个能量profile，平滑过一个15个残基的窗口，并通过在每一个残基上运行的限制数量将其标准化。这一profile被
写入文件"TvLDH.profile"中，它能够被用做一个成像程序的输入。比如，它能够被GNUPLOT使用'plot "TvLDH.profile" using 1:42 with lines'命令来标绘。或者，你可以使用plot_profiles.py脚本，它包含在tutorialzip文件中用来通过Python matplotlib安装包来标绘profiles。
GA341分数，作为使用PROCHECK程序的额外分析，确认了TvLDH.B99990001.pdb模型是一个合理的模型。然而，被标绘的DOPE分数profile（以下）显示了目标序列C-末端的残基90和100和长螺旋之间的长活性位点环的相对高能量的区域。（注意，我们在模板profile之上叠加了模型profile-与序列对比中的缺口相一致，标绘中的缺口能够被看到。记住这一分数不是绝对的，所以我们不能在两者间做一个数量上的比较。然而，通过比较两个profile的大致形状我们能得到一个我们的输入序列对比的质量概念-如果一个是明显地转移到另一个，很有可能这一对比也是从正确的一个转移而来）
该长环和220-250区域相交互，它组成了另一半的活性位点。后面的部分在模板中被很好地解决并且可能正确地模拟了目标结构，但是由
于和90-100区域的不利非键相互作用，它也被DOPE报告为高能量。一般，DOPE所推测的一个可能错误不一定是一个真正的错误，尤其当他凸显一个活性位点或一个蛋白质-蛋白质相互作用界面。然而，在这种情况下，在模板结构中同一活性位点环拥有更好的profile，这加强了这一评估-这一个模型可能在这一活性位点区域是不准确的。这一问题在advanced modeling tutorial中解决
高级示例：
基于多模板，循环优化及用户指定限制的蛋白配体复合物的建模
在这个例子中我们将不会一步一步地描述所有的MODELLER命令。请参考tutorial中的basic-example以详细了解。
这一建模的重要目的是帮助了解所要建模蛋白的功能。检查1bdm:A模板结构发现环93-100，是该酶功能上最重要的部分之一，是无序的并且在PDB结构中不存在。极有可能该长活性位点在没有配体的情况下是灵活的并且在衍射图中看不到。模板坐标的不可靠性和MODELLER不能够模拟长的插入片段，使得该环在TvLDH中被不准确建模，如同DOPE profile所显示的那样。
既然我们有兴趣了解这两个相似蛋白质的不同特异性，我们需要建立精准的模型。因此，我们需要找到增加这些模型精确性的新策略。在这个例子中，我们将探讨三个不同的方法：
使用多模板
使用ab-initio法对这个环建模
使用一个已知配体结合到结合位点来建模
Multiple templates
苹果酸脱氢酶1bdm的结构已经被整合到DBAli database 中的家族fm00495中，有4个成员（2mdh:A, 2mdh:B. 1b8p:A and 1bdm:A）。在MODELLER中用salign()命令创建的多比对被用在DBAli中创建家族的多结构比对。这个比对可以从DBAli数据库中下载，或者你可以用`salign.py'文件在你的电脑上计算。
用append_model命令从PDB文件中阅读所有的序列，接着用多次salign命令建立最初的粗略的比对，再通过更多的信息来提高它。比对结果用PIR和PAP两种格式输出，再一次使用salign命令来计算质量得分。
检查多结构比对结果之后，明显的发现2mdh中的B链含有一些与众不同的赖氨酸LYS残基。优化过程中PDB文件的标题表明蛋白的序列式未知的，并且在结构上很难识别大多数的残基。因此，把2mdh:B 从多结构比对结果中移除出去。
如同在tutorial中的basicexample一样，接下来我们需要将我们的查询序列同模板结构进行比对。我们再一次使用salign()命令（文件`align2d_mult.py'）来进行上述任务。我们设置align_block参数等于在模板比对中的结构的数目len(aln), (i.e. 3)，并且需要一个双序列比对，因为我们不想改变模板之间已经存在的比对。通过设置gap_function，我们需要用依赖结构的空位处罚，为这3个序列使用结构信息。只有序列信息用于最终的TvLDH序列。
接下来，使用`model_mult.py' 文件建立基于多模板比对的目标序列TvLDH的新模型
最后，我们用DOPE电势来评估新模型的坐标，使用文件`evaluate_model.py'。
模型评估结果显示不确定环通过使用多结构模板比对已经有所改进。模型总DOPE得分从-38999.7降到-39164.4。MODELLER也可以利用三个模板环区的变异性来创建一个更精确的环构象。但是，在该序列C末端的275残基附近的一个环构象的DOPE得分比用单模板得到的模型的DOPE得分还高
我们将用MODELLER中的loopmodel分类法优化273到2283残基环的构象。我们将用前一个例子中的模型1作为初始结构来优化环。你可以在loop_modeling子目录中重新命名为`TvDLH_mult.pdb'中发现这个结构。
Loop refining
环最优化方法依赖一个得分作用和适合于环模建的最优化方案。如果你用loopmodel而不是automodel命令，它将自动的优化比较模型，看下图的例子。
在这个例子中，loopmodel类用来优化已存在的坐标文件的区域。要注意的是这个例子重新定义了loopmodel.select_loop_atoms程序。在这种情况下是必要的，因为缺省选择程序选择比对结果中的所有缺口进行优化，而在这种情况下没有比对结果是可用的。即使你有一个比对结果，如果你想选择一个不同的区域进行优化，仍然要重新定义程序。注意为了节省时间，我们将仅建立10个不同的独立的最优化环结构，通过设置loop.ending_model参数为10。下图展示了10种环模型构造区域的重叠。我们事先标记的蓝色、绿色和红色环分别代表了通过DOPE得分评估的最好和最坏的环构造
文件`model_energies.py' 通过Python命令为所有建立的适于环的模型计算DOPE得分。优化最好的能量环是第8个模型，DOPE总得分39099.1。用`evaluate_model.py'计算的能量曲线图如下：
用ab-initio优化环得到的总体DOPE得分仅有非常小的增加。然而环状区域的DOPE得分有少量的增加。因此，我们用最优的结构继续进行下一步，在配位体目录中最优结构重命名为`TvLDH-loop.pdb'。环优化一个最准确的方法需要模拟上百的独立的构象和他们的分类，从中选择最具代表性的环结构。
Modeling ligands in the binding site
模拟结合位点的配基
1emd（来自大肠杆菌的苹果酸脱氢酶）在PDB中被证实。尽管1emd 序列和TvLDH序列只有32%的一致性，但是环区的活性位点和他们所处的环境是更加保守的。1emd环的90到100残基结构已经很好的解决了。此外，1emd溶解在柠檬酸盐基底类似物和NADH辅因子中。新的比对结果在PAP格式中展示如下（文件`TvLDH-1emd_bs.pap'）
依据一个编辑过的1emd结构修饰的比对序列(1emd_bs)作为第二个模板。比对是对应于一个模板，该模板基于在1emd_bs环的活性位点和TvLDH_model（上一步中最好的模型）折叠区的剩余部位。活性