建模后特定残基的进一步优化

lrf1980

问题：在modeller的官方网站上，有针对loop的优化脚本，当我们建好模型后，我们可以进一步对loop做优化，来完善我们的模型，但是在很多情况下，可能活性口袋并不在loop里，又或者说，对于要研究的蛋白质，通过文献调研和数据收集，手头上已经有比较多的突变数据，指出某些特定的残基对蛋白质活性影响很大。我们如何能正对这些特定的残基做进一步的优化，或者说在自己建的模型中，能够尽量多的针对这些残基产生尽可能多的构象，用于后续评估？
思考：在modeller的网站manual 里，有这么一个部分，叫做refining only part of the model.从这部分的功能解释上来看，应该可以用于解决上面提到的问题。http://salilab.org/modeller/manual/node23.html

The automodel class contains a automodel.select_atoms() functionwhich selects the atoms to be moved during optimization. By default, theroutine selects all atoms, but you can redefine it to select any subset ofatoms, and then onlythose atoms will be refined. (To redefine the routine, it is necessary tocreate a `subclass' of automodel, here called MyModel, whichhas the modified routine within it. We then use MyModel in place ofautomodel. The select_atoms routine should return aselection object; see Section 6.9 for furtherinformation.)

下面是实现该功能脚本 model-segment.py的内容：

# Comparative modeling by the automodel class
#
# Demonstrates how to refine only a part of the model.
#
# You may want to use the more exhaustive "loop" modeling routines instead.
#
from modeller import *
from modeller.automodel import *    # Load the automodel class
log.verbose()

# Override the 'select_atoms' routine in the 'automodel' class:
# (To build an all-hydrogen model, derive from allhmodel rather than automodel
# here.)
class MyModel(automodel):   
def select_atoms(self):        
# Select residues 1 and 2 (PDB numbering)        
return selection(self.residue_range('1:', '2:'))      
# The same thing from chain A (required for multi-chain models):        
# return selection(self.residue_range('1:A', '2:A'))        
# Residues 4, 6, 10:        
# return selection(self.residues['4'], self.residues['6'],        
#                  self.residues['10'])        
# All residues except 1-5:        
# return selection(self) - selection(self.residue_range('1', '5'))
env = environ()
# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']
# selected atoms do not feel the neighborhood
env.edat.nonbonded_sel_atoms = 2
# Be sure to use 'MyModel' rather than 'automodel' here!
a = MyModel(env,  
          alnfile  = 'alignment.ali',     # alignment filename  
          knowns   = '5fd1',              # codes of the templates           
          sequence = '1fdx')              # code of the target
a.starting_model= 3                # index of the first model
a.ending_model  = 3                # index of the last model                             
                                   # (determines how many models to calculate)
a.make()                           # do comparative modeling
从该脚本的内容来看，我们可以用这个脚本来实现：
1. 对选定范围内的氨基酸refine，比方说从20号残基到50号残基。
return selection(self.residue_range('20:', '50:'))
2. 对指定的某些氨基酸refine， 比方说，12， 15， 32， 34， 45。
return selection(self.residues['12'], self.residues['15'], self.residues['32']，self.residues['34'], self.residues['45'])
3. 再着就是除某些氨基酸外的其它所有氨基酸，或者说建模时，这几个氨基酸构象固定的refine，比方说，对10-15号氨基酸残基外的其它氨基酸refine。
return selection(self) - selection(self.residue_range('10', '15'))

当我们了解这些规则时，我们来看看是否可以实现我们的目的。我选4DXD这个蛋白质来做例子。这个蛋白质在pdb已经有晶体结构，我就假设这个pdb
是我们建好的模板，或者说是我们的模板。

1.首先我将我们目标蛋白质模型命名4DXD.pdb，原来的4DXD命名成4DXD1作为模板，然后我生成一个ali文件，这时候我的模板跟我的目标蛋白
的序列应该是一模一样的。（如果是直接从模板开始做这一步，那就不可能是100%一致）

2.然后我们再用align2d.py这个脚本生成一个4DXD-4DXD1.ali 和4DXD-4DXD1.pap文件。

3. 然后我们需定需要refine的氨基酸，这里我选的是9pc口袋下的一段氨基酸序列，189-204号残基

4. 下面是我的model-segment.py脚本内容：
# Comparative modeling by the automodel class
#
# Demonstrates how to refine only a part of the model.
#
# You may want to use the more exhaustive "loop" modeling routines instead.
#
from modeller import *
from modeller.automodel import *    # Load the automodel class

log.verbose()

# Override the 'select_atoms' routine in the 'automodel' class:
# (To build an all-hydrogen model, derive from allhmodel rather than automodel
# here.)
class MyModel(automodel):
    def select_atoms(self):
        # Select residues 1 and 2 (PDB numbering)
        return selection(self.residue_range('189:', '204:'))

        # The same thing from chain A (required for multi-chain models):
        # return selection(self.residue_range('1:A', '2:A'))

        # Residues 4, 6, 10:
        # return selection(self.residues['4'], self.residues['6'],
        #                  self.residues['10'])

        # All residues except 1-5:
        # return selection(self) - selection(self.residue_range('1', '5'))


env = environ()
# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']
# selected atoms do not feel the neighborhood
env.edat.nonbonded_sel_atoms = 2

# Be sure to use 'MyModel' rather than 'automodel' here!
a = MyModel(env,
            alnfile  = '4DXD-4DXD1.ali',     # alignment filename
            knowns   = '4DXD1',              # codes of the templates
            sequence = '4DXD',
            assess_methods=(assess.DOPE,
                              #soap_protein_od.Scorer(),
                              assess.GA341))              # code of the target

a.starting_model= 1                # index of the first model
a.ending_model  = 100                # index of the last model
                                   # (determines how many models to calculate)
a.make()                           # do comparative modeling

我在这个脚本里加了对每个模型的DOPE 和GA341评估。

5. 运行这个脚本，然后用pymol把100个模型重叠在一起，我们就可以观察到所有的模型里，只有189-205这个区域的氨基酸构象是变化的。然后在根据实际情况评估和选择合适的一个模型做后续的研究。

6. 附件中是所有的相关文件，有兴趣的同志可以自己下载下来运行一下。

秋水若寒

请教一下：用上述所说的，是否能够解决以下问题：
1、评估模型时，有一个区域中氨基酸得分低，是否可以进行结构优化

这些程序看不太明白。是否还有其他的优化方法？

秋水若寒

又看了一遍编程，还是没明白该如何操作，设置

lrf1980

秋水若寒 发表于 2015-1-28 22:25
又看了一遍编程，还是没明白该如何操作，设置

把你的脚本贴出来给我看看，还有就是在SAVES上评价的结果也贴出来