Protein Structure Modeling with Modeller

Leave a Comment / Computation Biochemistry / By

Introduction

Modeller는 Andrej Sali 연구실에서 개발된 대표적인 상동 모델링(homology modeling) 도구이다. 본 실습에서는 해당 연구실에서 제공하는 튜토리얼을 통해 단백질 구조 모델링의 전반적인 과정을 학습한다.

Tutorial 요약

여기서 연습하는 과장은 아래와 같이 구성되어 있다.

  1. Search for structures related to TvLDH
  2. Selecting a template
  3. Aligning TvLDF with the template
  4. Model building
  5. Model evaluation
  6. Visualization

Tutorial 각 과정에 대한 자세한 설명은 아래 링크를 참고한다.

https://salilab.org/modeller/tutorial/basic.html

여기서는 편의상 해당 tutorial 의 예제에 기반하여, directory 구조만 약간 변형한 예제를 이용한다.

아래 파일을 다운 받고, tutorial을 진행할 directory에서 압축을 풀고 tutorial을 진행한다.

basic-example-mod.tarDownload

1. Search for structures related to TvLDH

특정 단백질의 서열을 이용해 3차원 구조를 모델링하기 위해서는, 해당 서열과 유사한 구조를 가진 단백질을 찾아야 한다. 이를 위해서는 구조가 밝혀진 단백질로 구성된 데이터베이스가 필요하다. 대표적인 예로는 Protein Data Bank가 있다. 구조 데이터베이스에서 모델링하고자 하는 서열과 유사한 단백질을 찾기 위해서는 서열 비교를 수행해야 한다. 본 예제에서는 TvLDH의 서열을 구조 데이터베이스 내 단백질의 서열과 비교하는 과정이 필요하다.

이를 위해 본 예제에서는 서열 유사도가 95% 이하인 단백질만을 포함한 데이터베이스(database/pdb_95.pir)를 활용한다. 즉, 서열 유사도가 95% 이상인 단백질들 중에서는 하나를 대표 구조로 선택하여 데이터베이스를 구성한 것이다.

Modeller 는 서열 검색 기능을 제공한다. 아래는 Modeller를 이용해 서열 검색하는 Python 코드이다. (예제 파일: 1_ build_profile.py)

from modeller import *

log.verbose()
env = Environ()

#-- Prepare the input files

#-- Read in the sequence database
sdb = SequenceDB(env)
sdb.read(seq_database_file='database/pdb_95.pir', seq_database_format='PIR',
         chains_list='ALL', minmax_db_seq_len=(30, 4000), clean_sequences=True)

#-- Write the sequence database in binary form
sdb.write(seq_database_file='database/pdb_95.bin', seq_database_format='BINARY',
          chains_list='ALL')

#-- Now, read in the binary database
sdb.read(seq_database_file='database/pdb_95.bin', seq_database_format='BINARY',
         chains_list='ALL')

#-- Read in the target sequence/alignment
aln = Alignment(env)
aln.append(file='input/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='results_aln/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='results_aln/build_profile.ali', alignment_format='PIR')

2. Select a template

데이터베이스 검색 결과에서 가장 정렬이 잘 된 서열을 찾는다. 그 기준은 구조에 정렬된 서열 부위의 길이와 정렬된 부분의 서열 유사성이다. 서열의 대부분이 정렬되고 높은 서열 유사성을 가진 단백질을 템플릿(template)으로 선택한다

이를 위해 Modeller의 기능을 일부 활용할 수 있다. 예제에서는 2_compare.py를 이용해 정렬된 서열의 관계를 서로 비교할 수도 있다. 

from modeller import *

env = Environ()
aln = Alignment(env)
for (pdb, chain) in (('1b8p', 'A'), ('1bdm', 'A'), ('1civ', 'A'),
                     ('5mdh', 'A'), ('7mdh', 'A'), ('1smk', 'A')):
    m = Model(env, file='pdb/'+pdb, model_segment=('FIRST:'+chain, 'LAST:'+chain))
    aln.append_model(m, atom_files='pdb/'+pdb, align_codes=pdb+chain)
aln.malign()
aln.malign3d()
aln.compare_structures()
aln.id_table(matrix_file='results_aln/family.mat')
env.dendrogram(matrix_file='results_aln/family.mat', cluster_cut=-1.0)

3. Aligning TvLDF with the template

Template을 선택한 후, 서열의 어떤 부분이 단백질 구조에 해당하는지를 확인하기 위해 추가적인 정렬을 수행한다. 이 과정은 예제의 3_align2d.py 명령을 통해 실행된다. 정렬이 완료되면, 정렬된 서열 정보가 담긴 *.ali 및 *.pap 파일이 생성된다.

from modeller import *

env = Environ()
aln = Alignment(env)
mdl = Model(env, file='pdb/1bdm', model_segment=('FIRST:A','LAST:A'))
aln.append_model(mdl, align_codes='1bdmA', atom_files='pdb/1bdm.pdb')
aln.append(file='input/TvLDH.ali', align_codes='TvLDH')
aln.align2d(max_gap_length=50)
aln.write(file='results_aln/TvLDH-1bdmA.ali', alignment_format='PIR')
aln.write(file='results_aln/TvLDH-1bdmA.pap', alignment_format='PAP')

4. Model building

서열을 템플릿에 정렬한 후, 해당 서열의 3차원 구조를 모델링한다. 이 과정에서는 템플릿의 아미노산을 대상 서열의 아미노산으로 치환하고, 구조적 에너지 계산을 통해 에너지가 낮은 구조를 탐색한다. 모델링에는 무작위적인 요소가 포함되므로, 양질의 구조를 얻기 위해 여러 번의 모델링을 수행하는 것이 일반적이다. 본 예제에서는 5회의 구조 모델링을 진행한다.

from modeller import *
from modeller.automodel import *

env = Environ()
a = AutoModel(env, alnfile='results_aln/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()

import os
os.system('mv TvLDH* results_models')

5. Model evaluation

모델링된 구조의 정확성을 평가하기 위해, 모델의 에너지 특성을 분석한다. 각 잔기(residue)의 안정성을 확인하기 위해 DOPE(Distance-scaled, Optimized Protein Energy) 스코어를 계산하며, 이를 템플릿 구조의 DOPE 스코어와 비교하여 모델의 품질을 판단한다.

from modeller import *
from modeller.scripts import complete_pdb

log.verbose()    # request verbose output
env = Environ()
env.libs.topology.read(file='$(LIB)/top_heav.lib') # read topology
env.libs.parameters.read(file='$(LIB)/par.lib') # read parameters

# read model file
mdl = complete_pdb(env, 'results_models/TvLDH.B99990001.pdb')

# Assess with DOPE:
s = Selection(mdl)   # all atom selection
s.assess_dope(output='ENERGY_PROFILE NO_REPORT', file='results_models/TvLDH.profile',
              normalize_profile=True, smoothing_window=15)

# read model file
mdl = complete_pdb(env, 'pdb/1bdm.pdb')


from modeller import *
from modeller.scripts import complete_pdb

log.verbose()    # request verbose output
env = Environ()
env.libs.topology.read(file='$(LIB)/top_heav.lib') # read topology
env.libs.parameters.read(file='$(LIB)/par.lib') # read parameters

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

# read model file
mdl = complete_pdb(env, 'pdb/1bdm.pdb', model_segment=('FIRST:A', 'LAST:A'))

s = Selection(mdl)
s.assess_dope(output='ENERGY_PROFILE NO_REPORT', file='results_models/1bdmA.profile',
              normalize_profile=True, smoothing_window=15)


import matplotlib.pyplot as plt
import modeller

def r_enumerate(seq):
    """Enumerate a sequence in reverse order"""
    # Note that we don't use reversed() since Python 2.3 doesn't have it
    num = len(seq) - 1
    while num >= 0:
        yield num, seq[num]
        num -= 1

def get_profile(profile_file, seq):
    """Read `profile_file` into a Python array, and add gaps corresponding to
       the alignment sequence `seq`."""
    # Read all non-comment and non-blank lines from the file:
    f = open(profile_file)
    vals = []
    for line in f:
        if not line.startswith('#') and len(line) > 10:
            spl = line.split()
            vals.append(float(spl[-1]))
    # Insert gaps into the profile corresponding to those in seq:
    for n, res in r_enumerate(seq.residues):
        for gap in range(res.get_leading_gaps()):
            vals.insert(n, None)
    # Add a gap at position '0', so that we effectively count from 1:
    vals.insert(0, None)
    return vals

e = modeller.Environ()
a = modeller.Alignment(e, file='results_aln/TvLDH-1bdmA.ali')

template = get_profile('results_models/1bdmA.profile', a['1bdmA'])
model = get_profile('results_models/TvLDH.profile', a['TvLDH'])

# Plot the template and model profiles in the same plot for comparison:
fig, ax = plt.subplots()
ax.set_xlabel('Alignment position')
ax.set_ylabel('DOPE per-residue score')
ax.plot(model, color='red', linewidth=2, label='Model')
ax.plot(template, color='green', linewidth=2, label='Template')
fig.legend()
fig.savefig('dope_profile.png', dpi=200)

6. Visualization

마지막으로 PyMOL을 이용해 템플릿 구조와 모델을 시각화하고, 두 구조를 비교 분석한다.

Views: 323

Leave a Comment

Your email address will not be published. Required fields are marked *

Scroll to Top