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