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被盗。自己已经上不了了。
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版本的。。
>>> 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却不识别呢？
请你指导一下！谢谢！ 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
你知道SOAP对模型打分评价吗怎么操作吗？教程说下载SOAP-Protein-OD library file，但是下载下来依旧不能打分评估。求解答、、、 woniuzhixing 发表于 2014-11-18 09:16
你知道SOAP对模型打分评价吗怎么操作吗？教程说下载SOAP-Protein-OD library file，但是下载下来依旧不能打 ...
我还不知道啊，没有学到这里。问问川大－灰太狼看看。 不是啊，你已经做了呀，单分子建模的，assess_methods=(assess.DOPE,
#soap_protein_od.Scorer(),
assess.GA341))，中间打着#就是啊，不过#的作用是不进行该操作，教程讲去掉#，下载SOAP-Protein-OD library file文件即可，但是我下载了，结果显示modeller无法打开该文件，求解啊？
woniuzhixing 发表于 2014-11-18 20:56
不是啊，你已经做了呀，单分子建模的，assess_methods=(assess.DOPE,
#soap ...
还是没有做嘛，因为我是直接拷贝里面的脚本，我改天试试，要是能成功，我再分享。如何。
好呀，大家都试试，做出互相分享。。 多模板建模：
The structure of the Malate Dehydrogenase 1bdm has been clustered in the DBAli database within the family fm00495 of 4 members .具体是怎么做的？用DBAli工具？？ 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/.")