Modeller 学习记录（二）

lrf1980

前面讲完了安装及运行，那基本上就到了实际操作，建立自己要做的蛋白质模型了。开始建模型前我做的准备工作：
1. 到salilab网站的tutorial去下载basic modeling example， advanced modeling example。因为我不会python，所以我只能依葫芦画瓢用这些例子里的脚本来做自己的工作。
2. 到uniprot上去找自己要做的蛋白质序列，把需要建模的那段序列的ali文件建好。怎么建？很简单，以葫芦画瓢，下面的例子是example里的ali
>P1;TvLDH
sequence:TvLDH:::::::0.00: 0.00
MSEAAHVLITGAAGQIGYILSHWIASGELYGDRQVYLHLLDIPPAMNRLTALTMELEDCAFPHLAGFVATTDPKA
AFKDIDCAFLVASMPLKPGQVRADLISSNSVIFKNTGEYLSKWAKPSVKVLVIGNPDNTNCEIAMLHAKNLKPEN
FSSLSMLDQNRAYYEVASKLGVDVKDVHDIIVWGNHGESMVADLTQATFTKEGKTQKVVDVLDHDYVFDTFFKKI
GHRAWDILEHRGFTSAASPTKAAIQHMKAWLFGTAPGEVLSMGIPVPEGNPYGIKPGVVFSFPCNVDKEGKIHVV
EGFKVNDWLREKLDFTEKDLFHEKEIALNHLAQGG*

用你的蛋白质序列替换TvLDH的序列，当然也需要用你的蛋白质名称替换TvLDH。然后存成ali后缀格式文件就好。
3. 准备好最新的pdb数据库文件，pir和bin。前面讲过了，这两个文件怎么生成。
4. 建一个以蛋白质文件命名的目录，将所有相关的文件都放这个目录，然后所有的工作都在这个目录里进行。（*.ali, *.pir, *.bin, build_profile.py,及从basic modeling example目录里拷贝过来的各种脚本）

建模开始：
1. build_profile.py的修改，下面贴出原始的build_profile.py文件，从basic example里来的。改动几个地方，就可以用这个来查找相似蛋白质信息。

from modeller import *

log.verbose()
env = environ()

#-- Prepare the input files

#-- Read in the sequence database
sdb = sequence_db(env)
sdb.read(seq_database_file='pdb_95.pir', seq_database_format='PIR',
         chains_list='ALL', minmax_db_seq_len=(30, 4000), clean_sequences=True)
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
我用了pdball.pir来替换这个pdb_95.pir
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
#-- Write the sequence database in binary form
sdb.write(seq_database_file='pdb_95.bin', seq_database_format='BINARY',
          chains_list='ALL')
#-- Now, read in the binary database
sdb.read(seq_database_file='pdb_95.bin', seq_database_format='BINARY',
         chains_list='ALL')
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
同上，我用了pdball.bin来替换这个pdb_95.bin
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

#-- Read in the target sequence/alignment
aln = alignment(env)
aln.append(file='TvLDH.ali', alignment_format='PIR', align_codes='ALL')
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
用自己的蛋白质ali文件替代TvLDH.ali
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
#-- 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 the alignment file
aln.write(file='build_profile.ali', alignment_format='PIR')


最后在命令行下运行mod9.14 build_profile.py
完成之后会在目录里生成三个文件：
build_profile.ali
build_profile.log
build_profile.prf

我只用log文件去抓蛋白质相似性信息。其它两个文件不怎么用。抓出下面这部分的信息就可以，主要的目的是为了挑选出相似性高的蛋白质来做进一步的分析。
Dynamically allocated memory at     amaxprofile [B,KiB,MiB]:      7496268    7320.574     7.149
> 1a5z                       1     190    6900     312     335   28.14    0.85E-08       2   164    75   242    63   229
> 1b8pA                      1     951   29400     327     335   42.50     0.0           3   316     7   331     6   325
> 1bdmA                      1    1058   30700     318     335   44.52     0.0           4   309     1   325     1   310
> 1t2dA                      1    1717    5200     315     335   24.70    0.67E-04       5   238     5   256     4   250
> 1civA                      1    1784   23950     374     335   34.66     0.0           6   325     6   334    33   358
> 2cmd                       1    1858    5900     312     335   27.24    0.17E-05       7   289     7   320     3   303
> 1o6zA                      1    2060    5800     303     335   26.32    0.27E-05       8   278     7   320     3   287
> 1ur5A                      1    3954    4500     299     335   30.67    0.25E-02       9   158    13   191     9   171
> 1guzA                      1    3955    7100     305     335   25.27    0.29E-08      10   265    13   301     8   280
> 1gv0A                      1    3962    5350     301     335   25.53    0.29E-04      11   274    13   323     8   289
> 1hyeA                      1    4405    6800     307     335   29.28    0.14E-07      12   173     7   191     3   183
> 1i0zA                      1    4433    5650     332     335   24.64    0.67E-05      13   207    85   300    94   304
> 1i10A                      1    4435    5600     331     335   26.21    0.87E-05      14   196    85   295    93   298
> 1ldnA                      1    5824    5000     316     335   25.76    0.19E-03      15   214    78   298    73   301
> 6ldh                       1    5831    4600     329     335   23.08    0.17E-02      16   244    47   301    56   302
> 2ldx                       1    5833    5400     331     335   25.83    0.25E-04      17   227    66   306    67   306
> 5ldh                       1    5835    5800     333     335   25.59    0.30E-05      18   207    85   300    94   304
> 9ldtA                      1    5836    6000     331     335   25.94    0.11E-05      19   207    85   301    93   304
> 1llc                       1    5896    5000     321     335   25.82    0.20E-03      20   164    64   239    53   234
> 1lldA                      1    5897    6650     313     335   30.67    0.32E-07      21   216    13   242     9   233
> 5mdhA                      1    6084   32900     333     335   44.41     0.0          22   328     2   332     1   331
> 7mdhA                      1    6085   22600     351     335   33.74     0.0          23   325     6   334    14   339
> 1mldA                      1    6188    6100     313     335   25.93    0.58E-06      24   183     5   198     1   189
> 1oc4A                      1    6736    5450     315     335   27.87    0.18E-04      25   174     5   191     4   186
> 1ojuA                      1    6799    5700     294     335   27.98    0.44E-05      26   218    78   320    68   285
> 1pzgA                      1    7253    6350     327     335   30.00    0.16E-06      27   114    74   191    71   190
> 1smkA                      1    7850    8750     313     335   33.85     0.0          28   188     7   202     4   198
> 1sovA                      1    7980    4700     316     335   26.59    0.94E-03      29   160    81   256    76   248
> 1y6jA                      1   10912    5750     289     335   32.73    0.33E-05      30   109    77   191    58   167


把这部分信息整理，选高相似度及E值尽可能小的蛋白质来做下一步的compare分析。