蛋白结构预测生物信息杂谈

#分子模拟#同源建模从入门到精通(四)

2017-05-12  本文已影响0人  生信杂谈

  上一期刚刚讲了一种处理部分模板结构质量差的解决方法,不知道各位看官有没有看过瘾。今天我们要将两个相对比较简单的方法进行模型质量的优化提升。

loop环精炼

loop环优化依赖于评分功能,并且优化算法仅适用于loop环,可以发现上期处理完以后蛋白273-283附近结果相对不理想,检查是否能够再次进行优化。(chimera显示的二级结构有一些不准确,我们在这里不做探讨,实际上273-283是完整的loop环)

在这里我们使用loopmodel,而不是automodel,具体如下:
<font color="red" >相信这么多天的推文,大家对同源建模已经有了一定的基础,所以后面我们就仅进行例子查看,并讲解主要要注意的地方,不再写一个伪代码,一个示例代码</font>

# Loop refinement of an existing model
from modeller import *
from modeller.automodel import *

log.verbose()
env = environ()

# directories for input atom files
env.io.atom_files_directory = './:../atom_files'

# Create a new class based on 'loopmodel' so that we can redefine
# select_loop_atoms (necessary)
class MyLoop(loopmodel):
    # This routine picks the residues to be refined by loop modeling
    def select_loop_atoms(self):
        # 10 residue insertion 
        return selection(self.residue_range('273', '283'))

m = MyLoop(env,
           inimodel='TvLDH-mult.pdb', # initial model of the target
           sequence='TvLDH')          # code of the target

m.loop.starting_model= 1           # index of the first loop model 
m.loop.ending_model  = 10          # index of the last loop model
m.loop.md_level = refine.very_fast # loop refinement method; this yields
                                   # models quickly but of low quality;
                                   # use refine.slow for better models

m.make()

File: loop_modeling/loop_refine.py

这个例子我们创建了一个Myloop的子类,用来继承loopmodel类,并返回要精炼的loop环范围,上面的示例为273至283的区域. refine.very_fast使得建模更快,但是相对质量不佳,建议使用注释后面的refine.slow。最后再进行DOPE打分。

基于配体建模

个人理解基于配体建模最佳的方式是目标蛋白与模板蛋白有相似或者相同的小分子配体,这样建出来的模型效果相对较好。如果我的理解不对,欢迎指出问题:
首先我们按照同源建模入门到精通(二)的方法对loop环优化后的配体进行模板比对,将带有配体的1emd比对上去,或者手动比对上去也可以。结果如下:

 _aln.pos            10        20        30        40        50        60
TvLDH        MSEAAHVLITGAAGQIGYILSHWIASGELYGDRQVYLHLLDIPPAMNRLTALTMELEDCAFPHLA 
TvLDH_model  MSEAAHVLITGAAGQIGYILSHWIASGELYGDRQVYLHLLDIPPAMNRLTALTMELEDCAFPHLA 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos       70        80        90       100       110       120       130
TvLDH        GFVATTDPKAAFKDIDCAFLVASMPLKPGQVRADLISSNSVIFKNTGEYLSKWAKPSVKVLVIGN 
TvLDH_model  GFVATTDPKAAFKDIDCAFLVASMPLKPGQVRADLISSNSVIFKNTGEYLSKWAKPSVKVLVIGN 
1emd         ----------------------GVRRKPGMDRSDLFNVN-------------------------- 
 _consrvd                              ***  * **   *


 _aln.pos           140       150       160       170       180       190
TvLDH        PDNTNCEIAMLHAKNLKPENFSSLSMLDQNRAYYEVASKLGVDVKDVHDIIVWGNHGESMVADLT 
TvLDH_model  PDNTNCEIAMLHAKNLKPENFSSLSMLDQNRAYYEVASKLGVDVKDVHDIIVWGNHGESMVADLT 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos      200       210       220       230       240       250       260
TvLDH        QATFTKEGKTQKVVDVLDHDYVFDTFFKKIGHRAWDILEHRGFTSAASPTKAAIQHMKAWLFGTA 
TvLDH_model  QATFTKEGKTQKVVDVLDHDYVFDTFFKKIGHRAWDILEHRGFTSAASPTKAAIQHMKAWLFGTA 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos           270       280       290       300       310       320
TvLDH        PGEVLSMGIPVPEGNPYGIKPGVVFSFPCNVDKEGKIHVVEGFKVNDWLREKLDFTEKDLFHEKE 
TvLDH_model  PGEVLSMGIPVPEGNPYGIKPGVVFSFPCNVDKEGKIHVVEGFKVNDWLREKLDFTEKDLFHEKE 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos      330
TvLDH        IALNHLAQGG/.. 
TvLDH_model  IALNHLAQGG/-- 
1emd         ----------/.. 
 _consrvd

比对好后我们在目标序列与带有配体的序列后加"/",表示序列结束,然后加"."表示匹配配体,该示例后面加了两个"."表示匹配两个配体。

最后进行建模:

from modeller import *
from modeller.automodel import *

class MyModel(automodel):
    def special_restraints(self, aln):
        rsr = self.restraints
        for ids in (('NH1:161:A', 'O1A:336:B'),
                    ('NH2:161:A', 'O1B:336:B'),
                    ('NE2:186:A', 'O2:336:B')):
            atoms = [self.atoms[i] for i in ids]
            rsr.add(forms.upper_bound(group=physical.upper_distance,
                                      feature=features.distance(*atoms),
                                      mean=3.5, stdev=0.1))

env = environ()
env.io.hetatm = True
a = MyModel(env, alnfile='TvLDH-1emd_bs.ali',
              knowns=('TvLDH_model','1emd'), sequence='TvLDH')
a.starting_model = 1
a.ending_model = 5
a.make()

File: ligand/model-multiple-hetero.py

建模脚本需要注意的是``env.io.hetatm = True`表示非标准”残基"也可以被读取,这样才能读取到配体。

分享到这里其实同源建模用的最普遍的内容已经讲完,接下来准备和大家一起学习迭代建模和建模更深的一些常用函数。这一系列就将先告一段落,等到时候研究rosetta的时候或者看到关于建模的一些好文章再继续撰写该类型。
在这里我们使用loopmodel,而不是automodel,具体如下:
<font color="red" >相信这么多天的推文,大家对同源建模已经有了一定的基础,所以后面我们就仅进行例子查看,并讲解主要要注意的地方,不再写一个伪代码,一个示例代码</font>

# Loop refinement of an existing model
from modeller import *
from modeller.automodel import *

log.verbose()
env = environ()

# directories for input atom files
env.io.atom_files_directory = './:../atom_files'

# Create a new class based on 'loopmodel' so that we can redefine
# select_loop_atoms (necessary)
class MyLoop(loopmodel):
    # This routine picks the residues to be refined by loop modeling
    def select_loop_atoms(self):
        # 10 residue insertion 
        return selection(self.residue_range('273', '283'))

m = MyLoop(env,
           inimodel='TvLDH-mult.pdb', # initial model of the target
           sequence='TvLDH')          # code of the target

m.loop.starting_model= 1           # index of the first loop model 
m.loop.ending_model  = 10          # index of the last loop model
m.loop.md_level = refine.very_fast # loop refinement method; this yields
                                   # models quickly but of low quality;
                                   # use refine.slow for better models

m.make()

File: loop_modeling/loop_refine.py

这个例子我们创建了一个Myloop的子类,用来继承loopmodel类,并返回要精炼的loop环范围,上面的示例为273至283的区域. refine.very_fast使得建模更快,但是相对质量不佳,建议使用注释后面的refine.slow。最后再进行DOPE打分。

基于配体建模

个人理解基于配体建模最佳的方式是目标蛋白与模板蛋白有相似或者相同的小分子配体,这样建出来的模型效果相对较好。如果我的理解不对,欢迎指出问题:
首先我们按照同源建模入门到精通(二)的方法对loop环优化后的配体进行模板比对,将带有配体的1emd比对上去,或者手动比对上去也可以。结果如下:

 _aln.pos            10        20        30        40        50        60
TvLDH        MSEAAHVLITGAAGQIGYILSHWIASGELYGDRQVYLHLLDIPPAMNRLTALTMELEDCAFPHLA 
TvLDH_model  MSEAAHVLITGAAGQIGYILSHWIASGELYGDRQVYLHLLDIPPAMNRLTALTMELEDCAFPHLA 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos       70        80        90       100       110       120       130
TvLDH        GFVATTDPKAAFKDIDCAFLVASMPLKPGQVRADLISSNSVIFKNTGEYLSKWAKPSVKVLVIGN 
TvLDH_model  GFVATTDPKAAFKDIDCAFLVASMPLKPGQVRADLISSNSVIFKNTGEYLSKWAKPSVKVLVIGN 
1emd         ----------------------GVRRKPGMDRSDLFNVN-------------------------- 
 _consrvd                              ***  * **   *


 _aln.pos           140       150       160       170       180       190
TvLDH        PDNTNCEIAMLHAKNLKPENFSSLSMLDQNRAYYEVASKLGVDVKDVHDIIVWGNHGESMVADLT 
TvLDH_model  PDNTNCEIAMLHAKNLKPENFSSLSMLDQNRAYYEVASKLGVDVKDVHDIIVWGNHGESMVADLT 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos      200       210       220       230       240       250       260
TvLDH        QATFTKEGKTQKVVDVLDHDYVFDTFFKKIGHRAWDILEHRGFTSAASPTKAAIQHMKAWLFGTA 
TvLDH_model  QATFTKEGKTQKVVDVLDHDYVFDTFFKKIGHRAWDILEHRGFTSAASPTKAAIQHMKAWLFGTA 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos           270       280       290       300       310       320
TvLDH        PGEVLSMGIPVPEGNPYGIKPGVVFSFPCNVDKEGKIHVVEGFKVNDWLREKLDFTEKDLFHEKE 
TvLDH_model  PGEVLSMGIPVPEGNPYGIKPGVVFSFPCNVDKEGKIHVVEGFKVNDWLREKLDFTEKDLFHEKE 
1emd         ----------------------------------------------------------------- 
 _consrvd


 _aln.pos      330
TvLDH        IALNHLAQGG/.. 
TvLDH_model  IALNHLAQGG/-- 
1emd         ----------/.. 
 _consrvd

比对好后我们在目标序列与带有配体的序列后加"/",表示序列结束,然后加"."表示匹配配体,该示例后面加了两个"."表示匹配两个配体。

最后进行建模:

from modeller import *
from modeller.automodel import *

class MyModel(automodel):
    def special_restraints(self, aln):
        rsr = self.restraints
        for ids in (('NH1:161:A', 'O1A:336:B'),
                    ('NH2:161:A', 'O1B:336:B'),
                    ('NE2:186:A', 'O2:336:B')):
            atoms = [self.atoms[i] for i in ids]
            rsr.add(forms.upper_bound(group=physical.upper_distance,
                                      feature=features.distance(*atoms),
                                      mean=3.5, stdev=0.1))

env = environ()
env.io.hetatm = True
a = MyModel(env, alnfile='TvLDH-1emd_bs.ali',
              knowns=('TvLDH_model','1emd'), sequence='TvLDH')
a.starting_model = 1
a.ending_model = 5
a.make()

File: ligand/model-multiple-hetero.py

建模脚本需要注意的是``env.io.hetatm = True`表示非标准”残基"也可以被读取,这样才能读取到配体。
分享到这里其实同源建模用的最普遍的内容已经讲完,接下来准备和大家一起学习迭代建模和建模更深的一些常用函数。这一系列就将先告一段落,等到时候研究rosetta的时候或者看到关于建模的一些好文章再继续撰写该类型。
往期教程:
#分子模拟#同源建模从入门到精通(一)
#分子模拟#同源建模从入门到精通(二)
#分子模拟#同源建模从入门到精通(三)

上一篇下一篇

猜你喜欢

热点阅读