Trajectory Analysis Quick Reference

You can execute various trajectory analyzers in three different environments:

  1. Maestro’s export movie, movie
  2. The analyze_simulation.py script, which reads an st2 file containing a list of analyzers and their parameters
  3. The Schrödinger Python API

 

Required for all analyzers

To execute any of the analyzers, first: 

  • Import output CMS into Maestro
  • Double-click ‘T’ next to the entry name

 

For every st2 file, include the following wrapper: 

Keywords = [
{Analyzer1}
{Analyzer2}
]

 

Python API

from schrodinger.application.desmond.packages import topo
from schrodinger.application.desmond.packages import traj_util
from schrodinger.application.desmond.packages import analysis

msys_model, cms_model, tr =
     traj_util.read_cms_and_traj("example-out.cms")

 

Analyzer

Trajectory Player

st2 file contents

Python API

Distance
  • Measure distance in Maestro Workspace
  • Plot → Measurements → Currently in Workspace
 
{Distance = { 
	a1 = 24
	a2 = 115
  }
}
 
dist = analysis.Distance(msys_model, cms_model, 24, 115)
results = analysis.analyze(tr, dist)

Angle
  • Measure angle in Maestro Workspace
  • Plot → Measurements → Currently in Workspace
 
{Angle = { 
	a1 = 24
	a2 = 115
	a3 = 127
  }
}
 
angle = analysis.Angle(msys_model, cms_model, 24, 115, 127)
results = analysis.analyze(tr, angle)

Dihedral
  • Measure dihedral in Maestro Workspace
  • Plot → Measurements → Currently in Workspace
 
{Torsion = {
	a1 = 24
	a2 = 115
	a3 = 127
	a4 = 131
  }
}
 
dihedral = analysis.Torsion(msys_model, cms_model,
				24, 115, 127, 131)
results = analysis.analyze(tr, dihedral)

Ligand RMSD
  • Select ligand in Workspace
  • Plot → Descriptors → RMSD
 
{RMSD = {
	ASL = "ligand"
	Frame = 0
	Type = "Ligand"
   }
}
 
ref_frame = 0
ligand_aids = cms_model.select_atom("ligand")
ligand_gids = topo.aids2gids(cms_model, ligand_aids,
		include_pseudoatoms=False)
ligand_ref_pos = tr[ref_frame].pos(ligand_gids)
lig_rmsd = analysis.LigandRMSD(msys_model, cms_model,
		ligand_aids, ligand_ref_pos)
results = analysis.analyze(tr, lig_rmsd)

Ligand RMSD; align on protein backbone
  • Select ligand in Workspace
  • Plot → Descriptors → RMSD Settings → set Other atoms to “protein and backbone”
  • Plot → Descriptors → RMSD
 
{RMSD = {
	ASL = "ligand"
	FitBy = "protein
		and backbone"
	Frame = 0
	Type = "Ligand"
   }
}
 
ref_frame = 0
ligand_aids = cms_model.select_atom("ligand")
ligand_gids = topo.aids2gids(cms_model, ligand_aids,
		include_pseudoatoms=False)
ligand_ref_pos = tr[ref_frame].pos(ligand_gids)
protein_aids = cms_model.select_atom("protein and backbone")
protein_gids = topo.aids2gids(cms_model, protein_aids,
		include_pseudoatoms=False)
protein_ref_pos = tr[ref_frame].pos(protein_gids)
lig_rmsd = analysis.LigandRMSD(msys_model, cms_model,
		ligand_aids, ligand_ref_pos,
		fit_aids=protein_aids,
		fit_ref_pos=protein_ref_pos)
results = analysis.analyze(tr, lig_rmsd)

Protein backbone RMSD
  • Select protein backbone in Workspace
  • Plot → Descriptors → RMSD
 
{RMSD = {
	ASL = "backbone"
	Frame = 0
   }
}
 
ref_frame = 0
protein_aids = cms_model.select_atom("backbone")
protein_gids = topo.aids2gids(cms_model, protein_aids,
		include_pseudoatoms=False)
protein_ref_pos = tr[ref_frame].pos(protein_gids)
protein_rmsd = analysis.RMSD(msys_model, cms_model,
		protein_aids, protein_ref_pos)
results = analysis.analyze(tr, protein_rmsd)

Protein backbone RMSF
  • Select protein backbone in Workspace
  • Plot → RMSF → Per Residue
 
{RMSF = {
	ASL = "backbone"
	Frame = 0
	Tab = "p_rmsf_tab"
   }
}
 
ref_frame = 0
protein_aids = cms_model.select_atom("backbone")
protein_gids = topo.aids2gids(cms_model, protein_aids,
		include_pseudoatoms=False)
protein_ref_pos = tr[ref_frame].pos(protein_gids)
protein_rmsf = analysis.ProteinRMSF(msys_model, cms_model,
		protein_aids,
		protein_aids,
		protein_ref_pos)
results = analysis.analyze(tr, protein_rmsf)

Ligand atoms RMSF
  • Select protein backbone in Workspace
  • Plot → RMSF → Per Atom
 
{RMSF = {
	ASL = "ligand"
	Frame = 0
   }
}
 
ref_frame = 0
ligand_aids = cms_model.select_atom("ligand")
lig_rmsf = analysis.RMSF(msys_model, cms_model, ligand_aids,
		None, None)
results = analysis.analyze(tr, lig_rmsf)

Radius of gyration of a ligand
  • Select ligand in Workspace
  • Plot → Descriptors → Radius of Gyration
 
{Rad_Gyration = {
	ASL = "ligand"
   }
}
 
rgyr = analysis.Gyradius(msys_model, cms_model,
		asl="ligand")
results = analysis.analyze(tr, rgyr)

Molecular surface area of a ligand
  • Select ligand in Workspace
  • Plot → Descriptors → Molecular Surface Area
 
{Molecular_Surface_Area = {
	ASL = "ligand"
   }
}
 
molsa = analysis.MolecularSurfaceArea(msys_model, cms_model,
		"ligand")
results = analysis.analyze(tr, molsa)

Polar surface area of a ligand
  • Select ligand in Workspace
  • Plot → Descriptors → Polar Surface Area
 
{Polar_Surface_Area = {
	ASL = "ligand"
   }
}
 
psa = analysis.PolarSurfaceArea(msys_model, cms_model,
		"ligand")
results = analysis.analyze(tr, psa)

Solvent accessible surface area of a ligand
  • Select ligand in Workspace
  • Plot → Descriptors → Solvent Accessible Surface Area
 
{SA_Surface_Area = {
	ASL = "ligand"
	Exclude_ASL = "water"
   }
}
 
sasa = analysis.SolventAccessibleSurfaceArea(
		msys_model, cms_model, "ligand", exclude_asl="water")
results = analysis.analyze(tr, sasa)

Hydrogen bonds between protein and ligand
  • Plot → Interaction Counts→ Hydrogen Bonds
 
{HBonds = {
	ASL1 = "protein"
	ASL2 = "ligand"
   }
}
 
protein_aids = cms_model.select_atom("protein")
ligand_aids = cms_model.select_atom("ligand")
hbonds = analysis.HydrogenBondFinder(msys_model, cms_model,
		protein_aids, ligand_aids)
results = analysis.analyze(tr, hbonds)

Planar angle
  • Measure planar angle in Maestro Workspace
  • Plot → Measurements → Planar Angle
 
{Planar_Angle = {
	a1 = 24
	a2 = 115
	a3 = 127
	a4 = 131
	a5 = 142
	a6 = 156
   }
}
 
planar = analysis.PlanarAngle(msys_model, cms_model,
		24, 115, 127, 131, 142, 156)
results = analysis.analyze(tr, planar)

Fitted planar angle

 
{Fitted_Planar_Angle = {
	ASL1 = "res.num 10"
	ASL2 = "res.num 20"
   }
}
 
aids1 = cms_model.select_atom("res.num 10")
aids2 = cms_model.select_atom("res.num 20")
fp = analysis.FittedPlanarAngle(msys_model, cms_model,
		xids0=aids1, xids1=aids2)
results = analysis.analyze(tr, fp)

Ramachandran (Phi/Psi torsions)
  • Plot → Phi-Psi Angles (Selected Residues)
 
{Ramachandran = {
	ASL = "backbone"
   }
}
 
rama = analysis.Ramachandran(msys_model, cms_model,
		"backbone")
results = analysis.analyze(tr, rama)

Secondary structure

 
{SecondaryStructure = {
	ASL = backbone
   }
}
 
ca_aids = cms_model.select_atom(
		'fillres(backbone) AND (atom.ptype " CA ")')
ss = analysis.SecondaryStructure(msys_model, cms_model,
		aids=ca_aids)
results = analysis.analyze(tr, ss)

Ligand hydrogen bonds

 
{LigandHBonds = {
	ASL1 = "ligand"
   }
}
 
ligand_aids = cms_model.select_atom("ligand")
hbonds = analysis.HydrogenBondFinder(msys_model, cms_model,
		aids1=ligand_aids, aids2=ligand_aids)
results = analysis.analyze(tr, hbonds)

Protein-ligand interactions

 
{ProtLigInter = {
	ASL = "protein"
	LigandASL = "ligand"
	MetalASL = "metals"
   }
}
 
pli = analysis.ProtLigInter(msys_model, cms_model,
		"protein", "ligand")
results = analysis.analyze(tr, pli)

Protein-protein interactions

 
{ProtProtInter = {
	ASL = "protein"
   }
}
 
ppi = analysis.ProtProtInter(msys_model, cms_model,
		"protein")
results = analysis.analyze(tr, ppi)

SASA by protein residue

 
{SASA_By_Protein_Residue = {
	ASL = solute
	Resolution = 0.2
   }
}
 
sasa_res = analysis.SolventAccessibleSurfaceAreaByResidue(
		msys_model, cms_model, "protein")
results = analysis.analyze(tr, sasa_res)

Molecule energy
  • Plot → Energy Calculations
 
{Energy = {
	Molecule_Number = 1
	Type = "Potential"
	Type2 = "Total"
   }
}
 
from schrodinger.application.desmond.packages \
		import energygroup
me = energygroup.MoleculeEnergy(cms_model,
		molecule_number=1,
		type="Potential", type2="Total")

Interaction energy

 
{InteractionEnergy = {
	ASL1 = "protein"
	ASL2 = "ligand"
   }
}
 
from schrodinger.application.desmond.packages \
		import energygroup
ie = energygroup.InteractionEnergy(cms_model,
		"protein", "ligand")

Molecule energy matrix

 
{MoleculeEnergyMatrix = {
	Molecule_Numbers = [1, 2]
   }
}
 
from schrodinger.application.desmond.packages \
		import energygroup
mem = energygroup.MoleculeEnergyMatrix(cms_model,
		molecule_numbers=[1, 2])

Selection energy matrix

 
{SelectionEnergyMatrix = {
	Asl_Selections = [
		"protein"
		"ligand"
	]
   }
}
 
from schrodinger.application.desmond.packages \
		import energygroup
sem = energygroup.SelectionEnergyMatrix(cms_model,
		asl_selections=["protein", "ligand"])

Partial cohesive energy

 
{PartialCohesiveEnergy = {
	ASL = "res.num 1-100"
   }
}
 
from schrodinger.application.desmond.packages \
		import energygroup
pce = energygroup.PartialCohesiveEnergy(cms_model,
		"res.num 1-100")

Bulk properties

 
{Bulk = {
	Type = "Potential"
	SpecificHeatWithLastNPercent = 70
   }
}
 
from schrodinger.application.desmond.packages \
		import energygroup
bulk = energygroup.Bulk(cms_model,
		type="Potential", last_n_percent=70)

Area per lipid

 
{Area_Per_Lipid = {
	MembraneASL = "membrane"
   }
}
 
apl = analysis.AreaPerLipid(msys_model, cms_model,
		membrane_asl="membrane")
results = analysis.analyze(tr, apl)

Membrane thickness

 
{Membrane_Thickness = {
	MembraneASL = "membrane"
	ProxyASL = "atom.ele P"
   }
}
 
mt = analysis.MembraneThickness(msys_model, cms_model,
		membrane_asl="membrane",
		proxy_asl="atom.ele P")
results = analysis.analyze(tr, mt)

Membrane density profile

 
{Membrane_Density_Profile = {
	MembraneASL = "membrane"
	DensityASL = "lipids"
	Zmin = -40
	Zmax = 40
	SliceHeight = 1.0
   }
}
 
mdp = analysis.MembraneDensityProfile(msys_model,
		cms_model, membrane_asl="membrane",
		density_sel_asl="lipids",
		z_min=-40, z_max=40, slice_height=1.0)
results = analysis.analyze(tr, mdp)

Order parameter

 
{Order_Parameter = {
	ASL = "membrane"
	DirectorType = "Z"
	DescriptorType = "lipid_CH"
   }
}
 
director = analysis.AxisDirector("Z")
descriptor = analysis.LipidDirector(msys_model,
		cms_model, "membrane", "CH")
op = analysis.OrderParameter(director, descriptor,
		analysis.reduce_lipid_vec_list)
results = analysis.analyze(tr, op)

Amorphous crystal interactions

 
{AmorphousCrystalInter = {
	ASL = "res.num 1"
   }
}
 
aci = analysis.AmorphousCrystalInter(msys_model,
		cms_model, "res.num 1")
results = analysis.analyze(tr, aci)