modeller补全残基

PDB中经常存在缺失loop信息，可以使用Schrodinger软件的Protein Prepare Wizard 进行loop的修复。PDB中不存在缺失loop信息时，可以使用Schrodinger软件的Crosslink proteins进行loop的修复。但是当缺失loop残基数超过20时，这两个模块都无法正常运行，这时可以选择使用modeller修复loop。

值得注意的是，教程给出了两个脚本，一个是以缺失的蛋白PDB结构为模板进行同源模建修复loop，这样得到的蛋白骨架结构会和PDB原始结构存在小幅差异；第二个脚本是只行进loop的修复，而保持蛋白其他部分结构不变。

一个例子是PDB代码1qg8，它缺少残基134-136和218-231（见PDB文件中的REMARK 465行）。我们可以使用Modeller来 "填补 "这些缺失的残基，将原始结构（没有缺失的残基）作为一个模板，并使用完整的序列建立一个比较模型。

首先，我们获得具有已知结构的残基序列：

from modeller import *
# Get the sequence of the 1qg8 PDB file, and write to an alignment file
code = '1qg8'

e = Environ()
m = Model(e, file=code)
aln = Alignment(e)
aln.append_model(m, align_codes=code)
aln.write(file=code+'.seq')

这将产生一个序列文件，1qg8.seq。

>P1;1qg8
structureX:1qg8:   2 :A: 256 :A:undefined:undefined:-1.00:-1.00
PKVSVIMTSYNKSDYVAKSISSILSQTFSDFELFIMDDNSNEETLNVIRPFLNDNRVRFYQSDISGVKERTEKTR
YAALINQAIEMAEGEYITYATDDNIYMPDRLLKMVRELDTHPEKAVIYSASKTYHLNDIVKETVRPAAQVTWNAP
CAIDHCSVMHRYSVLEKVKEKFGSYWDESPAFYRIGDARFFWRVNHFYPFYPLDEELDLNYITEFVRNLPPQRNC
RELRESLKKLGMG*

从PDB的REMARKs或SEQRES记录中，我们知道了缺失的残基，所以现在我们可以在原始的1qg8结构（作为模板）和完整的序列之间做一个对齐，其中有对应于缺失残基的空白字符。我们将其放入一个新的对齐文件，alignment.ali。

>P1;1qg8
structureX:1qg8:   2 :A: 256 :A:undefined:undefined:-1.00:-1.00
PKVSVIMTSYNKSDYVAKSISSILSQTFSDFELFIMDDNSNEETLNVIRPFLNDNRVRFYQSDISGVKERTEKTR
YAALINQAIEMAEGEYITYATDDNIYMPDRLLKMVRELDTHPEKAVIYSASKTYHLN---DIVKETVRPAAQVTW
NAPCAIDHCSVMHRYSVLEKVKEKFGSYWDESPAFYRIGDARFFWRVNHFYPFYPLDEELDLNYIT---------
-----EFVRNLPPQRNCRELRESLKKLGMG*
>P1;1qg8_fill
sequence:::::::::
PKVSVIMTSYNKSDYVAKSISSILSQTFSDFELFIMDDNSNEETLNVIRPFLNDNRVRFYQSDISGVKERTEKTR
YAALINQAIEMAEGEYITYATDDNIYMPDRLLKMVRELDTHPEKAVIYSASKTYHLNENRDIVKETVRPAAQVTW
NAPCAIDHCSVMHRYSVLEKVKEKFGSYWDESPAFYRIGDARFFWRVNHFYPFYPLDEELDLNYITDQSIHFQLF
ELEKNEFVRNLPPQRNCRELRESLKKLGMG*

现在我们可以使用标准的Modeller 'LoopModel'class 来生成一个包含所有残基的模型，然后细化循环区域。

from modeller import *
from modeller.automodel import *    # Load the AutoModel class

log.verbose()
env = Environ()

# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']

a = LoopModel(env, alnfile = 'alignment.ali',
              knowns = '1qg8', sequence = '1qg8_fill')
a.starting_model= 1
a.ending_model  = 1

a.loop.starting_model = 1
a.loop.ending_model   = 2
a.loop.md_level       = refine.fast

a.make()

如果你不想要循环细化，只需使用 AutoModel class而不是LoopModel，并删除设置a.loop参数的三行。

请注意，循环建模默认情况下只对两个循环中较短的那个循环进行细化。您可以修改select_loop_atoms例程来完善两个环路，但是对于这种长的插入，您不可能得到好的结果。在这种情况下，你可能应该尝试为序列的这一部分找到另一个模板，或者如果你有理由相信这个插入不是一个环路的话，可以施加二级结构约束。

由于 AutoModel 或 LoopModel 将以您输入的 PDB 为模板建立一个比较模型，因此最终模型中的所有原子都有可能移动。如果你真的不希望非缺失的残基移动，你可以覆盖select_atoms方法，用类似下面的脚本只选择缺失的残基（注意残基编号偏离1，因为Modeller从A链的1开始给模型编号，而原始PDB从2开始编号）。

from modeller import *
from modeller.automodel import *    # Load the AutoModel class

log.verbose()
env = Environ()

# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']

class MyModel(AutoModel):
    def select_atoms(self):
        return Selection(self.residue_range('133:A', '135:A'),
                         self.residue_range('217:A', '230:A'))

a = MyModel(env, alnfile = 'alignment.ali',
            knowns = '1qg8', sequence = '1qg8_fill')
a.starting_model= 1
a.ending_model  = 1

a.make()