IMP logo
IMP Reference Guide  2.21.0
The Integrative Modeling Platform
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
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()
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
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
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
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")
89 xldb.create_set_from_file(tf.name)
90 os.remove(tf.name)
91 
92 # 3) Add the restraint
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):
110  cr.add_to_model()
111  output_objects.append(cr)
112 
113 # Excluded volume - one for each state (they don't interact)
115  included_objects=(m1A, m1B, m1C))
116 evr1.add_to_model()
117 output_objects.append(evr1)
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
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()