Hi,
I have a model of a tetramer, and I would like to evaluate/show only the DOPE profile for one chain, because the tetramer is symmetric.
Here is the script I used for calculating the DOPE profile:
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, 'model.B99990001.pdb')
# Assess with DOPE: s = selection(mdl) # all atom selection s.assess_dope(output='ENERGY_PROFILE NO_REPORT', file='model.profile', normalize_profile=True, smoothing_window=15)
and for plotting:
import pylab import modeller
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 = file(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 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='3zrs-align-nterm.ali')
model = get_profile('3zrs-nterm.profile', a['3zrs-nterm']) template = get_profile('3zrs-nterm.profile', a['3zrs-nterm'])
# Plot the template and model profiles in the same plot for comparison: pylab.figure(1, figsize=(10,6)) pylab.xlabel('Alignment position') pylab.ylabel('DOPE per-residue score') pylab.plot(model, color='red', linewidth=2, label='Model') pylab.plot(template, color='green', linewidth=2, label='Template') pylab.legend() pylab.savefig('dope_profile.png', dpi=65)
These scripts are from the tutorial on the website.
Best,
Matthias