1 """@namespace IMP.pmi.nonmaintained
5 from __future__
import print_function
11 def __init__(self, m):
14 self.rigid_bodies = []
15 self.floppy_bodies = []
16 self.maxtrans_rb = 2.0
17 self.maxtrans_fb = 2.0
20 def add_protein(self, name, res):
21 (firstres, lastres) = res
22 from math
import pi, cos, sin
25 nres = lastres - firstres
26 radius = (nres) * 5 / 2 / pi
28 for res
in range(firstres, lastres):
29 alpha = 2 * pi / nres * (res - firstres)
30 x = radius * cos(alpha)
31 y = radius * sin(alpha)
36 d.set_coordinates_are_optimized(
True)
38 self.hier.add_child(h)
40 def get_hierarchy(self):
43 def set_rod(self, chainname, res):
44 (firstres, lastres) = res
49 residue_indexes=list(range(
52 ps = sel.get_selected_particles()
54 self.rigid_bodies.append(rb)
56 def get_particles_to_sample(self):
61 ps[
"Floppy_Bodies_Rods"] = (self.floppy_bodies, self.maxtrans_fb)
62 ps[
"Rigid_Bodies_Rods"] = (
71 def __init__(self, m):
76 self.floppy_bodies = []
77 self.maxtrans_fb = 0.2
78 self.particle_database = {}
80 def add_bead(self, radius, label="None", color=None):
84 self.particle_database[label] = p
85 self.floppy_bodies.append(p)
89 d.set_coordinates_are_optimized(
True)
96 self.hier.add_child(p)
98 self.set_color(label, color)
99 return self.particle_database[label]
101 def set_color(self, label, value):
102 p = self.particle_database[label]
106 def set_floppy_bodies_max_trans(self, maxtrans):
107 self.maxtrans_fb = maxtrans
109 def get_hierarchy(self):
112 def get_bead(self, label):
113 return self.particle_database[label]
115 def set_maxtrans_fb(self, maxtrans_fb):
116 self.maxtrans_fb = maxtrans_fb
118 def get_particles_to_sample(self):
123 ps[
"Floppy_Bodies_Beads"] = (self.floppy_bodies, self.maxtrans_fb)
127 class MultipleStates(object):
129 def __init__(self, nstates, m):
130 global itertools, tools, restraints
136 self.floppy_bodies = []
137 self.rigid_bodies = []
138 for ncopy
in range(nstates):
139 self.floppy_bodies.append([])
140 self.rigid_bodies.append([])
142 self.rigid_bodies_are_sampled =
True
143 self.floppy_bodies_are_sampled =
True
146 self.prot_lowres = {}
147 self.nstates = nstates
151 self.xyzmodellist = []
153 self.maxtrans_rb = 0.15
154 self.maxrot_rb = 0.03
155 self.maxtrans_fb = 0.15
161 def set_label(self, label):
164 def set_rigid_bodies_are_sampled(self, input=True):
165 self.rigid_bodies_are_sampled = input
167 def set_floppy_bodies_are_sampled(self, input=True):
168 self.floppy_bodies_are_sampled = input
170 def get_rigid_bodies(self):
171 return self.rigid_bodies
173 def set_rigid_bodies_max_trans(self, maxtrans):
174 self.maxtrans_rb = maxtrans
176 def set_rigid_bodies_max_rot(self, maxrot):
177 self.maxrot_rb = maxrot
179 def set_floppy_bodies_max_trans(self, maxtrans):
180 self.maxtrans_fb = maxtrans
182 def get_hierarchies(self):
185 def destroy_residues(self, segments):
188 for prot
in self.prot:
189 for segment
in segments:
191 if (segment[0] == -1
or segment[1] == -1):
197 residue_indexes=list(range(segment[0],
199 for p
in s.get_selected_particles():
201 print(
"MultipleStates: one particle was not destroied because it was a RigidMember.")
211 def add_residues_to_chains(
214 residue_type=IMP.atom.LYS):
218 for rc
in residuechainlist:
224 atom_type=IMP.atom.AT_CA)
226 print(s.get_selected_particles())
227 if len(s.get_selected_particles()) == 0:
228 for prot
in self.prot:
229 print(
"adding " + str(rc))
234 d.set_coordinates_are_optimized(
True)
243 print(tools.get_residue_index_and_chain_from_particle(a))
246 p = s.get_selected_particles()[0]
248 print(rc, s.get_selected_particles()[0])
255 atom_type=IMP.atom.AT_CA)
257 print(s.get_selected_particles())
259 def add_beads(self, segments, xyzs=None, radii=None, colors=None):
261 this method generate beads in missing portions.
262 The segments argument must be a list of selections
263 in the form [(firstres,lastres,chain)]
264 each selection will generate a bead
275 for n, s
in enumerate(segments):
280 for prot
in self.prot:
281 for prot
in self.prot:
282 cps = IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
285 if chain.get_id() == chainid:
288 rindexes = list(range(firstres, lasteres + 1))
289 f.set_residue_indexes(rindexes)
291 "Fragment_" +
'%i-%i' %
294 mass = len(rindexes) * 110.0
296 if n + 1 > len(radii):
297 mass = len(rindexes) * 110.0
299 radius = (3 * vol / math.pi) ** (1 / 3)
303 if n + 1 > len(xyzs):
312 if n + 1 <= len(colors):
316 d = IMP.atom.XYZR.setup_particle(
319 def renumber_residues(self, chainid, newfirstresiduenumber):
320 for prot
in self.prot:
321 cps = IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
324 ps = c.get_children()
327 offs = newfirstresiduenumber - ri
331 r.set_index(ri + offs)
333 def destroy_everything_but_the_residues(self, segments):
335 for prot
in self.prot:
337 for segment
in segments:
340 if (segment[0] == -1
or segment[1] == -1):
346 residue_indexes=list(range(segment[0],
348 pstokeep += s.get_selected_particles()
351 if p
not in pstokeep:
353 print(
"MultipleStates: one particle was not destroied because it was a RigidMember.")
362 def generate_linkers_restraint_and_floppy_bodies(self, segment):
364 this methods automatically links the particles consecutively
365 according to the sequence. The restraint applied is a harmonic upper bound,
366 with a distance that is proportional to the number of residues
371 linker_restraint_objects = []
372 for ncopy, prot
in enumerate(self.prot):
373 if (segment[0] == -1
or segment[1] == -1):
379 residue_indexes=list(range(segment[0],
382 for p
in s.get_selected_particles():
384 (r, c) = tools.get_residue_index_and_chain_from_particle(p)
389 (r, c) = tools.get_residue_index_and_chain_from_particle(p)
390 p.set_name(str(r) +
":" + str(c))
391 tools.set_floppy_body(p)
392 self.floppy_bodies[ncopy].append(p)
394 residue_indexes.append((r, Floppy, c, p))
396 residue_indexes.sort()
398 pruned_residue_list = []
399 r0 = residue_indexes[0]
400 pruned_residue_list.append(r0)
404 for i
in range(1, len(residue_indexes)):
405 r = residue_indexes[i]
409 pruned_residue_list.append(r0)
410 pruned_residue_list.append(r)
412 elif r[1] != r0[1]
and r0[1] ==
False:
413 pruned_residue_list.append(r0)
414 pruned_residue_list.append(r)
416 elif r[1] == r0[1]
and r0[1]:
417 pruned_residue_list.append(r)
419 elif r[1] != r0[1]
and r[1] ==
False:
420 pruned_residue_list.append(r)
423 r0 = pruned_residue_list[0]
426 for i
in range(1, len(pruned_residue_list)):
427 r = pruned_residue_list[i]
431 linkdomaindef.append((r0[0], r[0], r[2]))
434 print(
" creating linker between atoms defined by: " + str(linkdomaindef))
436 ld = restraints.LinkDomains(prot, linkdomaindef, 1.0, 3.0)
437 ld.set_label(str(ncopy))
439 linker_restraint_objects.append(ld)
442 return linker_restraint_objects
444 def get_ref_hierarchies(self):
447 def get_number_of_states(self):
450 def get_rigid_bodies(self):
452 for rbl
in self.rigid_bodies:
457 def get_floppy_bodies(self):
459 for fbl
in self.floppy_bodies:
464 def set_rigid_bodies(self, rigid_body_list):
465 if len(self.prot) == 0:
466 print(
"MultipleStates.set_rigid_bodies: hierarchy was not initialized")
468 for ncopy, prot
in enumerate(self.prot):
470 for element
in rigid_body_list:
472 for interval
in element:
476 if (interval[0] == -1
or interval[1] == -1):
482 residue_indexes=list(range(interval[0],
484 for p
in s.get_selected_particles():
488 for key
in self.prot_lowres:
489 if (interval[0] == -1
or interval[1] == -1):
491 self.prot_lowres[key][ncopy],
495 self.prot_lowres[key][
496 ncopy], chains=interval[2],
497 residue_indexes=list(range(interval[0], interval[1] + 1)))
498 for p
in s.get_selected_particles():
504 rb.set_name(str(element))
507 print(
"MultipleStates.set_rigid_bodies: selection " + str(interval) +
" has zero elements")
508 self.rigid_bodies[ncopy] += rbl
510 def set_floppy_bodies(self, floppy_body_list):
513 if len(self.prot) == 0:
514 print(
"MultipleStates: hierarchy was not initialized")
517 for ncopy, prot
in enumerate(self.prot):
519 for element
in floppy_body_list:
521 for interval
in element:
525 if (interval[0] == -1
or interval[1] == -1):
531 residue_indexes=list(range(interval[0],
533 for p
in s.get_selected_particles():
534 (r, c) = tools.get_residue_index_and_chain_from_particle(
536 tools.set_floppy_body(p)
537 p.set_name(str(r) +
":" + str(c))
539 self.floppy_bodies[ncopy] += atoms
541 def get_particles_to_sample(self):
546 rblist = self.get_rigid_bodies()
547 fblist = self.get_floppy_bodies()
548 if self.rigid_bodies_are_sampled:
549 ps[
"Rigid_Bodies_MultipleStates"] = (
553 if self.floppy_bodies_are_sampled:
554 ps[
"Floppy_Bodies_MultipleStates"] = (fblist, self.maxtrans_fb)
557 def set_hierarchy_from_pdb(self, pdblistoflist):
558 "the input is a list of list of pdbs"
559 "one list for each copy"
561 for copy
in range(0, self.nstates):
562 prot = self.read_pdbs(pdblistoflist[copy])
563 self.prot.append(prot)
565 self.xyzmodellist.append(xyz)
567 def set_ref_hierarchy_from_pdb(self, pdblistoflist):
568 "the input is a list of list of pdbs"
569 "one list for each copy"
571 for copy
in range(0, self.nstates):
572 prot = self.read_pdbs(pdblistoflist[copy])
573 self.refprot.append(prot)
575 self.xyzreflist.append(xyz)
577 def read_pdbs(self, list_pdb_file):
578 """read pdbs from an external list file
579 create a simplified representation
580 if the pdbs are given a individual strings, it will read the
581 pdbs and give the chain name as specified in the pdb
582 If it s a tuple like (filename,chainname) it will read
583 the pdb and assing a name chainname
588 for pdb
in list_pdb_file:
596 #destroy CA atoms, for the future
597 for p in IMP.atom.get_leaves(h):
598 coor=IMP.core.XYZ(p).get_coordinates()
599 r=IMP.atom.Hierarchy(p).get_parent()
600 IMP.core.XYZ.setup_particle(r,coor)
604 cps = IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE)
607 #consolidate the chains
610 s0=IMP.atom.Selection(hier, chains=cid)
612 p=s0.get_selected_particles()[0]
613 re=IMP.atom.Residue(IMP.atom.Atom(p).get_parent()
614 ch=IMP.atom.Chain(re).get_parent())
621 if type(pdb) == tuple:
628 #destroy CA atoms, for the future
629 for p in IMP.atom.get_leaves(h):
630 coor=IMP.core.XYZ(p).get_coordinates()
631 r=IMP.atom.Hierarchy(p).get_parent()
632 IMP.core.XYZ.setup_particle(r,coor)
636 cps = IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE)
643 def recenter(self, prot):
644 "recenter the hierarchy"
646 center = IMP.algebra.get_zero_vector_3d()
652 d.set_coordinates(d.get_coordinates() - center)
653 d.set_coordinates_are_optimized(
True)
656 # bug generating code: keeping it for history
658 rb=IMP.atom.create_rigid_body(prot)
659 rbcoord=rb.get_coordinates()
660 rot=IMP.algebra.get_identity_rotation_3d()
661 tmptrans=IMP.algebra.Transformation3D(rot,rbcoord)
662 trans=tmptrans.get_inverse()
663 IMP.core.transform(rb,trans)
664 IMP.core.RigidBody.teardown_particle(rb)
665 self.m.remove_particle(rb)
668 def shuffle_configuration(self, bounding_box_length):
669 "shuffle configuration, used to restart the optimization"
670 "it only works if rigid bodies were initialized"
671 if len(self.rigid_bodies) == 0:
672 print(
"MultipleStates: rigid bodies were not intialized")
673 hbbl = bounding_box_length / 2
674 for rbl
in self.rigid_bodies:
684 rb.set_reference_frame(
687 def generate_simplified_hierarchy(self, nres):
689 self.prot_lowres[nres] = []
690 for prot
in self.prot:
698 self.prot_lowres[nres].append(sh)
700 def get_simplified_hierarchy(self, nres):
701 return self.prot_lowres[nres]
703 def calculate_drms(self):
706 if len(self.xyzmodellist) == 0:
707 print(
"MultipleStates: hierarchies were not intialized")
709 if len(self.xyzreflist) == 0:
710 print(
"MultipleStates: reference hierarchies were not intialized")
713 for i
in range(len(self.xyzreflist)):
714 for j
in range(len(self.xyzmodellist)):
717 self.xyzmodellist[j],
720 drmsdval = tools.get_drmsd(
721 self.xyzmodellist[j],
723 drmsd[
"MultipleStates_DRMSD_" +
724 str(i) +
"-Model_" + str(j)] = drmsdval
728 for assign
in itertools.permutations(list(range(len(self.xyzreflist)))):
730 for i, j
in enumerate(assign):
731 s += drmsd[
"MultipleStates_DRMSD_" +
732 str(j) +
"-Model_" + str(i)]
735 drmsd[
"MultipleStates_Total_DRMSD"] = min(min_drmsd)
738 def get_output(self):
740 if len(self.refprot) != 0:
741 drms = self.calculate_drms()
743 output[
"MultipleStates_Total_Score_" +
744 self.label] = str(self.m.evaluate(
False))
750 class LinkDomains(object):
752 def __init__(self, prot, resrangelist, kappa, length=5.0):
761 self.resrangelist = resrangelist
765 self.m = self.prot.get_model()
768 for pair
in self.resrangelist:
778 atom_type=IMP.atom.AT_CA)
779 p0 = s0.get_selected_particles()[0]
788 atom_type=IMP.atom.AT_CA)
789 p1 = s1.get_selected_particles()[0]
795 dist0 = float(pair[1] - pair[0]) * self.length
800 "LinkDomains_" + str(pair[0]) +
"-" + str(pair[1]) +
"_" + str(pair[2]))
801 self.rs.add_restraint(pr)
802 self.pairs.append((p0, p1, r0, c0, r1, c1, pr))
804 def set_label(self, label):
807 def add_to_model(self):
808 self.m.add_restraint(self.rs)
813 def get_hierarchy(self):
819 def get_residue_ranges(self):
820 return self.resrangelist
822 def get_length(self):
825 def get_restraint(self):
828 def get_restraints(self):
830 for r
in self.rs.get_restraints():
831 rlist.append(IMP.core.PairRestraint.get_from(r))
834 def print_chimera_pseudobonds(self, filesuffix, model=0):
835 f = open(filesuffix +
".chimera",
"w")
838 s =
"#" + str(model) +
":" + str(p[2]) +
"." + p[3] +
"@" + atype + \
839 " #" + str(model) +
":" + str(p[4]) +
"." + p[5] +
"@" + atype
843 def get_output(self):
846 score = self.rs.unprotected_evaluate(
None)
847 output[
"_TotalScore"] = str(score)
848 output[
"LinkDomains_" + self.label] = str(score)
849 for rst
in self.rs.get_restraints():
853 output[
"LinkDomains_" + rst.get_name() +
854 "_" + self.label] = rst.unprotected_evaluate(
None)
856 for i
in range(len(self.pairs)):
858 p0 = self.pairs[i][0]
859 p1 = self.pairs[i][1]
860 r0 = self.pairs[i][2]
861 c0 = self.pairs[i][3]
862 r1 = self.pairs[i][4]
863 c1 = self.pairs[i][5]
865 label = str(r0) +
":" + c0 +
"_" + str(r1) +
":" + c1
869 output[
"LinkDomains_Distance_" + label +
"_" +
878 class UpperBound(object):
880 def __init__(self, prot, respairs, kappa, length=5.0):
890 self.respairs = respairs
894 self.m = self.prot.get_model()
896 for pair
in self.respairs:
901 residue_index=pair[0],
902 atom_type=IMP.atom.AT_CA)
903 p0 = s0.get_selected_particles()[0]
911 residue_index=pair[2],
912 atom_type=IMP.atom.AT_CA)
913 p1 = s1.get_selected_particles()[0]
924 "UpperBound_" + str(pair[0]) +
"-" + str(pair[1]) +
"_" + str(pair[2]))
925 self.rs.add_restraint(pr)
927 def set_label(self, label):
930 def add_to_model(self):
931 self.m.add_restraint(self.rs)
933 def get_hierarchy(self):
939 def get_residue_ranges(self):
942 def get_length(self):
945 def get_restraint(self):
948 def get_output(self):
951 score = self.rs.unprotected_evaluate(
None)
952 output[
"_TotalScore"] = str(score)
953 output[
"UpperBound_" + self.label] = str(score)
959 class ExcludedVolumeResidue(object):
961 def __init__(self, prot, kappa):
966 self.m = self.prot.get_model()
968 atoms = IMP.atom.get_by_type(self.prot, IMP.atom.ATOM_TYPE)
976 lsa.add_particles(atoms)
979 self.rs.add_restraint(evr)
981 def set_label(self, label):
984 def add_excluded_particle_pairs(self, excluded_particle_pairs):
987 lpc.add_particle_pairs(excluded_particle_pairs)
989 IMP.core.ExcludedVolumeRestraint.get_from(
990 self.rs.get_restraints()[0]).add_pair_filter(icpf)
992 def add_to_model(self):
993 self.m.add_restraint(self.rs)
995 def get_hierarchy(self):
1001 def get_restraint(self):
1004 def get_output(self):
1007 score = self.rs.unprotected_evaluate(
None)
1008 output[
"_TotalScore"] = str(score)
1009 output[
"ExcludedVolumeResidue_" + self.label] = str(score)
1015 class BipartiteExcludedVolumeResidue(object):
1017 def __init__(self, prot1, prot2, kappa):
1023 self.m = self.prot.get_model()
1025 atoms1 = IMP.atom.get_by_type(prot1, IMP.atom.ATOM_TYPE)
1027 ls1.add_particles(atoms1)
1035 atoms2 = IMP.atom.get_by_type(prot2, IMP.atom.ATOM_TYPE)
1037 ls2.add_particles(atoms2)
1052 self.rs.add_restraint(evr3)
1053 self.m.add_restraint(rs)
1055 def set_label(self, label):
1058 def add_to_model(self):
1059 self.m.add_restraint(self.rs)
1061 def get_hierarchy(self):
1062 return self.prot1, self.prot2
1064 def get_kappa(self):
1067 def get_restraint(self):
1070 def get_output(self):
1073 score = self.rs.unprotected_evaluate(
None)
1074 output[
"_TotalScore"] = str(score)
1075 output[
"BipartiteExcludedVolumeResidue_" + self.label] = str(score)
1081 class TemplateRestraint(object):
1083 def __init__(self, ps1, ps2, cutoff=6.5, kappa=1.0, forcerb=False):
1084 self.m = ps1[0].get_model()
1086 self.cutoff = cutoff
1089 self.forcerb = forcerb
1100 if dist <= self.cutoff:
1104 self.rs.add_restraint(pr)
1106 def set_label(self, label):
1109 def add_to_model(self):
1110 self.m.add_restraint(self.rs)
1112 def get_cutoff(self):
1115 def get_kappa(self):
1118 def get_restraint(self):
1121 def get_output(self):
1124 score = self.rs.unprotected_evaluate(
None)
1125 output[
"_TotalScore"] = str(score)
1126 output[
"TemplateRestraint_" + self.label] = str(score)
1132 class MarginalChi3Restraint(object):
1134 def __init__(self, part1, part2):
1135 global impisd2, tools
1136 import IMP.isd2
as impisd2
1139 self.m = part1.get_model()
1142 self.sigmamaxtrans = 0.1
1146 self.sigma = tools.SetupNuisance(
1154 for i
in range(len(self.ps1)):
1155 mc = impisd2.MarginalChi3Restraint(
1159 self.rs.add_restraint(mc)
1161 def set_label(self, label):
1164 def add_to_model(self):
1165 self.m.add_restraint(self.rs)
1167 def get_restraint(self):
1170 def get_particles_to_sample(self):
1172 ps[
"Nuisances_MarginalChi3Restraint_Sigma_" +
1173 self.label] = ([self.sigma], self.sigmamaxtrans)
1176 def get_output(self):
1179 score = self.rs.unprotected_evaluate(
None)
1180 output[
"_TotalScore"] = str(score)
1181 output[
"MarginalChi3Restraint_" + self.label] = str(score)
1182 output[
"MarginalChi3Restraint_Sigma_" +
1183 self.label] = str(self.sigma.get_scale())
1192 this class initializes a CrossLinkMS restraint and contains
1193 all useful information, such as the cross-link database, contained in self.pairs
1194 If restraint_file=None, it will proceed creating simulated data
1197 def __init__(self, prots,
1198 listofxlinkertypes=[
"BS3",
"BS2G",
"EGS"], map_between_protein_names_and_chains=
None,
1199 sigmamin=1.0, sigmamax=1.0, sigmagrid=1, sigmaissampled=
False, typeofprofile=
"gofr"):
1200 global impisd2, tools
1201 import IMP.isd2
as impisd2
1204 if map_between_protein_names_and_chains
is None:
1205 map_between_protein_names_and_chains = {}
1212 self.m = self.prots[0].get_model()
1213 self.sigmamin = sigmamin
1214 self.sigmamax = sigmamax
1215 self.sigmagrid = sigmagrid
1216 self.sigmaissampled = sigmaissampled
1218 self.sigmatrans = 0.1
1219 self.sigmaglobal = tools.SetupNuisance(self.m, self.sigmamin,
1220 self.sigmamin, self.sigmamax, self.sigmaissampled).get_particle()
1221 self.outputlevel =
"low"
1222 self.listofxlinkertypes = listofxlinkertypes
1224 self.reaction_rates =
None
1225 self.allpairs_database =
None
1226 self.residue_list =
None
1228 self.crosslinker_dict = self.get_crosslinker_dict(typeofprofile)
1231 self.mbpnc = map_between_protein_names_and_chains
1235 def get_crosslinker_dict(self, typeofprofile="gofr"):
1240 disttuple = (0.0, 200.0, 1000)
1241 omegatuple = (1.0, 1000.0, 30)
1242 sigmatuple = (self.sigmamin, self.sigmamax, self.sigmagrid)
1244 crosslinker_dict = {}
1245 if "BS3" in self.listofxlinkertypes:
1246 crosslinker_dict[
"BS3"] = tools.get_cross_link_data(
"bs3l",
1247 "pmf_bs3l_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1248 don=
None, doff=
None, prior=0, type_of_profile=typeofprofile)
1249 if "BS2G" in self.listofxlinkertypes:
1250 crosslinker_dict[
"BS2G"] = tools.get_cross_link_data(
"bs2gl",
1251 "pmf_bs2gl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1252 don=
None, doff=
None, prior=0, type_of_profile=typeofprofile)
1253 if "EGS" in self.listofxlinkertypes:
1254 crosslinker_dict[
"EGS"] = tools.get_cross_link_data(
"egl",
1255 "pmf_egl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1256 don=
None, doff=
None, prior=0, type_of_profile=typeofprofile)
1257 if "Short" in self.listofxlinkertypes:
1259 crosslinker_dict[
"Short"] = tools.get_cross_link_data_from_length(
1264 return crosslinker_dict
1270 if restraint_file
is None:
1272 restraint_list = self.allpairs_database
1276 f = open(restraint_file)
1277 restraint_list = f.readlines()
1281 self.added_pairs_list = []
1282 self.missing_residues = []
1283 for line
in restraint_list:
1286 force_restraint =
False
1288 if restraint_file
is None:
1289 if line[
"Is_Detected"]:
1290 crosslinker = line[
"Crosslinker"]
1291 (r1, c1) = line[
"Identified_Pair1"]
1292 (r2, c2) = line[
"Identified_Pair2"]
1298 tokens = line.split()
1300 if (tokens[0] ==
"#"):
1306 crosslinker = tokens[4]
1309 self.index = int(tokens[5])
1313 if (tokens[len(tokens) - 1] ==
"F"):
1314 force_restraint =
True
1318 totallist = eval(line)
1319 self.add_crosslink_according_to_new_file(totallist)
1323 print(
'''CrossLinkMS: attempting to add restraint between
1324 residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
1334 atom_type=IMP.atom.AT_CA)
1335 p1 = (s1.get_selected_particles()[0])
1337 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1344 atom_type=IMP.atom.AT_CA)
1345 p2 = (s2.get_selected_particles()[0])
1347 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1350 for copy
in self.prots:
1355 atom_type=IMP.atom.AT_CA)
1356 p1s.append(s1.get_selected_particles()[0])
1361 atom_type=IMP.atom.AT_CA)
1362 p2s.append(s2.get_selected_particles()[0])
1369 print(
'''CrossLinkMS: WARNING> residue %d of chain %s and
1370 residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
1375 if (p1s[0], p2s[0], crosslinker)
in self.added_pairs_list:
1376 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1378 if (p2s[0], p1s[0], crosslinker)
in self.added_pairs_list:
1379 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1382 print(
"CrossLinkMS: added pair %d %s %d %s" % (r1, c1, r2, c2))
1383 self.added_pairs_list.append((p1s[0], p2s[0], crosslinker))
1385 rs_name =
'restraint_' + str(index)
1387 ln = impisd2.CrossLinkMSRestraint(
1389 self.crosslinker_dict[crosslinker])
1390 for i
in range(len(p1s)):
1391 ln.add_contribution(p1s[i], p2s[i])
1393 for i
in range(len(self.prots)):
1394 self.pairs.append((p1s[i], p2s[i],
1395 crosslinker, rs_name,
1396 100, 100, (r1, c1, i), (r2, c2, i),
1397 crosslinker, i, ln))
1399 self.rs.add_restraint(ln)
1402 self.rs2.add_restraint(
1403 impisd2.UniformPrior(self.sigmaglobal,
1405 self.sigmaglobal.get_upper() - 1.0,
1406 self.sigmaglobal.get_lower() + 0.1))
1407 print(
"CrossLinkMS: missing residues")
1408 for ms
in self.missing_residues:
1409 print(
"CrossLinkMS:missing " + str(ms))
1412 def add_crosslink_according_to_new_file(self, totallist):
1413 force_restraint =
False
1414 ambiguous_list = totallist[0]
1415 crosslinker = totallist[1]
1416 if (totallist[2] ==
"F"):
1417 force_restraint =
True
1426 for pair
in ambiguous_list:
1430 c1 = self.mbpnc[pair[0][0]]
1432 "CrossLinkMS: WARNING> protein name " + \
1433 pair[0][0] +
" was not defined"
1436 c2 = self.mbpnc[pair[1][0]]
1438 "CrossLinkMS: WARNING> protein name " + \
1439 pair[1][0] +
" was not defined"
1441 r1 = int(pair[0][1])
1442 r2 = int(pair[1][1])
1444 print(
'''CrossLinkMS: attempting to add restraint between
1445 residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
1452 atom_type=IMP.atom.AT_CA)
1453 p1 = (s1.get_selected_particles()[0])
1455 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1457 self.missing_residues.append((r1, c1))
1463 atom_type=IMP.atom.AT_CA)
1464 p2 = (s2.get_selected_particles()[0])
1466 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1468 self.missing_residues.append((r2, c2))
1476 atom_type=IMP.atom.AT_CA)
1477 p1 = s1.get_selected_particles()[0]
1482 atom_type=IMP.atom.AT_CA)
1483 p2 = s2.get_selected_particles()[0]
1486 if (p1, p2, crosslinker)
in self.added_pairs_list:
1487 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1489 if (p2, p1, crosslinker)
in self.added_pairs_list:
1490 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1498 print(
'''CrossLinkMS: WARNING> residue %d of chain %s and
1499 residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
1509 print(
"CrossLinkMS: added pair %d %s %d %s" % (r1, c1, r2, c2))
1511 self.added_pairs_list.append((p1, p2, crosslinker))
1514 rs_name =
'{:05}'.format(self.index % 100000)
1516 ln = impisd2.CrossLinkMSRestraint(
1518 self.crosslinker_dict[crosslinker])
1519 for i
in range(len(p1s)):
1521 ln.add_contribution(p1s[i], p2s[i])
1522 self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name,
1523 100, 100, (r1s[i], c1s[i], i), (r2s[i], c2s[i], i), crosslinker, i, ln))
1524 self.rs.add_restraint(ln)
1528 def simulate_data(self, crosslinker, weights, sensitivity_threshold=0.1,
1529 false_positive_half=0.02, elapsed_time=0.01,
1530 ratemin=100, ratemax=100):
1532 from random
import choice
1533 from random
import randrange
1534 from itertools
import combinations
1535 from numpy.random
import binomial
1537 self.weights = weights
1538 self.sensitivity_threshold = sensitivity_threshold
1539 self.false_positive_half = false_positive_half
1540 self.elapsed_time = elapsed_time
1541 self.ratemin = ratemin
1542 self.ratemax = ratemax
1546 if self.reaction_rates
is None:
1547 self.reaction_rates = {}
1552 atom_type=IMP.atom.AT_CA)
1553 self.residue_list = []
1555 for p1
in s0.get_selected_particles():
1556 (r1, c1) = tools.get_residue_index_and_chain_from_particle(p1)
1557 self.residue_list.append((r1, c1))
1558 if self.ratemin != self.ratemax:
1559 self.reaction_rates[(
1561 c1)] = randrange(self.ratemin,
1565 self.reaction_rates[(r1, c1)] = self.ratemax
1567 if self.allpairs_database
is None:
1568 self.allpairs_database = []
1571 allcomb = list(combinations(self.residue_list, 2))
1572 for ((r1, c1), (r2, c2))
in allcomb:
1577 for copy
in self.prots:
1579 chains=c1, atom_type=IMP.atom.AT_CA)
1580 p1s.append(s1.get_selected_particles()[0])
1582 chains=c2, atom_type=IMP.atom.AT_CA)
1583 p2s.append(s2.get_selected_particles()[0])
1585 ln = impisd2.CrossLinkMSRestraint(
1587 self.crosslinker_dict[crosslinker])
1588 for i
in range(len(p1s)):
1589 ln.add_contribution(p1s[i], p2s[i])
1593 reactionrate1 = self.reaction_rates[(r1, c1)]
1594 reactionrate2 = self.reaction_rates[(r2, c2)]
1595 prob = ln.get_marginal_probabilities()[i]
1596 effrate = float(reactionrate1 * reactionrate2) / \
1597 (reactionrate1 + reactionrate2)
1598 probt = self.weights[i] * \
1599 (1 - exp(-effrate * prob * elapsed_time))
1600 falsepositiveprob = exp(-probt / false_positive_half)
1601 falsepositivebool =
False
1602 falsepositive = binomial(n=1, p=falsepositiveprob)
1603 if (falsepositive == 1):
1604 falsepositivebool =
True
1605 randompair = choice(allcomb)
1606 randpair1 = randompair[0]
1607 randpair2 = randompair[1]
1609 randpair1 = (r1, c1)
1610 randpair2 = (r2, c2)
1611 if (probt > sensitivity_threshold):
1614 detectedbool =
False
1616 self.allpairs_database.append({})
1617 self.allpairs_database[-1][
"Particle1"] = p1s[i]
1618 self.allpairs_database[-1][
"Particle2"] = p2s[i]
1619 self.allpairs_database[-1][
"Distance"] = dist
1620 self.allpairs_database[-1][
"Crosslinker"] = crosslinker
1621 self.allpairs_database[-1][
"IMPRestraint"] = ln
1622 self.allpairs_database[-1][
"IMPRestraint_Probability"] = prob
1623 self.allpairs_database[-1][
"Reaction_Rate1"] = reactionrate1
1624 self.allpairs_database[-1][
"Reaction_Rate2"] = reactionrate2
1625 self.allpairs_database[-1][
"Effective_Rate"] = effrate
1626 self.allpairs_database[-1][
"CrossLink_Fraction"] = probt
1627 self.allpairs_database[
1628 -1][
"Resid1_Chainid1_Copy1"] = (r1, c1, i)
1629 self.allpairs_database[
1630 -1][
"Resid2_Chainid2_Copy2"] = (r2, c2, i)
1631 self.allpairs_database[
1632 -1][
"Is_False_Positive"] = falsepositivebool
1633 self.allpairs_database[-1][
"Identified_Pair1"] = randpair1
1634 self.allpairs_database[-1][
"Identified_Pair2"] = randpair2
1635 self.allpairs_database[-1][
"Is_Detected"] = detectedbool
1637 def set_hierarchy(self, prots):
1641 def initialize_simulated_database(self):
1643 self.allpairs_database =
None
1645 def get_number_detected_inter(self, xl_type):
1649 for el
in self.allpairs_database:
1650 if el[
"Is_Detected"]
and \
1651 ( el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1] )
and \
1652 el[
"Crosslinker"] == xl_type:
1656 def get_number_detected_inter_false_positive(self, xl_type):
1660 for el
in self.allpairs_database:
1661 if el[
"Is_Detected"]
and \
1662 ( el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1] )
and \
1663 el[
"Crosslinker"] == xl_type
and el[
"Is_False_Positive"]:
1667 def show_simulated_data(self, what="Inter"):
1669 if not self.allpairs_database
is None:
1671 for el
in self.allpairs_database:
1673 if el[
"Is_Detected"]:
1674 p1 = el[
"Identified_Pair1"]
1675 p2 = el[
"Identified_Pair2"]
1676 isf = el[
"Is_False_Positive"]
1678 el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1])
1679 cl = el[
"Crosslinker"]
1680 detectedlist.append((p1, p2, isf, cl, isinter))
1682 if el[
"Is_Detected"]
and what ==
"Detected":
1684 if el[
"Is_Detected"]
and el[
"Is_False_Positive"]
and what ==
"FalsePositive":
1686 if el[
"Is_Detected"]
and el[
"Is_False_Positive"] ==
False and what ==
"TruePositive":
1688 if el[
"Is_Detected"]
and what ==
"Intra" and \
1689 (el[
"Identified_Pair1"][1] == el[
"Identified_Pair2"][1]):
1691 if el[
"Is_Detected"]
and what ==
"Inter" and \
1692 (el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1]):
1698 print(
"Residue1: %6s, chainid1: %6s, copy1: %6d" % el[
"Resid1_Chainid1_Copy1"])
1699 print(
"Residue2: %6s, chainid2: %6s, copy2: %6d" % el[
"Resid2_Chainid2_Copy2"])
1700 keylist = list(el.keys())
1703 print(
"----", k, el[k])
1704 print(
"r1 c1 r2 c2 FP XL Inter")
1705 for d
in detectedlist:
1706 print(d[0][0], d[0][1], d[1][0], d[1][1], d[2], d[3], d[4])
1708 print(
"CrossLinkMS: Simulated data not initialized")
1711 def dump_simulated_data(self, filename="simulated_cross_link.dat"):
1713 sclf = open(filename,
"w")
1714 for el
in self.allpairs_database:
1719 def write_simulated_data(self, filename="detected_cross_link.dat"):
1721 sclf = open(filename,
"w")
1723 for el
in self.allpairs_database:
1724 if el[
"Is_Detected"]:
1726 p1 = el[
"Identified_Pair1"]
1727 p2 = el[
"Identified_Pair2"]
1728 isf = el[
"Is_False_Positive"]
1729 cl = el[
"Crosslinker"]
1743 def set_label(self, label):
1746 def add_to_model(self):
1747 self.m.add_restraint(self.rs)
1748 self.m.add_restraint(self.rs2)
1750 def get_hierarchies(self):
1754 return self.sigmaglobal
1756 def get_restraint_sets(self):
1757 return self.rs, self.rs2
1759 def get_restraint(self):
1761 tmprs.add_restraint(self.rs)
1762 tmprs.add_restraint(self.rs2)
1765 def set_output_level(self, level="low"):
1767 self.outputlevel = level
1769 def print_chimera_pseudobonds(self, filesuffix, model=0):
1770 f = open(filesuffix +
".chimera",
"w")
1772 for p
in self.pairs:
1773 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
1774 " #" + str(model) +
":" + \
1775 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
1779 def get_particles_to_sample(self):
1781 if self.sigmaissampled:
1782 ps[
"Nuisances_CrossLinkMS_Sigma_" +
1783 self.label] = ([self.sigmaglobal], self.sigmatrans)
1786 def get_output(self):
1791 score = self.rs.unprotected_evaluate(
None)
1792 score2 = self.rs2.unprotected_evaluate(
None)
1793 output[
"_TotalScore"] = str(score + score2)
1795 output[
"CrossLinkMS_Likelihood_" + self.label] = str(score)
1796 output[
"CrossLinkMS_Prior_" + self.label] = str(score2)
1797 output[
"CrossLinkMS_Sigma"] = str(self.sigmaglobal.get_scale())
1799 if self.outputlevel ==
"high":
1800 for i
in range(len(self.pairs)):
1802 p0 = self.pairs[i][0]
1803 p1 = self.pairs[i][1]
1804 crosslinker = self.pairs[i][2]
1805 ln = self.pairs[i][10]
1806 index = self.pairs[i][9]
1807 rsname = self.pairs[i][3]
1808 resid1 = self.pairs[i][6][0]
1809 chain1 = self.pairs[i][6][1]
1810 copy1 = self.pairs[i][6][2]
1811 resid2 = self.pairs[i][7][0]
1812 chain2 = self.pairs[i][7][1]
1813 copy2 = self.pairs[i][7][2]
1814 label_copy = str(rsname) +
":" + str(index) +
"-" + str(resid1) + \
1815 ":" + chain1 +
"_" +
"-" + \
1816 str(resid2) +
":" + chain2 +
"_" + crosslinker
1817 output[
"CrossLinkMS_Partial_Probability_" +
1818 label_copy] = str(ln.get_marginal_probabilities()[index])
1821 label = str(resid1) +
":" + chain1 + \
1822 "_" + str(resid2) +
":" + chain2
1824 output[
"CrossLinkMS_Score_" +
1825 str(rsname) +
"_" + crosslinker +
"_" + label] = str(ln.unprotected_evaluate(
None))
1829 output[
"CrossLinkMS_Distance_" +
1837 class BinomialXLMSRestraint(object):
1839 def __init__(self, m, prots,
1840 listofxlinkertypes=[
"BS3",
"BS2G",
"EGS"], map_between_protein_names_and_chains=
None, typeofprofile=
'pfes'):
1842 if map_between_protein_names_and_chains
is None:
1843 map_between_protein_names_and_chains = {}
1845 global impisd2, tools, exp
1846 import IMP.isd2
as impisd2
1858 self.weightmaxtrans = 0.05
1859 self.weightissampled =
False
1861 self.sigmainit = 5.0
1863 self.sigmaminnuis = 0.0
1864 self.sigmamax = 10.0
1865 self.sigmamaxnuis = 11.0
1867 self.sigmaissampled =
True
1868 self.sigmamaxtrans = 0.1
1874 self.betaissampled =
True
1875 print(
"BinomialXLMSRestraint: beta is sampled")
1877 self.betaissampled =
False
1878 print(
"BinomialXLMSRestraint: beta is NOT sampled")
1879 self.betamaxtrans = 0.01
1882 self.deltainit=0.001
1885 self.deltaissampled=False
1886 self.deltamaxtrans=0.001
1890 self.lamminnuis=0.00001
1892 self.lammaxnuis=100.0
1893 self.lamissampled=False
1894 self.lammaxtrans=0.1
1898 self.psi_dictionary = {}
1900 self.sigma = tools.SetupNuisance(self.m, self.sigmainit,
1901 self.sigmaminnuis, self.sigmamaxnuis, self.sigmaissampled).get_particle()
1903 self.beta = tools.SetupNuisance(self.m, self.betainit,
1904 self.betamin, self.betamax, self.betaissampled).get_particle()
1907 self.delta=tools.SetupNuisance(self.m,self.deltainit,
1908 self.deltamin,self.deltamax,self.deltaissampled).get_particle()
1910 self.lam=tools.SetupNuisance(self.m,self.laminit,
1911 self.lamminnuis,self.lammaxnuis,self.lamissampled).get_particle()
1913 self.weight=tools.SetupWeight(m,False).get_particle()
1915 for n in range(len(self.prots)):
1916 self.weight.add_weight()
1918 self.outputlevel =
"low"
1919 self.listofxlinkertypes = listofxlinkertypes
1921 self.reaction_rates =
None
1922 self.allpairs_database =
None
1923 self.residue_list =
None
1925 self.crosslinker_dict = self.get_crosslinker_dict(typeofprofile)
1928 self.mbpnc = map_between_protein_names_and_chains
1931 def create_psi(self, index, value):
1934 self.psiissampled =
True
1935 print(
"BinomialXLMSRestraint: psi " + str(index) +
" is sampled")
1937 self.psiinit = value
1938 self.psiissampled =
False
1939 print(
"BinomialXLMSRestraint: psi " + str(index) +
" is NOT sampled")
1940 self.psiminnuis = 0.0000001
1941 self.psimaxnuis = 0.4999999
1944 self.psitrans = 0.01
1945 self.psi = tools.SetupNuisance(self.m, self.psiinit,
1946 self.psiminnuis, self.psimaxnuis, self.psiissampled).get_particle()
1947 self.psi_dictionary[index] = (
1952 def get_psi(self, index, value):
1953 if not index
in self.psi_dictionary:
1954 self.create_psi(index, value)
1955 return self.psi_dictionary[index]
1957 def get_crosslinker_dict(self, typeofprofile):
1962 disttuple = (0.0, 200.0, 500)
1963 omegatuple = (0.01, 1000.0, 30)
1964 sigmatuple = (self.sigmamin, self.sigmamax, self.nsigma)
1966 crosslinker_dict = {}
1967 if "BS3" in self.listofxlinkertypes:
1968 crosslinker_dict[
"BS3"] = tools.get_cross_link_data(
"bs3l",
1969 "pmf_bs3l_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=
None, doff=
None, prior=1, type_of_profile=typeofprofile)
1970 if "BS2G" in self.listofxlinkertypes:
1971 crosslinker_dict[
"BS2G"] = tools.get_cross_link_data(
"bs2gl",
1972 "pmf_bs2gl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=
None, doff=
None, prior=1, type_of_profile=typeofprofile)
1973 if "EGS" in self.listofxlinkertypes:
1974 crosslinker_dict[
"EGS"] = tools.get_cross_link_data(
"egl",
1975 "pmf_egl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=
None, doff=
None, prior=1, type_of_profile=typeofprofile)
1976 if "Short" in self.listofxlinkertypes:
1978 crosslinker_dict[
"Short"] = tools.get_cross_link_data_from_length(
1983 return crosslinker_dict
1988 f = open(restraint_file)
1989 restraint_list = f.readlines()
1993 self.added_pairs_list = []
1994 self.missing_residues = []
1995 for line
in restraint_list:
1998 force_restraint =
False
2001 totallist = eval(line)
2002 self.add_crosslink_according_to_new_file(totallist)
2004 self.rs2.add_restraint(
2005 impisd2.UniformPrior(
2012 for psiindex
in self.psi_dictionary:
2013 if self.psi_dictionary[psiindex][2]:
2014 psip = self.psi_dictionary[psiindex][0]
2017 print(
"BinomialXLMSRestraint: setup 0, adding BinomialJeffreysPrior to psi particle " + str(psiindex))
2018 self.rs2.add_restraint(impisd2.BinomialJeffreysPrior(psip))
2021 self.rs2.add_restraint(
2022 impisd2.UniformPrior(
2028 def add_crosslink_according_to_new_file(self, totallist, constructor=0):
2032 force_restraint =
False
2033 ambiguous_list = totallist[0]
2034 crosslinker = totallist[1]
2035 if (totallist[2] ==
"F"):
2036 force_restraint =
True
2047 for pair
in ambiguous_list:
2051 c1 = self.mbpnc[pair[0][0]]
2053 print(
"BinomialXLMSRestraint: WARNING> protein name " + pair[0][0] +
" was not defined")
2056 c2 = self.mbpnc[pair[1][0]]
2058 print(
"BinomialXLMSRestraint: WARNING> protein name " + pair[1][0] +
" was not defined")
2061 r1 = int(pair[0][1])
2062 r2 = int(pair[1][1])
2063 psi = float(pair[2])
2065 psivalue = float(pair[3])
2069 print(
'''CrossLinkMS: attempting to add restraint between
2070 residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
2077 atom_type=IMP.atom.AT_CA)
2078 p1 = (s1.get_selected_particles()[0])
2080 print(
"BinomialXLMSRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
2082 self.missing_residues.append((r1, c1))
2088 atom_type=IMP.atom.AT_CA)
2089 p2 = (s2.get_selected_particles()[0])
2091 print(
"BinomialXLMSRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
2093 self.missing_residues.append((r2, c2))
2101 atom_type=IMP.atom.AT_CA)
2102 p1 = s1.get_selected_particles()[0]
2107 atom_type=IMP.atom.AT_CA)
2108 p2 = s2.get_selected_particles()[0]
2111 if (p1, p2, crosslinker)
in self.added_pairs_list:
2112 print(
"BinomialXLMSRestraint: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2114 if (p2, p1, crosslinker)
in self.added_pairs_list:
2115 print(
"BinomialXLMSRestraint: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2123 print(
'''BinomialXLMSRestraint: WARNING> residue %d of chain %s and
2124 residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
2134 psivalues.append(psivalue)
2136 print(
"BinomialXLMSRestraint: added pair %d %s %d %s" % (r1, c1, r2, c2))
2138 self.added_pairs_list.append((p1, p2, crosslinker))
2141 rs_name =
'{:05}'.format(self.index % 100000)
2144 print(
"BinomialXLMSRestraint: constructor 0")
2145 ln = impisd2.BinomialCrossLinkMSRestraint(
2149 self.crosslinker_dict[crosslinker])
2152 print(
"BinomialXLMSRestraint: constructor 1")
2153 ln = impisd2.BinomialCrossLinkMSRestraint(
2158 self.crosslinker_dict[crosslinker])
2160 for i
in range(len(p1s)):
2161 ln.add_contribution()
2163 psi = self.get_psi(psis[i], psivalues[i])
2165 ln.add_particle_pair(
2170 self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name,
2171 self.index, 100, (r1s[i], c1s[i], i), (r2s[i], c2s[i], i), crosslinker, i, ln))
2177 self.rs2.add_restraint(pr)
2179 self.rs.add_restraint(ln)
2181 print(
"BinomialXLMSRestraint: missing residues")
2182 for ms
in self.missing_residues:
2183 print(
"BinomialXLMSRestraint:missing " + str(ms))
2189 def set_label(self, label):
2192 def add_to_model(self):
2193 self.m.add_restraint(self.rs)
2194 self.m.add_restraint(self.rs2)
2196 def get_hierarchies(self):
2200 return self.sigmaglobal
2202 def get_restraint_sets(self):
2203 return self.rs, self.rs2
2205 def get_restraint(self):
2207 tmprs.add_restraint(self.rs)
2208 tmprs.add_restraint(self.rs2)
2211 def set_output_level(self, level="low"):
2213 self.outputlevel = level
2215 def print_chimera_pseudobonds(self, filesuffix, model=0):
2216 f = open(filesuffix +
".chimera",
"w")
2218 for p
in self.pairs:
2219 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
2220 " #" + str(model) +
":" + \
2221 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
2225 def print_chimera_pseudobonds_with_psiindexes(self, filesuffix, model=0):
2227 f = open(filesuffix +
".chimera",
"w")
2229 for p
in self.pairs:
2230 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
2231 " #" + str(model) +
":" + \
2232 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
2236 def get_output(self):
2241 score = self.rs.unprotected_evaluate(
None)
2242 score2 = self.rs2.unprotected_evaluate(
None)
2243 output[
"_TotalScore"] = str(score + score2)
2244 output[
"CrossLinkMS_Likelihood_" + self.label] = str(score)
2245 output[
"CrossLinkMS_Prior_" + self.label] = str(score2)
2247 if self.outputlevel ==
"high":
2248 for i
in range(len(self.pairs)):
2250 p0 = self.pairs[i][0]
2251 p1 = self.pairs[i][1]
2252 crosslinker = self.pairs[i][2]
2253 psiindex = self.pairs[i][4]
2254 ln = self.pairs[i][10]
2255 index = self.pairs[i][9]
2256 rsname = self.pairs[i][3]
2257 resid1 = self.pairs[i][6][0]
2258 chain1 = self.pairs[i][6][1]
2259 copy1 = self.pairs[i][6][2]
2260 resid2 = self.pairs[i][7][0]
2261 chain2 = self.pairs[i][7][1]
2262 copy2 = self.pairs[i][7][2]
2263 label_copy = str(rsname) +
":" + str(index) +
"-" + str(resid1) + \
2264 ":" + chain1 +
"_" +
"-" + \
2265 str(resid2) +
":" + chain2 +
"_" + crosslinker
2268 label = str(resid1) +
":" + chain1 + \
2269 "_" + str(resid2) +
":" + chain2
2271 output[
"CrossLinkMS_Score_" + str(rsname) +
"_" + str(index) +
"_" +
2272 crosslinker +
"_" + label] = str(ln.unprotected_evaluate(
None))
2276 output[
"CrossLinkMS_Distance_" +
2279 output[
"CrossLinkMS_Sigma_" + self.label] = str(self.sigma.get_scale())
2281 output["CrossLinkMS_Delta_"+self.label]=str(self.delta.get_scale())
2282 output["CrossLinkMS_Lambda_"+self.label]=str(self.lam.get_scale())
2284 output[
"CrossLinkMS_Beta_" + self.label] = str(self.beta.get_scale())
2285 for psiindex
in self.psi_dictionary:
2286 output[
"CrossLinkMS_Psi_" +
2287 str(psiindex) +
"_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2289 for n in range(self.weight.get_number_of_states()):
2290 output["CrossLinkMS_weights_"+str(n)+"_"+self.label]=str(self.weight.get_weight(n))
2294 def get_particles_to_sample(self):
2296 if self.sigmaissampled:
2297 ps[
"Nuisances_CrossLinkMS_Sigma_" +
2298 self.label] = ([self.sigma], self.sigmamaxtrans)
2300 if self.deltaissampled:
2301 ps["Nuisances_CrossLinkMS_Delta_"+self.label]=([self.delta],self.deltamaxtrans)
2302 if self.lamissampled:
2303 ps["Nuisances_CrossLinkMS_Lambda_"+self.label]=([self.lam],self.lammaxtrans)
2305 if self.betaissampled:
2306 ps[
"Nuisances_CrossLinkMS_Beta_" +
2307 self.label] = ([self.beta], self.betamaxtrans)
2309 for psiindex
in self.psi_dictionary:
2310 if self.psi_dictionary[psiindex][2]:
2311 ps[
"Nuisances_CrossLinkMS_Psi_" +
2312 str(psiindex) +
"_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2315 if self.weightissampled:
2316 ps["Weights_CrossLinkMS_"+self.label]=([self.weight],self.weightmaxtrans)
2323 class CrossLinkMSSimple(object):
2325 def __init__(self, prot, restraints_file, TruncatedHarmonic=True):
2326 """read crosslink restraints between two residue
2327 of different chains from an external text file
2328 syntax: part_name_1 part_name_2 distance error
2329 example: 0 1 1.0 0.1"""
2333 self.m = self.prot.get_model()
2339 if TruncatedHarmonic:
2352 addedd_pairs_list = []
2353 for line
in open(restraints_file):
2357 force_restraint =
True
2358 tokens = line.split()
2361 if (tokens[0] ==
"#"):
2367 crosslinker = tokens[4]
2370 index = int(tokens[5])
2374 if (tokens[len(tokens) - 1] ==
"F"):
2375 force_restraint =
True
2377 print(
"attempting to add restraint between residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
2388 atom_type=IMP.atom.AT_CA)
2389 p1 = (s1.get_selected_particles()[0])
2391 print(
"WARNING> residue %d of chain %s is not there" % (r1, c1))
2398 atom_type=IMP.atom.AT_CA)
2399 p2 = (s2.get_selected_particles()[0])
2401 print(
"WARNING> residue %d of chain %s is not there" % (r2, c2))
2413 print(
"attempting to add restraint between residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
2418 print(
"WARNING> residue %d of chain %s and residue %d of chain %s belong to the same rigid body" % (r1, c1, r2, c2))
2423 if (p1, p2, crosslinker)
in addedd_pairs_list:
2424 print(
"WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2426 if (p2, p1, crosslinker)
in addedd_pairs_list:
2427 print(
"WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2430 print(
"added pair %d %s %d %s" % (r1, c1, r2, c2))
2432 addedd_pairs_list.append((p1, p2, crosslinker))
2434 rs_name =
'restraint_' + str(index)
2437 ln.set_name(
"CrossLinkMSSimple_" + str(r1)
2438 +
":" + str(c1) +
"-" + str(r2) +
":" + str(c2))
2441 self.rs.add_restraint(ln)
2445 self.rs.add_restraint(pr)
2462 def set_label(self, label):
2465 def add_to_model(self):
2466 self.m.add_restraint(self.rs)
2468 def get_restraint(self):
2471 def print_chimera_pseudobonds(self, filesuffix, model=0):
2472 f = open(filesuffix +
".chimera",
"w")
2474 for p
in self.pairs:
2475 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
2476 " #" + str(model) +
":" + \
2477 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
2481 def get_output(self):
2486 score = self.rs.unprotected_evaluate(
None)
2487 output[
"_TotalScore"] = str(score)
2488 output[
"CrossLinkMSSimple_Score_" + self.label] = str(score)
2489 for i
in range(len(self.pairs)):
2491 p0 = self.pairs[i][0]
2492 p1 = self.pairs[i][1]
2493 crosslinker = self.pairs[i][2]
2494 ln = self.pairs[i][9]
2495 pr = self.pairs[i][10]
2496 resid1 = self.pairs[i][6][0]
2497 chain1 = self.pairs[i][6][1]
2498 resid2 = self.pairs[i][7][0]
2499 chain2 = self.pairs[i][7][1]
2501 label = str(resid1) +
":" + chain1 +
"_" + \
2502 str(resid2) +
":" + chain2
2503 output[
"CrossLinkMSSimple_Score_" + crosslinker +
2504 "_" + label] = str(ln.unprotected_evaluate(
None))
2505 output[
"CrossLinkMSSimple_Score_Linear_" + crosslinker +
2506 "_" + label] = str(pr.unprotected_evaluate(
None))
2509 output[
"CrossLinkMSSimple_Distance_" +
2515 def get_residue_index_and_chain_from_particle(p):
2522 def select_calpha_or_residue(
2527 SelectResidue=
False):
2531 residue_index=resid, atom_type=IMP.atom.AT_CA)
2533 ps = s.get_selected_particles()
2539 print(ObjectName +
" multiple residues selected for selection residue %s chain %s " % (resid, chain))
2543 residue_index=resid)
2544 ps = s.get_selected_particles()
2550 print(ObjectName +
" multiple residues selected for selection residue %s chain %s " % (resid, chain))
2553 print(ObjectName +
" residue %s chain %s does not exist" % (resid, chain))
2559 class SetupMembranePoreRestraint(object):
2565 selection_tuples_outside=
None,
2566 selection_tuples_membrane=
None,
2567 selection_tuples_inside=
None,
2568 center=(0.0,0.0,0.0),
2571 membrane_tickness=40.0,
2577 self.m = representation.prot.get_model()
2579 self.representation=representation
2582 self.z_tickness=z_tickness
2583 self.membrane_tickness=membrane_tickness
2589 def create_representation(self):
2592 h.set_name(
"Membrane_Pore"+self.label)
2593 self.representation.prot.add_child(h)
2595 inner=self.z_tickness/2-self.membrane_tickness/2
2596 outer=self.z_tickness/2+self.membrane_tickness/2
2600 particles_surface_1=[]
2601 particles_surface_2=[]
2605 v=self.math.pi/10.0*float(i)+self.math.pi/2.0
2606 u=2.0*self.math.pi/10.0*float(j)
2607 x=(self.radius+inner*self.math.cos(v))*self.math.cos(u)
2608 y=(self.radius+inner*self.math.cos(v))*self.math.sin(u)
2609 z=inner*self.math.sin(v)
2613 d.set_coordinates((x,y,z))
2617 surface_1.append((x,y,z))
2618 particles_surface_1.append(p2)
2620 x=(self.radius+outer*self.math.cos(v))*self.math.cos(u)
2621 y=(self.radius+outer*self.math.cos(v))*self.math.sin(u)
2622 z=outer*self.math.sin(v)
2626 d.set_coordinates((x,y,z))
2630 surface_2.append((x,y,z))
2631 particles_surface_2.append(p2)
2635 x=3.0*self.radius/10.0*float(i)-self.radius*1.5
2636 y=3.0*self.radius/10.0*float(j)-self.radius*1.5
2637 if(self.math.sqrt(x**2+y**2)<self.radius):
2640 for n,z
in enumerate([self.z_tickness/2+self.membrane_tickness/2,
2641 self.z_tickness/2-self.membrane_tickness/2,
2642 -self.z_tickness/2+self.membrane_tickness/2,
2643 -self.z_tickness/2-self.membrane_tickness/2]):
2647 d.set_coordinates((x,y,z))
2651 if n == 0
or n == 3:
2652 surface_1.append((x,y,z))
2653 particles_surface_1.append(p2)
2654 if n == 1
or n == 2:
2655 surface_2.append((x,y,z))
2656 particles_surface_2.append(p2)
2658 from scipy.spatial
import Delaunay
2659 tri = Delaunay(surface_1)
2666 """Add a line between the i-th and j-th points, if not in the list already"""
2667 if (i, j) in edges or (j, i) in edges:
2671 edge_points.append(points[ [i, j] ])
2673 # loop over triangles:
2674 # ia, ib, ic = indices of corner points of the triangle
2675 for ia, ib, ic in tri.vertices:
2684 for ia, ib, ic, id in tri.simplices:
2690 p1=particles_surface_1[e[0]]
2691 p2=particles_surface_1[e[1]]
2692 print(p1,p2,e[0],e[1])
2693 IMP.atom.Bonded.setup_particle(p1)
2694 IMP.atom.Bonded.setup_particle(p2)
2695 IMP.atom.create_bond(IMP.atom.Bonded(p1),IMP.atom.Bonded(p2),1)
2698 for i
in range(len(particles_surface_1)-1):
2699 p1=particles_surface_1[i]
2700 p2=particles_surface_1[i+1]
A filter which returns true if a container containers the Pair.
double get_volume_from_residue_type(ResidueType rt)
Return an estimate for the volume of a given residue.
A function that is harmonic over an interval.
void show_molecular_hierarchy(Hierarchy h)
Print out the molecular hierarchy.
static Atom setup_particle(Model *m, ParticleIndex pi, Atom other)
static Fragment setup_particle(Model *m, ParticleIndex pi)
Upper bound harmonic function (non-zero when feature > mean)
static XYZR setup_particle(Model *m, ParticleIndex pi)
A decorator for a particle which has bonds.
Select atoms which are selected by both selectors.
A class to store an fixed array of same-typed values.
double get_drmsd(const Vector3DsOrXYZs0 &m0, const Vector3DsOrXYZs1 &m1)
Calculate distance the root mean square deviation between two sets of 3D points.
Rotation3D get_random_rotation_3d(const Rotation3D ¢er, double distance)
Pick a rotation at random near the provided one.
Color get_rgb_color(double f)
Return the color for f from the RGB color map.
double get_ball_radius_from_volume_3d(double volume)
Return the radius of a sphere with a given volume.
ParticlesTemp get_particles(Model *m, const ParticleIndexes &ps)
static Residue setup_particle(Model *m, ParticleIndex pi, ResidueType t, int index, int insertion_code)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Return all close ordered pairs of particles taken from the two SingletonContainers.
static XYZ setup_particle(Model *m, ParticleIndex pi)
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
void read_pdb(TextInput input, int model, Hierarchy h)
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
Object used to hold a set of restraints.
static Molecule setup_particle(Model *m, ParticleIndex pi)
Classes to handle different kinds of restraints.
Store a list of ParticleIndexPairs.
Select all non-alternative ATOM records.
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
double get_volume_from_mass(double m, ProteinDensityReference ref=ALBER)
Estimate the volume of a protein from its mass.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Store a list of ParticleIndexes.
A decorator for a particle representing an atom.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
static Bonded setup_particle(Model *m, ParticleIndex pi)
void add_restraints(RMF::FileHandle fh, const Restraints &hs)
A decorator for a particle with x,y,z coordinates.
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
this class initializes a CrossLinkMS restraint and contains all useful information, such as the cross-link database, contained in self.pairs If restraint_file=None, it will proceed creating simulated data
A decorator for a residue.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Hierarchy create_simplified_along_backbone(Chain input, const IntRanges &residue_segments, bool keep_detailed=false)
Prevent a set of particles and rigid bodies from inter-penetrating.
Class to handle individual model particles.
Store info for a chain of a protein.
Applies a PairScore to a Pair.
Select all CA ATOM records.
Functionality for loading, creating, manipulating and scoring atomic structures.
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Hierarchies get_leaves(const Selection &h)
Select hierarchy particles identified by the biological name.
static RigidBody setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor ps)
Applies a PairScore to each Pair in a list.
ResidueType get_residue_type(char c)
Get the residue type from the 1-letter amino acid code.
Harmonic function (symmetric about the mean)
A decorator for a particle with x,y,z coordinates and a radius.