IMP logo
IMP Reference Guide  2.21.0
The Integrative Modeling Platform
version 2.21.0
pmi/ambiguity.py
1 ## \example pmi/ambiguity.py
2 """This script shows how to create a system with multiple copies of the
3 same molecule.
4 We also create some cross-links which take into account the ambiguity.
5 The key to ambiguity is using the same molecule name for ambiguous copies.
6 That way when you perform Selection it automatically finds all relevant
7 molecules.
8 """
9 
10 import IMP
11 import IMP.atom
12 import IMP.algebra
13 import IMP.rmf
14 import IMP.pmi
15 import IMP.pmi.topology
16 import IMP.pmi.dof
17 import IMP.pmi.macros
18 import IMP.pmi.restraints.stereochemistry
19 import IMP.pmi.restraints.crosslinking
20 import ihm.cross_linkers
21 import tempfile
22 import os
23 import sys
24 
25 IMP.setup_from_argv(sys.argv, "ambiguity example")
26 
27 # ##################### SYSTEM SETUP #####################
28 
29 # Setup multistate system
30 mdl = IMP.Model()
31 s = IMP.pmi.topology.System(mdl)
32 st1 = s.create_state()
33 st2 = s.create_state()
34 
35 # For each state add some molecules - we'll make some bead only structures
36 # State 1: ProtA (chainA), ProtA (chainB), ProtC (chainC)
37 sequence = 'A'*10
38 m1A = st1.create_molecule('ProtA', sequence, chain_id='A')
39 m1A.add_representation(m1A, resolutions=[1])
40 # create_clone() will copy name/structure/representation
41 # You cannot edit it!
42 # There is also a function create_copy() which
43 # only copies the name, then you can change reps
44 m1B = m1A.create_clone(chain_id='B')
45 m1C = st1.create_molecule('ProtC', sequence, chain_id='C')
46 m1C.add_representation(m1C, resolutions=[1])
47 
48 # State 2: ProtA (chainA), ProtC (chainC)
49 m2A = st2.create_molecule('ProtA', sequence, chain_id='A')
50 m2A.add_representation(m2A, resolutions=[1])
51 m2C = st2.create_molecule('ProtC', sequence, chain_id='C')
52 m2C.add_representation(m2C, resolutions=[1])
53 root_hier = s.build()
54 
55 # Display all the states, molecules, representations
56 IMP.atom.show_with_representations(root_hier)
57 
58 # Setup all molecules to move as flexible beads and super rigid bodies
59 # "Super rigid bodies" aren't really rigid, it's just a mover that
60 # moves the whole body together
61 dof = IMP.pmi.dof.DegreesOfFreedom(mdl)
62 for mol in (m1A, m1B, m1C, m2A, m2C):
63  dof.create_flexible_beads(mol, max_trans=0.1)
64  dof.create_super_rigid_body(mol)
65 
66 
67 # ##################### RESTRAINTS #####################
68 output_objects = [] # keep a list of functions that need to be reported
69 
70 # Crosslinks setup
71 # 1) Create file. This one XL has 3 ambiguity options: State1 has 2,
72 # State2 has 1
73 lines = '''id,mol1,res1,mol2,res2,score
74 1,ProtA,3,ProtC,9,1.0
75 '''
76 tf = tempfile.NamedTemporaryFile(delete=False, mode='w')
77 tf.write(lines)
78 tf.close()
79 
80 # 2) Define the columns
81 kw = IMP.pmi.io.crosslink.CrossLinkDataBaseKeywordsConverter()
82 kw.set_unique_id_key("id")
83 kw.set_protein1_key("mol1")
84 kw.set_protein2_key("mol2")
85 kw.set_residue1_key("res1")
86 kw.set_residue2_key("res2")
87 kw.set_id_score_key("score")
88 xldb = IMP.pmi.io.crosslink.CrossLinkDataBase(kw)
89 xldb.create_set_from_file(tf.name)
90 os.remove(tf.name)
91 
92 # 3) Add the restraint
93 xlr = IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint(
94  root_hier=root_hier,
95  database=xldb,
96  length=21,
97  label="XL",
98  resolution=1,
99  slope=0.01,
100  linker=ihm.cross_linkers.dss)
101 xlr.add_to_model()
102 output_objects.append(xlr)
103 # needed to sample the nuisance particles (noise params)
104 dof.get_nuisances_from_restraint(xlr)
105 
106 # Connectivity keeps things connected along the backbone
107 crs = []
108 for mol in (m1A, m1B, m1C, m2A, m2C):
109  cr = IMP.pmi.restraints.stereochemistry.ConnectivityRestraint(mol)
110  cr.add_to_model()
111  output_objects.append(cr)
112 
113 # Excluded volume - one for each state (they don't interact)
114 evr1 = IMP.pmi.restraints.stereochemistry.ExcludedVolumeSphere(
115  included_objects=(m1A, m1B, m1C))
116 evr1.add_to_model()
117 output_objects.append(evr1)
118 evr2 = IMP.pmi.restraints.stereochemistry.ExcludedVolumeSphere(
119  included_objects=(m2A, m2C))
120 evr2.add_to_model()
121 output_objects.append(evr2)
122 
123 # ##################### SAMPLING #####################
124 # randomize particles a bit
125 IMP.pmi.tools.shuffle_configuration(root_hier, max_translation=20)
126 
127 # Shift state 2
128 # Even though the two states don't interact,
129 # it'll be easier to look at the RMF if we separate them
130 trans = IMP.algebra.Transformation3D([50, 0, 0])
131 for fb in IMP.core.get_leaves(m2A.get_hierarchy()) + \
132  IMP.core.get_leaves(m2C.get_hierarchy()):
133  IMP.core.transform(IMP.core.XYZ(fb), trans)
134 
135 # Run replica exchange Monte Carlo sampling
136 rex = IMP.pmi.macros.ReplicaExchange(
137  mdl,
138  root_hier=root_hier, # pass the root hierarchy
139  monte_carlo_sample_objects=dof.get_movers(), # pass MC movers
140  global_output_directory='ambiguity_output/',
141  output_objects=output_objects,
142  monte_carlo_steps=10,
143  number_of_frames=1) # increase number of frames to get better results!
144 rex.execute_macro()