Lesson 5 of 14 · 12 min

Building a model with SWISS-MODEL and MODELLER

Comparative modeling infers a target's 3D structure from a homologous template of known structure. Here we model Trypanosoma vivax lactate dehydrogenase (TvLDH) against the malate dehydrogenase template 1BDM chain A two ways: the fully automated SWISS-MODEL server and a scripted MODELLER run. We then rank the MODELLER candidates and compare the two outputs.

SWISS-MODEL web pipeline

SWISS-MODEL takes a target sequence, searches template libraries (or accepts a template you choose), builds the model with ProMod3, and scores it with GMQE and per-residue QMEANDisCo. The web form is point-and-click, but the same job runs headlessly through the REST API using a personal access token from your account page.

bash

export SM_TOKEN="your_access_token"

# Submit the target for automated template selection and modeling.
# target_sequences is a JSON array (one entry per chain).
PROJECT=$(curl -s https://swissmodel.expasy.org/automodel \
  -H "Authorization: Token ${SM_TOKEN}" \
  -H "Content-Type: application/json" \
  -d '{"target_sequences": ["MSEAAHVLITGAAGQIGYILSHWIASGELYGDRQVYLHLLDIPPAMNRLTALTMELEDCAFPHLAGFVATTDPKAAFKDIDCAFLVASMPLKPGQVRADLISSNSVIFKNTGEYLSKWAKPSVKVLVIGNPDNTNCEIAMLHAKNLKPENFSSLSMLDQNRAYYEVASKLGVDVKDVHDIIVWGNHGESMVADLTQATFTKEGKTQKVVDVLDHDYVFDTFFKKIGHRAWDILEHRGFTSAASPTKAAIQHMKAWLFGTAPGEVLSMGIPVPEGNPYGIKPGVVFSFPCNVDKEGKIHVVEGFKVNDWLREKLDFTEKDLFHEKEIALNHLAQGG"], "project_title": "TvLDH"}' \
  | python3 -c 'import sys,json; print(json.load(sys.stdin)["project_id"])')

# The summary endpoint returns the job status; repeat this call until it is COMPLETED
curl -s "https://swissmodel.expasy.org/project/${PROJECT}/models/summary/" \
  -H "Authorization: Token ${SM_TOKEN}"

{
  "status": "COMPLETED",
  "models": [
    {
      "model_id": 1,
      "gmqe": 0.74,
      "qmeandisco_global": 0.71,
      "template": "1bdm.1.A",
      "seq_id": 51.2,
      "coverage": 0.97,
      "coordinates_url": "https://swissmodel.expasy.org/project/aB3xY/models/01.pdb"
    }
  ]
}

Align with MODELLER

python

from modeller import *

env = Environ()

# Read the template structure (chain A of 1BDM) and start an alignment
aln = Alignment(env)
mdl = Model(env, file='1bdm', model_segment=('FIRST:A', 'LAST:A'))
aln.append_model(mdl, align_codes='1bdmA', atom_files='1bdm.pdb')

# Add the target sequence (PIR record in TvLDH.ali) and align to the structure
aln.append(file='TvLDH.ali', align_codes='TvLDH')
aln.align2d(max_gap_length=50)   # structure-aware gap placement

aln.write(file='TvLDH-1bdmA.ali', alignment_format='PIR')
aln.write(file='TvLDH-1bdmA.pap', alignment_format='PAP')

Build and rank models

python

from modeller import *
from modeller.automodel import *

env = Environ()

a = AutoModel(env,
              alnfile='TvLDH-1bdmA.ali',
              knowns='1bdmA',
              sequence='TvLDH',
              assess_methods=(assess.DOPE, assess.normalized_dope, assess.GA341))
a.starting_model = 1
a.ending_model = 5          # generate five candidate models
a.make()

# Keep successful models and rank by normalized DOPE (lower is better)
ok = [m for m in a.outputs if m['failure'] is None]
ok.sort(key=lambda m: m['Normalized DOPE score'])
best = ok[0]
print('Best: %s  nDOPE=%.3f  GA341=%.4f'
      % (best['name'], best['Normalized DOPE score'], best['GA341 score'][0]))

>> Summary of successfully produced models:
Filename                          molpdf    DOPE score  Normalized DOPE score  GA341 score
------------------------------------------------------------------------------------------
TvLDH.B99990001.pdb           1963.75073  -38121.855               -1.31245      1.00000
TvLDH.B99990002.pdb           1893.44714  -38323.715               -1.42088      1.00000
TvLDH.B99990003.pdb           2050.19507  -37962.113               -1.22765      1.00000
TvLDH.B99990004.pdb           1978.68042  -38210.660               -1.35401      1.00000
TvLDH.B99990005.pdb           1925.31104  -38401.980               -1.46333      1.00000

Best: TvLDH.B99990005.pdb  nDOPE=-1.463  GA341=1.0000

A normalized DOPE (z-DOPE) below -1 flags a model with native-like packing, and unlike raw DOPE it does not scale with chain length, so it compares fairly across targets. GA341 runs from 0 to 1, and values above 0.7 mark a reliable fold. Rank on normalized DOPE first, break ties with GA341 and molpdf, then superpose the winner on the SWISS-MODEL coordinates and read QMEANDisCo against your nDOPE.

Try it yourself: Fetch the template with wget -O 1bdm.pdb https://files.rcsb.org/download/1BDM.pdb and place the TvLDH sequence in a PIR file TvLDH.ali. Run align2d, generate five models with AutoModel, and pick the lowest normalized-DOPE candidate. Submit the same sequence to SWISS-MODEL through the API, superpose the two best models with TM-align or PyMOL's align command, and report the C-alpha RMSD. Then compare QMEANDisCo against normalized DOPE to decide which model you trust.

Homology modeling: principles and the pipeline Threading and fold recognition