lrf1980 发表于 2014-11-15 22:04:46

Modeller 学习记录(五)

本帖最后由 lrf1980 于 2014-11-15 22:07 编辑

接上面的内容,这里将涉及多模版建模,单模版包含配体建模,及多模版包含配体建模。

如果对前面四篇都能够完全走下来,现在就可以开始来着手做多模版建模了,多模版建模跟单模版建模差不多,只是用到的align脚本不一样。
对多模版建模,当我们选定了需要的模版,我们就可以用salign.py这个脚本。advanced example里的salign.py脚本内容:
# Illustrates the SALIGN multiple structure/sequence alignment

from modeller import *

log.verbose()
env = environ()
env.io.atom_files_directory = './:../atom_files/'

aln = alignment(env)
for (code, chain) in (('2mdh', 'A'), ('1bdm', 'A'), ('1b8p', 'A')):
    mdl = model(env, file=code, model_segment=('FIRST:'+chain, 'LAST:'+chain))
    aln.append_model(mdl, atom_files=code, align_codes=code+chain)
---------------------------------------
基本上需要做的也就是将2mdh,1bdm,1b8p等替换成选定的模版code及相应的chain code就好
---------------------------------------
for (weights, write_fit, whole) in (((1., 0., 0., 0., 1., 0.), False, True),
                                    ((1., 0.5, 1., 1., 1., 0.), False, True),
                                    ((1., 1., 1., 1., 1., 0.), True, False)):
    aln.salign(rms_cutoff=3.5, normalize_pp_scores=False,
               rr_file='$(LIB)/as1.sim.mat', overhang=30,
               gap_penalties_1d=(-450, -50),
               gap_penalties_3d=(0, 3), gap_gap_score=0, gap_residue_score=0,
               dendrogram_file='fm00495.tree',
               alignment_type='tree', # If 'progresive', the tree is not
                                    # computed and all structues will be
                                    # aligned sequentially to the first
               feature_weights=weights, # For a multiple sequence alignment only
                                        # the first feature needs to be non-zero
               improve_alignment=True, fit=True, write_fit=write_fit,
               write_whole_pdb=whole, output='ALIGNMENT QUALITY')
--------------------------------------
对于这设置,我一直想搞懂都代表什么,无奈能力有限,还没搞清楚,有谁了解的,分享一下吧
--------------------------------------
aln.write(file='fm00495.pap', alignment_format='PAP')
aln.write(file='fm00495.ali', alignment_format='PIR')

aln.salign(rms_cutoff=1.0, normalize_pp_scores=False,
         rr_file='$(LIB)/as1.sim.mat', overhang=30,
         gap_penalties_1d=(-450, -50), gap_penalties_3d=(0, 3),
         gap_gap_score=0, gap_residue_score=0, dendrogram_file='1is3A.tree',
         alignment_type='progressive', feature_weights=*6,
         improve_alignment=False, fit=False, write_fit=True,
         write_whole_pdb=False, output='QUALITY')

这个salig.py运行完成,就可以来用下一个脚本align2d_mult.py来生成align的ali文件。

from modeller import *

log.verbose()
env = environ()

env.libs.topology.read(file='$(LIB)/top_heav.lib')

# Read aligned structure(s):
aln = alignment(env)
aln.append(file='fm00495.ali', align_codes='all')
aln_block = len(aln)

# Read aligned sequence(s):
aln.append(file='TvLDH.ali', align_codes='TvLDH')

# Structure sensitive variable gap penalty sequence-sequence alignment:
aln.salign(output='', max_gap_length=20,
         gap_function=True,   # to use structure-dependent gap penalty
         alignment_type='PAIRWISE', align_block=aln_block,
         feature_weights=(1., 0., 0., 0., 0., 0.), overhang=0,
         gap_penalties_1d=(-450, 0),
         gap_penalties_2d=(0.35, 1.2, 0.9, 1.2, 0.6, 8.6, 1.2, 0., 0.),
         similarity_flag=True)
------------------------------------
上面这段设置,也不不明白具体含义,有了解的同志来帮忙解读一下吧
------------------------------------
aln.write(file='TvLDH-mult.ali', alignment_format='PIR')


最后运行model_mult.py就可以建模了。

from modeller import *
from modeller.automodel import *

env = environ()
a = automodel(env, alnfile='TvLDH-mult.ali',
            knowns=('1bdmA','2mdhA','1b8pA'), sequence='TvLDH')
a.starting_model = 1
a.ending_model = 5
a.make()


在这里补充一点,可以在这个脚本里添加计算DOPE的语句,这样,生成的log文件里可以直接看到所有建好的模型的DOPE分数。就不用再去用evaluate_model.py去跑一遍
将下面这段代码夹在a.starting_model=1 与knowns=('1bdmA',........)之间就可以

assess_methods=(assess.DOPE,
                              #soap_protein_od.Scorer(),
                              assess.GA341))


对于单模版包含配体的建模,基本上步骤跟单模版建模一样。当我们选好模版运行好align2d.py这个脚本后,我们会得到一个目标蛋白和模版蛋白序列align文件,结尾为.ali, 在simple example里就是TvLDH_1bmdA.ali. 1bmd就是我们的模版。
为了说明问题,我摘录部分该文件的内容如下:
>;P1;1bdmA
structureX:1bdm.pdb:   0 :A:+318 :A:undefined:undefined:-1.00:-1.00
VEGLEINEFARKRMEITAQELLDEMEQVKAL--GLI*

>;P1;TvLDH
sequence:TvLDH:   : :   : ::: 0.00: 0.00
VEGFKVNDWLREKLDFTEKDLFHEKEIALNHLAQGG*


假设我们的模版1bdmA中包含一个配体,我们需要在建模的时候把这个配体包含到我们建的模型中,那么就在1bmdA序列的结尾将序列由
VEGLEINEFARKRMEITAQELLDEMEQVKAL--GLI* ->改成 VEGLEINEFARKRMEITAQELLDEMEQVKAL--GLI.*
目标蛋白序列也做相应的更改
VEGFKVNDWLREKLDFTEKDLFHEKEIALNHLAQGG* -> 改成 VEGFKVNDWLREKLDFTEKDLFHEKEIALNHLAQGG.*
这里可能还需要把318改成319。 因为我们添加了“.”到序列里。这个是否需要改,运行一下model_single.py就知道了。另外model_single.py里面需要添加下面的代码:
env = environ()
env.io.hetatm = True #这部分是添加进来的,告诉程序将模版里面的hetatm复制到模型中去
a = automodel(env, alnfile='HAhR-3H82.ali',
基本上,包含配体的单模版建模就是这么一个过程。最后当模型建好后,成功的标准就是你是否能够在你建的模型里面找到改配体。有该配体,就说明成功实现包含配体的建模,如果没有,那就是说明有地方没有弄好,需要仔细检查。
为了保险起见,个人感觉,需要稍微对模版蛋白质文件做简单的处理,比如说,去掉结构中的溶剂分子,将不需要的chain remove掉。尽量做到蛋白质文件里只包含用到的序列及配体,这样成功的机会会大很多。所以包含配体建模的语法就是:
1. 添加env.io.hetatm = True
2. 在align的序列末尾加".“,并以*作为结尾

多模版包含配体的做法,这里只针对多模版中只有一个模版包含配体的情况。操作步骤基本一致。先做完salign.py 及align2d_mult.py,然后对生成的ali文件修改。在包含配体的模版序列及目标蛋白序列后添加".",其它模版的序列结尾用"-"替代"."。基本上就可以了。
按照modeller的说明,如果是包含配体及水分子应该用"/.."放在序列结尾,目前还没有尝试过这个,所以没有实际应用的经验,在这里不做说明。大家可以试试。

基本上这就可以完成建模了。
PS.有任何问题及建议,请站内信件,勿在群中qq联系,qq被盗。自己已经上不了了。

woniuzhixing 发表于 2014-11-22 15:27:38

lrf1980 发表于 2014-11-21 07:44
SOAP 打分开启方法

从下面网站上下载SOAP相关文件放到目录 modeller-9.14/modlib下


我下载的是basic-example,modeller9.14,这两个文件夹均放在了C盘中,目录级别是一样的,并没有把谁放在谁里面,接着将soap_protein_od.hdf5文件下载下来,属性显示是180M大小,然后该文件放在了modeller-9.14/modlib/目录中,结果显示如下:
>> ENERGY; Differences between the model's features and restraints:
Number of all residues in MODEL                   :      335
Number of all, selected real atoms                :   2605    2605
Number of all, selected pseudo atoms            :      0       0
Number of all static, selected restraints         :    28624   28624
COVALENT_CYS                                    :      F
NONBONDED_SEL_ATOMS                               :      1
Number of non-bonded pairs (excluding 1-2,1-3,1-4):   523644
Dynamic pairs routine                           : 1, NATM x NATM double loop
Atomic shift for contacts update (UPDATE_DYNAMIC) :    0.390
LENNARD_JONES_SWITCH                              :    6.500   7.500
COULOMB_JONES_SWITCH                              :    6.500   7.500
RESIDUE_SPAN_RANGE                              :      1    9999
NLOGN_USE                                       :       15
CONTACT_SHELL                                     :   15.000
DYNAMIC_PAIRS,_SPHERE,_COULOMB,_LENNARD,_MODELLER :      T       F       F       F       T
SPHERE_STDV                                       :    0.050
RADII_FACTOR                                    :    0.820
Current energy                                    :      -38043.4922




<< end of ENERGY.
DOPE score               : -38043.492188
>> Model assessment by SOAP-Protein-OD score
hdf5err____E> unable to open file


Traceback (most recent call last):
File "C:\basic-example\model-single.py", line 13, in <module>
    a.make()
File "C:\Modeller9.14\modlib\modeller\automodel\automodel.py", line 121, in make
    self.multiple_models(atmsel)
File "C:\Modeller9.14\modlib\modeller\automodel\automodel.py", line 225, in multiple_models
    self.outputs.append(self.single_model(atmsel, num))
File "C:\Modeller9.14\modlib\modeller\automodel\automodel.py", line 322, in single_model
    self.model_analysis(atmsel, filename, out, num)
File "C:\Modeller9.14\modlib\modeller\automodel\automodel.py", line 369, in model_analysis
    assess_keys = self.assess(atmsel, self.assess_methods, out)
File "C:\Modeller9.14\modlib\modeller\automodel\automodel.py", line 417, in assess
    (key,value) = method(atmsel)
File "C:\Modeller9.14\modlib\modeller\terms.py", line 60, in __call__
    return (self.name, atmsel.assess(self))
File "C:\Modeller9.14\modlib\modeller\selection.py", line 655, in assess
    molpdf, terms = assessor._assess(self, output=output, **vars)
File "C:\Modeller9.14\modlib\modeller\terms.py", line 56, in _assess
    return atmsel.energy(edat=self._get_energy_data_cached(),
File "C:\Modeller9.14\modlib\modeller\terms.py", line 68, in _get_energy_data_cached
    self._edat = self._get_energy_data_all()
File "C:\Modeller9.14\modlib\modeller\terms.py", line 73, in _get_energy_data_all
    edat.energy_terms.append(self)
File "C:\Modeller9.14\modlib\modeller\util\modlist.py", line 155, in append
    self.insert(len(self), obj)
File "C:\Modeller9.14\modlib\modeller\util\modlist.py", line 167, in insert
    self._insfunc(indx, obj)
File "C:\Modeller9.14\modlib\modeller\terms.py", line 38, in _insfunc
    obj._add_term(self.__edat(), indx)
File "C:\Modeller9.14\modlib\modeller\soap_protein_od.py", line 21, in _add_term
    self.__library)
ModellerError: hdf5err____E> unable to open file
显示不能打开文件,跪求照葫芦画瓢为啥还是错??python是2.7版本的。。
>>>

zhuimeng0812 发表于 2015-5-21 17:03:33

lrf1980 发表于 2015-5-20 21:30
要包含配体的话,好像不能用*.fit.pdb做模版,你看看*_fit.pdb,这里面是没有配体信息的,你还需要改成 ...

你的意思是 只保留一个蛋白中有配体信息,其他的模板 把配体去掉!
我再向你请教一个问题,modeller 同源模建的时候,提供的模板是不是必须要有链的信息。比如:A、B、C......链等,而把这些表示链的信息A、B、C...... 去掉,modeller 就不识别了?对吗? 如下面的B 去掉:

ATOM   4731N   PRO B   9   111.97133.33097.8561.00 92.26         N
ATOM   4732CAPRO B   9   112.16034.79897.7601.00 92.31         C
ATOM   4733C   PRO B   9   111.12735.40796.8161.00 93.29         C
ATOM   4734O   PRO B   9   110.96634.95695.6771.00 92.89         O
ATOM   4735CBPRO B   9   113.57335.03697.2371.00 91.79         C
ATOM   4736CGPRO B   9   114.27033.72297.6241.00 93.83         C
ATOM   4737CDPRO B   9   113.19632.63597.4401.00 91.53         C


报的错误是:
modeller.ModellerError: rdpdb___303E> No atoms were read from the specified input PDB file, since the starting residue number and/or chain id in MODEL_SEGMENT (or the alignment file header) was not found; requested starting position: residue number " 6", chain " C"; atom file name:1zoy.pdb


让我奇怪的是:swiss-model 却是可以构建的,swiss-model 也是用的modeller 的模块,只是更便于操作而已,为什么 swiss-model可以,而modeller却不识别呢?

请你指导一下!谢谢!

lrf1980 发表于 2014-11-21 07:44:17

woniuzhixing 发表于 2014-11-18 20:56
不是啊,你已经做了呀,单分子建模的,assess_methods=(assess.DOPE,
                              #soap ...

SOAP 打分开启方法

从下面网站上下载SOAP相关文件放到目录 modeller-9.14/modlib下
http://salilab.org/SOAP/

例如我做的事protein 的soap评估,那么就在modlib下面有 这个文件:modeller-9.14/modlib//soap_protein_od.hdf5

以model-simple.py 这个脚本为例,开启SOAP
from modeller import *
from modeller.automodel import *
#from modeller import soap_protein_od

env = environ()
a = automodel(env, alnfile='TvLDH-1bdmA.ali',
            knowns='1bdmA', sequence='TvLDH',
            assess_methods=(assess.DOPE,
                              #soap_protein_od.Scorer(),
                              assess.GA341))
a.starting_model = 1
a.ending_model = 5
a.make()

将脚本中被注释掉的两行开启 (删除#就可以),然后运行就可以了。

我运行成功的例子结果

>> Summary of successfully produced models:
Filename                        molpdf   DOPE score    GA341 score SOAP-Protein-OD score
--------------------------------------------------------------------------------------------
***.B99990001.pdb             656.49341    -9962.69434      0.97499   -63212.68750
***.B99990002.pdb             631.16901    -9959.50488      0.73881   -63845.61719

woniuzhixing 发表于 2014-11-18 09:16:59

你知道SOAP对模型打分评价吗怎么操作吗?教程说下载SOAP-Protein-OD library file,但是下载下来依旧不能打分评估。求解答、、、

lrf1980 发表于 2014-11-18 20:14:21

woniuzhixing 发表于 2014-11-18 09:16
你知道SOAP对模型打分评价吗怎么操作吗?教程说下载SOAP-Protein-OD library file,但是下载下来依旧不能打 ...

我还不知道啊,没有学到这里。问问川大-灰太狼看看。

woniuzhixing 发表于 2014-11-18 20:56:01

不是啊,你已经做了呀,单分子建模的,assess_methods=(assess.DOPE,
                              #soap_protein_od.Scorer(),
                              assess.GA341)),中间打着#就是啊,不过#的作用是不进行该操作,教程讲去掉#,下载SOAP-Protein-OD library file文件即可,但是我下载了,结果显示modeller无法打开该文件,求解啊?

lrf1980 发表于 2014-11-18 22:03:17

woniuzhixing 发表于 2014-11-18 20:56
不是啊,你已经做了呀,单分子建模的,assess_methods=(assess.DOPE,
                              #soap ...

还是没有做嘛,因为我是直接拷贝里面的脚本,我改天试试,要是能成功,我再分享。如何。

woniuzhixing 发表于 2014-11-19 08:17:30

好呀,大家都试试,做出互相分享。。

puzhongji 发表于 2014-11-19 22:00:40

多模板建模:
The structure of the Malate Dehydrogenase 1bdm has been clustered in the DBAli database within the family fm00495 of 4 members .具体是怎么做的?用DBAli工具??

lrf1980 发表于 2014-11-19 23:31:00

woniuzhixing 发表于 2014-11-18 09:16
你知道SOAP对模型打分评价吗怎么操作吗?教程说下载SOAP-Protein-OD library file,但是下载下来依旧不能打 ...

这个应该可能是你要的脚本

# Example for: selection.assess(), soap_protein_od.Scorer()

from modeller import *
from modeller.scripts import complete_pdb
from modeller import soap_protein_od

env = environ()
env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

# Set up SOAP-Protein-OD scoring (note: if assessing multiple models, it is
# best to create 'sp' just once and keep it around, since reading in the
# potential from disk can take a long time).
sp = soap_protein_od.Scorer()

# Read a model previously generated by Modeller's automodel class
mdl = complete_pdb(env, '../atom_files/1fdx.B99990001.pdb')

# Select all atoms in the first chain
atmsel = selection(mdl.chains)

# Assess with the above Scorer
try:
    score = atmsel.assess(sp)
except ModellerError:
    print("The SOAP-Protein-OD library file is not included with MODELLER.")
    print("Please get it from http://salilab.org/SOAP/.")
页: [1] 2 3
查看完整版本: Modeller 学习记录(五)