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