利用MODLLER同源建模的基本步骤

最近开始系统的学习一下同源建模的过程,前段时间学习的时候简直就是一头雾水，不知道到底这玩艺儿怎么弄，现在终于理解了一点皮毛的东西，或许连皮毛也不算。

 

前期的准备工作（这里介绍并没有结构的蛋白序列的从头建立结构的过程）

首先：把你需要建立结构的蛋白序列写成modeller能够识别的文件格式

>P1;TvLDH

sequence:::::::::

MSEAAHVLITGAAGQIGYILSHWIASGELYGDRQVYLHLLDIPPAMNRLTALTMELEDCAFPHLAGFVATTDPKAAFKDIDCAFLVASMPLKPGQVRADLISSNSVIFKNTGEYLSKWAKPSVKVLVIGNPDNTNCEIAMLHAKNLKPENFSSLSMLDQNRAYYEVASKLGVDVKDVHDIIVWGNHGESMVADLTQATFTKEGKTQKVVDVLDHDYVFDTFFKKIGHRAWDILEHRGFTSAASPTKAAIQHMKAWLFGTAPGEVLSMGIPVPEGNPYGIKPGVVFSFPCNVDKEGKIHVVEGFKVNDWLREKLDFTEKDLFHEKEIALNHLAQGG*

第二行后面的字段基本没什么作用,只要把第一个字段sequence写上即可,注意，这一行一共10个字段用9个分号隔开。

第三行为将你要建立结构的序列的氨基酸写成单字母的形式，最后的 " * “ 符号表示序列的结束.

 

其次：将PDB库中的序列准备好（modeller会经常更新PDB库中的序列,这些序列会经常更新,你可以从modeller的官网上下载）下面是将PDB库中的序列转化为二进制的文件格式,可以加快搜索具有高度相似的序列的速度. 转换的脚本如下所示：

from modeller import *

log.verbose()

env = environ()

sdb = sequence_db(env)

sdb.convert(seq_database_file='pdb_95.pir', seq_database_format='PIR',  chains_list='ALL', minmax_db_seq_len=(30, 4000),  clean_sequences=True, outfile='pdb_95.hdf5')

 

第四行 sequence_db()对象用于创建一个存放大量的蛋白序列的库对象。

第五行调用对象的convert()方法，以PIR格式的方式读入pdb_95.pir,去除残基序列的数目小于30并且大于4000的序列，用minmx_db_seq_len参数指出，去除不标准的残基用clean_sequences指出，最终文件以二进制的形式写到pdb_95.hdf5文件中。

 

第三：利用上一步中生成的二进制文件进行搜索，找到和未知序列相似的结构序列。

from modeller import *

log.verbose()

env = environ()

# Read in the sequence database in binary format

sdb = sequence_db(env, seq_database_file='pdb_95.hdf5', seq_database_format='BINARY', chains_list='ALL')

# Read in the target sequence in PIR alignment format

aln = alignment(env)

aln.append(file='TvLDH.ali', alignment_format='PIR', align_codes='ALL')

# Convert the input sequence "alignment" into profile format

prf = aln.to_profile()

# Scan sequence database to pick up homologous sequences

prf.build(sdb, matrix_offset=-450, rr_file='${LIB}/blosum62.sim.mat',  gap_penalties_1d=(-500, -50), n_prof_iterations=1, check_profile=False, max_aln_evalue=0.01)

# Write out the profile in text format

prf.write(file='build_profile.prf', profile_format='TEXT')

# Convert the profile back to alignment format

aln = prf.to_alignment()

#- Write out a PIR alignment file

aln.write(file='build_profile.ali', alignment_format='PIR‘）

*****在build_profile.prf文件中，最重要的列为第2,10,11，12列，第二列指出了PDB库中的晶体结构的代码，第11列指出了待建结构与晶体结构的序列相似度的百分数，第12列指明了叠加之后的e值。

第四：从build_profile.prf这个输出文件中找到和待建序列的相似度高的模版，评价模版的结构和序列的相似性。

from modeller import *

env = environ()

env.io.atom_files_directory = ['.', '../atom_files']

# Make a simple 1:1 alignment of 7 template structures

aln = alignment(env)for (pdb, chain) in (('1b8p', 'A'), ('1y7t', 'A'), ('1civ', 'A'), ('5mdh', 'A'), ('7mdh', 'A'), ('3d5t', 'A'),   ('1smk', 'A')):   

 m = model(env, file=pdb, model_segment=('FIRST:'+chain, 'LAST:'+chain))  

 aln.append_model(m, atom_files=pdb, align_codes=pdb+chain)

# Sequence alignment

aln.malign()

# Structure alignment

aln.malign3d()

# Report details of the sequence/structure alignment

aln.compare_structures()

aln.id_table(matrix_file='family.mat')

env.dendrogram(matrix_file='family.mat', cluster_cut=-1.0)

log文件中的输出结果如下面的图所示：

图片中的第一部分展示了序列相似性的比较：

对角线中的内容为每个结构的残基数目，上三角为相同的残基的数目，下三角为两个结构的序列相似度的百分数。

我们选择了1y7t作为模版结构，因为这个结构的结晶分辨率为1.6埃以及它与待建结构的序列的相似度为45%。

第五：将待建结构与模版结构进行align, MODELLER利用这些align的结果来提取出构建结构时必要的限制。

最好的方法是利用align2d()的方法，尽管align2d()方法是基于一种动态的算法，但是它不同于标准的序列-序列的align

方法，因为这种方法在进行align的时候会考虑到模版的结构信息。

from modeller import *

env = environ()

env.io.atom_files_directory = ['.', '../atom_files']

aln = alignment(env)

# Read in the 1y7t template structure and add to alignment

mdl = model(env, file='1y7t', model_segment=('FIRST:A','LAST:A'))

aln.append_model(mdl, align_codes='1y7tA', atom_files='1y7t')

# Add in the TvLDH sequence

aln.append(file='../step1_search/TvLDH.ali', align_codes='TvLDH')

# Sequence/structure alignment

aln.align2d(max_gap_length=40)

# Write out resulting alignment in both PIR and PAP formats

aln.write(file='TvLDH-1y7tA.ali', alignment_format='PIR')

aln.write(file='TvLDH-1y7tA.pap', alignment_format='PAP')

 

第六：创建比较模型

可以使用MODELLER中的automodel类，只需要告诉automodel alignment文件，使用哪个代码以及要建立多少个模型。

from modeller import *

from modeller.automodel import *

env = environ()

env.io.atom_files_directory = ['.', '../atom_files']

a = automodel(env,alnfile='../step3_align/TvLDH1y7tA.ali',  knowns='1y7tA',sequence='TvLDH',assess_methods=assess.DOPE)

a.starting_model = 1

a.ending_model = 5

a.make()

 

第七：评价所建的模型

从上面的步骤中找到最好的模型，评价的方法有很多种，最简单的方法是molecule PDF (molpdf)方法，标准的MODELLER评分函数。这个评分结果会在每个PDB开头的REMARK中加入。其他的评价方法包括DOPE和GA341.

从上面产生的log文件中可以发现，第一个PDB文件模型是最好的，因为DOPE的值最小。

 

第八：Fit models into cryo-EM maps

通过将TvLDH的 cryo-EM数据导入来使产生的蛋白结构与这些数据相吻合，通过使用MODELLER中的Mod-EM方法。

from modeller import *

log.verbose()

env = environ()

struct='../step4_model/TvLDH.B99990001.pdb'

map='TvLDH.10A.mrc' 

resolution=10.0 

box_size=48

apix=1.88

x=-26.742; y=-9.5205; z=-10.375 #origin

steps=20

# Read in cryo-EM density map

den = density(env, file=map, em_density_format='MRC',  voxel_size=apix, resolution=resolution, em_map_size=box_size,  density_type='GAUSS', px=x,py=y,pz=z)

# Fit the PDB file into the map by MC simulated annealing

den.grid_search(em_density_format='MRC', num_structures=1,  em_pdb_name=struct, chains_num=[1],                start_type='CENTER', number_of_steps=steps,   angular_step_size=30., temperature=100.,                best_docked_models=1, translate_type='RANDOM',  em_fit_output_file='modem.log')