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