1 """@namespace IMP.pmi1.nonmaintained
5 from __future__
import print_function
10 "If you use any functionality from this module, please let the "
11 "IMP developers know. It will be removed in the next IMP release.")
15 def __init__(self, m):
18 self.rigid_bodies = []
19 self.floppy_bodies = []
20 self.maxtrans_rb = 2.0
21 self.maxtrans_fb = 2.0
24 def add_protein(self, name, res):
25 (firstres, lastres) = res
26 from math
import pi, cos, sin
29 nres = lastres - firstres
30 radius = (nres) * 5 / 2 / pi
32 for res
in range(firstres, lastres):
33 alpha = 2 * pi / nres * (res - firstres)
34 x = radius * cos(alpha)
35 y = radius * sin(alpha)
40 d.set_coordinates_are_optimized(
True)
42 self.hier.add_child(h)
44 def get_hierarchy(self):
47 def set_rod(self, chainname, res):
48 (firstres, lastres) = res
53 residue_indexes=list(range(
56 ps = sel.get_selected_particles()
58 self.rigid_bodies.append(rb)
60 def get_particles_to_sample(self):
65 ps[
"Floppy_Bodies_Rods"] = (self.floppy_bodies, self.maxtrans_fb)
66 ps[
"Rigid_Bodies_Rods"] = (
75 def __init__(self, m):
80 self.floppy_bodies = []
81 self.maxtrans_fb = 0.2
82 self.particle_database = {}
84 def add_bead(self, radius, label="None", color=None):
88 self.particle_database[label] = p
89 self.floppy_bodies.append(p)
93 d.set_coordinates_are_optimized(
True)
100 self.hier.add_child(p)
101 if not color
is None:
102 self.set_color(label, color)
103 return self.particle_database[label]
105 def set_color(self, label, value):
106 p = self.particle_database[label]
110 def set_floppy_bodies_max_trans(self, maxtrans):
111 self.maxtrans_fb = maxtrans
113 def get_hierarchy(self):
116 def get_bead(self, label):
117 return self.particle_database[label]
119 def set_maxtrans_fb(self, maxtrans_fb):
120 self.maxtrans_fb = maxtrans_fb
122 def get_particles_to_sample(self):
127 ps[
"Floppy_Bodies_Beads"] = (self.floppy_bodies, self.maxtrans_fb)
131 class MultipleStates(object):
133 def __init__(self, nstates, m):
134 global itertools, tools, restraints
140 self.floppy_bodies = []
141 self.rigid_bodies = []
142 for ncopy
in range(nstates):
143 self.floppy_bodies.append([])
144 self.rigid_bodies.append([])
146 self.rigid_bodies_are_sampled =
True
147 self.floppy_bodies_are_sampled =
True
150 self.prot_lowres = {}
151 self.nstates = nstates
155 self.xyzmodellist = []
157 self.maxtrans_rb = 0.15
158 self.maxrot_rb = 0.03
159 self.maxtrans_fb = 0.15
165 def set_label(self, label):
168 def set_rigid_bodies_are_sampled(self, input=True):
169 self.rigid_bodies_are_sampled = input
171 def set_floppy_bodies_are_sampled(self, input=True):
172 self.floppy_bodies_are_sampled = input
174 def get_rigid_bodies(self):
175 return self.rigid_bodies
177 def set_rigid_bodies_max_trans(self, maxtrans):
178 self.maxtrans_rb = maxtrans
180 def set_rigid_bodies_max_rot(self, maxrot):
181 self.maxrot_rb = maxrot
183 def set_floppy_bodies_max_trans(self, maxtrans):
184 self.maxtrans_fb = maxtrans
186 def get_hierarchies(self):
189 def destroy_residues(self, segments):
192 for prot
in self.prot:
193 for segment
in segments:
195 if (segment[0] == -1
or segment[1] == -1):
201 residue_indexes=list(range(segment[0],
203 for p
in s.get_selected_particles():
205 print(
"MultipleStates: one particle was not destroyed because it was a RigidMember.")
215 def add_residues_to_chains(
218 residue_type=IMP.atom.LYS):
222 for rc
in residuechainlist:
228 atom_type=IMP.atom.AT_CA)
230 print(s.get_selected_particles())
231 if len(s.get_selected_particles()) == 0:
232 for prot
in self.prot:
233 print(
"adding " + str(rc))
238 d.set_coordinates_are_optimized(
True)
247 print(tools.get_residue_index_and_chain_from_particle(a))
250 p = s.get_selected_particles()[0]
252 print(rc, s.get_selected_particles()[0])
259 atom_type=IMP.atom.AT_CA)
261 print(s.get_selected_particles())
263 def add_beads(self, segments, xyzs=None, radii=None, colors=None):
265 this method generate beads in missing portions.
266 The segments argument must be a list of selections
267 in the form [(firstres,lastres,chain)]
268 each selection will generate a bead
279 for n, s
in enumerate(segments):
284 for prot
in self.prot:
285 for prot
in self.prot:
286 cps = IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
289 if chain.get_id() == chainid:
292 rindexes = list(range(firstres, lasteres + 1))
293 f.set_residue_indexes(rindexes)
295 "Fragment_" +
'%i-%i' %
298 mass = len(rindexes) * 110.0
300 if n + 1 > len(radii):
301 mass = len(rindexes) * 110.0
303 radius = (3 * vol / math.pi) ** (1 / 3)
307 if n + 1 > len(xyzs):
316 if n + 1 <= len(colors):
320 d = IMP.atom.XYZR.setup_particle(
323 def renumber_residues(self, chainid, newfirstresiduenumber):
324 for prot
in self.prot:
325 cps = IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
328 ps = c.get_children()
331 offs = newfirstresiduenumber - ri
335 r.set_index(ri + offs)
337 def destroy_everything_but_the_residues(self, segments):
339 for prot
in self.prot:
341 for segment
in segments:
344 if (segment[0] == -1
or segment[1] == -1):
350 residue_indexes=list(range(segment[0],
352 pstokeep += s.get_selected_particles()
355 if p
not in pstokeep:
357 print(
"MultipleStates: one particle was not destroyed because it was a RigidMember.")
366 def generate_linkers_restraint_and_floppy_bodies(self, segment):
368 this methods automatically links the particles consecutively
369 according to the sequence. The restraint applied is a harmonic upper bound,
370 with a distance that is proportional to the number of residues
375 linker_restraint_objects = []
376 for ncopy, prot
in enumerate(self.prot):
377 if (segment[0] == -1
or segment[1] == -1):
383 residue_indexes=list(range(segment[0],
386 for p
in s.get_selected_particles():
388 (r, c) = tools.get_residue_index_and_chain_from_particle(p)
393 (r, c) = tools.get_residue_index_and_chain_from_particle(p)
394 p.set_name(str(r) +
":" + str(c))
395 tools.set_floppy_body(p)
396 self.floppy_bodies[ncopy].append(p)
398 residue_indexes.append((r, Floppy, c, p))
400 residue_indexes.sort()
402 pruned_residue_list = []
403 r0 = residue_indexes[0]
404 pruned_residue_list.append(r0)
408 for i
in range(1, len(residue_indexes)):
409 r = residue_indexes[i]
413 pruned_residue_list.append(r0)
414 pruned_residue_list.append(r)
416 elif r[1] != r0[1]
and r0[1] ==
False:
417 pruned_residue_list.append(r0)
418 pruned_residue_list.append(r)
420 elif r[1] == r0[1]
and r0[1]:
421 pruned_residue_list.append(r)
423 elif r[1] != r0[1]
and r[1] ==
False:
424 pruned_residue_list.append(r)
427 r0 = pruned_residue_list[0]
430 for i
in range(1, len(pruned_residue_list)):
431 r = pruned_residue_list[i]
435 linkdomaindef.append((r0[0], r[0], r[2]))
438 print(
" creating linker between atoms defined by: " + str(linkdomaindef))
440 ld = restraints.LinkDomains(prot, linkdomaindef, 1.0, 3.0)
441 ld.set_label(str(ncopy))
443 linker_restraint_objects.append(ld)
446 return linker_restraint_objects
448 def get_ref_hierarchies(self):
451 def get_number_of_states(self):
454 def get_rigid_bodies(self):
456 for rbl
in self.rigid_bodies:
461 def get_floppy_bodies(self):
463 for fbl
in self.floppy_bodies:
468 def set_rigid_bodies(self, rigid_body_list):
469 if len(self.prot) == 0:
470 print(
"MultipleStates.set_rigid_bodies: hierarchy was not initialized")
472 for ncopy, prot
in enumerate(self.prot):
474 for element
in rigid_body_list:
476 for interval
in element:
480 if (interval[0] == -1
or interval[1] == -1):
486 residue_indexes=list(range(interval[0],
488 for p
in s.get_selected_particles():
492 for key
in self.prot_lowres:
493 if (interval[0] == -1
or interval[1] == -1):
495 self.prot_lowres[key][ncopy],
499 self.prot_lowres[key][
500 ncopy], chains=interval[2],
501 residue_indexes=list(range(interval[0], interval[1] + 1)))
502 for p
in s.get_selected_particles():
508 rb.set_name(str(element))
511 print(
"MultipleStates.set_rigid_bodies: selection " + str(interval) +
" has zero elements")
512 self.rigid_bodies[ncopy] += rbl
514 def set_floppy_bodies(self, floppy_body_list):
517 if len(self.prot) == 0:
518 print(
"MultipleStates: hierarchy was not initialized")
521 for ncopy, prot
in enumerate(self.prot):
523 for element
in floppy_body_list:
525 for interval
in element:
529 if (interval[0] == -1
or interval[1] == -1):
535 residue_indexes=list(range(interval[0],
537 for p
in s.get_selected_particles():
538 (r, c) = tools.get_residue_index_and_chain_from_particle(
540 tools.set_floppy_body(p)
541 p.set_name(str(r) +
":" + str(c))
543 self.floppy_bodies[ncopy] += atoms
545 def get_particles_to_sample(self):
550 rblist = self.get_rigid_bodies()
551 fblist = self.get_floppy_bodies()
552 if self.rigid_bodies_are_sampled:
553 ps[
"Rigid_Bodies_MultipleStates"] = (
557 if self.floppy_bodies_are_sampled:
558 ps[
"Floppy_Bodies_MultipleStates"] = (fblist, self.maxtrans_fb)
561 def set_hierarchy_from_pdb(self, pdblistoflist):
562 "the input is a list of list of pdbs"
563 "one list for each copy"
565 for copy
in range(0, self.nstates):
566 prot = self.read_pdbs(pdblistoflist[copy])
567 self.prot.append(prot)
569 self.xyzmodellist.append(xyz)
571 def set_ref_hierarchy_from_pdb(self, pdblistoflist):
572 "the input is a list of list of pdbs"
573 "one list for each copy"
575 for copy
in range(0, self.nstates):
576 prot = self.read_pdbs(pdblistoflist[copy])
577 self.refprot.append(prot)
579 self.xyzreflist.append(xyz)
581 def read_pdbs(self, list_pdb_file):
582 """read pdbs from an external list file
583 create a simplified representation
584 if the pdbs are given a individual strings, it will read the
585 pdbs and give the chain name as specified in the pdb
586 If it s a tuple like (filename,chainname) it will read
587 the pdb and assign a name chainname
592 for pdb
in list_pdb_file:
600 #destroy CA atoms, for the future
601 for p in IMP.atom.get_leaves(h):
602 coor=IMP.core.XYZ(p).get_coordinates()
603 r=IMP.atom.Hierarchy(p).get_parent()
604 IMP.core.XYZ.setup_particle(r,coor)
608 cps = IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE)
611 #consolidate the chains
614 s0=IMP.atom.Selection(hier, chains=cid)
616 p=s0.get_selected_particles()[0]
617 re=IMP.atom.Residue(IMP.atom.Atom(p).get_parent()
618 ch=IMP.atom.Chain(re).get_parent())
625 if type(pdb) == tuple:
632 #destroy CA atoms, for the future
633 for p in IMP.atom.get_leaves(h):
634 coor=IMP.core.XYZ(p).get_coordinates()
635 r=IMP.atom.Hierarchy(p).get_parent()
636 IMP.core.XYZ.setup_particle(r,coor)
640 cps = IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE)
647 def recenter(self, prot):
648 "recenter the hierarchy"
650 center = IMP.algebra.get_zero_vector_3d()
656 d.set_coordinates(d.get_coordinates() - center)
657 d.set_coordinates_are_optimized(
True)
660 # bug generating code: keeping it for history
662 rb=IMP.atom.create_rigid_body(prot)
663 rbcoord=rb.get_coordinates()
664 rot=IMP.algebra.get_identity_rotation_3d()
665 tmptrans=IMP.algebra.Transformation3D(rot,rbcoord)
666 trans=tmptrans.get_inverse()
667 IMP.core.transform(rb,trans)
668 IMP.core.RigidBody.teardown_particle(rb)
669 self.m.remove_particle(rb)
672 def shuffle_configuration(self, bounding_box_length):
673 "shuffle configuration, used to restart the optimization"
674 "it only works if rigid bodies were initialized"
675 if len(self.rigid_bodies) == 0:
676 print(
"MultipleStates: rigid bodies were not initialized")
677 hbbl = bounding_box_length / 2
678 for rbl
in self.rigid_bodies:
688 rb.set_reference_frame(
691 def generate_simplified_hierarchy(self, nres):
693 self.prot_lowres[nres] = []
694 for prot
in self.prot:
702 self.prot_lowres[nres].append(sh)
704 def get_simplified_hierarchy(self, nres):
705 return self.prot_lowres[nres]
707 def calculate_drms(self):
710 if len(self.xyzmodellist) == 0:
711 print(
"MultipleStates: hierarchies were not initialized")
713 if len(self.xyzreflist) == 0:
714 print(
"MultipleStates: reference hierarchies were not initialized")
717 for i
in range(len(self.xyzreflist)):
718 for j
in range(len(self.xyzmodellist)):
721 self.xyzmodellist[j],
724 drmsdval = tools.get_drmsd(
725 self.xyzmodellist[j],
727 drmsd[
"MultipleStates_DRMSD_" +
728 str(i) +
"-Model_" + str(j)] = drmsdval
732 for assign
in itertools.permutations(list(range(len(self.xyzreflist)))):
734 for i, j
in enumerate(assign):
735 s += drmsd[
"MultipleStates_DRMSD_" +
736 str(j) +
"-Model_" + str(i)]
739 drmsd[
"MultipleStates_Total_DRMSD"] = min(min_drmsd)
742 def get_output(self):
744 if len(self.refprot) != 0:
745 drms = self.calculate_drms()
747 output[
"MultipleStates_Total_Score_" +
748 self.label] = str(self.m.evaluate(
False))
754 class LinkDomains(object):
756 def __init__(self, prot, resrangelist, kappa, length=5.0):
765 self.resrangelist = resrangelist
769 self.m = self.prot.get_model()
772 for pair
in self.resrangelist:
782 atom_type=IMP.atom.AT_CA)
783 p0 = s0.get_selected_particles()[0]
792 atom_type=IMP.atom.AT_CA)
793 p1 = s1.get_selected_particles()[0]
799 dist0 = float(pair[1] - pair[0]) * self.length
804 "LinkDomains_" + str(pair[0]) +
"-" + str(pair[1]) +
"_" + str(pair[2]))
805 self.rs.add_restraint(pr)
806 self.pairs.append((p0, p1, r0, c0, r1, c1, pr))
808 def set_label(self, label):
811 def add_to_model(self):
812 self.m.add_restraint(self.rs)
817 def get_hierarchy(self):
823 def get_residue_ranges(self):
824 return self.resrangelist
826 def get_length(self):
829 def get_restraint(self):
832 def get_restraints(self):
834 for r
in self.rs.get_restraints():
835 rlist.append(IMP.core.PairRestraint.get_from(r))
838 def print_chimera_pseudobonds(self, filesuffix, model=0):
839 f = open(filesuffix +
".chimera",
"w")
842 s =
"#" + str(model) +
":" + str(p[2]) +
"." + p[3] +
"@" + atype + \
843 " #" + str(model) +
":" + str(p[4]) +
"." + p[5] +
"@" + atype
847 def get_output(self):
850 score = self.rs.unprotected_evaluate(
None)
851 output[
"_TotalScore"] = str(score)
852 output[
"LinkDomains_" + self.label] = str(score)
853 for rst
in self.rs.get_restraints():
857 output[
"LinkDomains_" + rst.get_name() +
858 "_" + self.label] = rst.unprotected_evaluate(
None)
860 for i
in range(len(self.pairs)):
862 p0 = self.pairs[i][0]
863 p1 = self.pairs[i][1]
864 r0 = self.pairs[i][2]
865 c0 = self.pairs[i][3]
866 r1 = self.pairs[i][4]
867 c1 = self.pairs[i][5]
869 label = str(r0) +
":" + c0 +
"_" + str(r1) +
":" + c1
873 output[
"LinkDomains_Distance_" + label +
"_" +
882 class UpperBound(object):
884 def __init__(self, prot, respairs, kappa, length=5.0):
894 self.respairs = respairs
898 self.m = self.prot.get_model()
900 for pair
in self.respairs:
905 residue_index=pair[0],
906 atom_type=IMP.atom.AT_CA)
907 p0 = s0.get_selected_particles()[0]
915 residue_index=pair[2],
916 atom_type=IMP.atom.AT_CA)
917 p1 = s1.get_selected_particles()[0]
928 "UpperBound_" + str(pair[0]) +
"-" + str(pair[1]) +
"_" + str(pair[2]))
929 self.rs.add_restraint(pr)
931 def set_label(self, label):
934 def add_to_model(self):
935 self.m.add_restraint(self.rs)
937 def get_hierarchy(self):
943 def get_residue_ranges(self):
946 def get_length(self):
949 def get_restraint(self):
952 def get_output(self):
955 score = self.rs.unprotected_evaluate(
None)
956 output[
"_TotalScore"] = str(score)
957 output[
"UpperBound_" + self.label] = str(score)
963 class ExcludedVolumeResidue(object):
965 def __init__(self, prot, kappa):
970 self.m = self.prot.get_model()
972 atoms = IMP.atom.get_by_type(self.prot, IMP.atom.ATOM_TYPE)
980 lsa.add_particles(atoms)
983 self.rs.add_restraint(evr)
985 def set_label(self, label):
988 def add_excluded_particle_pairs(self, excluded_particle_pairs):
991 lpc.add_particle_pairs(excluded_particle_pairs)
993 IMP.core.ExcludedVolumeRestraint.get_from(
994 self.rs.get_restraints()[0]).add_pair_filter(icpf)
996 def add_to_model(self):
997 self.m.add_restraint(self.rs)
999 def get_hierarchy(self):
1002 def get_kappa(self):
1005 def get_restraint(self):
1008 def get_output(self):
1011 score = self.rs.unprotected_evaluate(
None)
1012 output[
"_TotalScore"] = str(score)
1013 output[
"ExcludedVolumeResidue_" + self.label] = str(score)
1019 class BipartiteExcludedVolumeResidue(object):
1021 def __init__(self, prot1, prot2, kappa):
1027 self.m = self.prot.get_model()
1029 atoms1 = IMP.atom.get_by_type(prot1, IMP.atom.ATOM_TYPE)
1031 ls1.add_particles(atoms1)
1039 atoms2 = IMP.atom.get_by_type(prot2, IMP.atom.ATOM_TYPE)
1041 ls2.add_particles(atoms2)
1056 self.rs.add_restraint(evr3)
1057 self.m.add_restraint(rs)
1059 def set_label(self, label):
1062 def add_to_model(self):
1063 self.m.add_restraint(self.rs)
1065 def get_hierarchy(self):
1066 return self.prot1, self.prot2
1068 def get_kappa(self):
1071 def get_restraint(self):
1074 def get_output(self):
1077 score = self.rs.unprotected_evaluate(
None)
1078 output[
"_TotalScore"] = str(score)
1079 output[
"BipartiteExcludedVolumeResidue_" + self.label] = str(score)
1085 class TemplateRestraint(object):
1087 def __init__(self, ps1, ps2, cutoff=6.5, kappa=1.0, forcerb=False):
1088 self.m = ps1[0].get_model()
1090 self.cutoff = cutoff
1093 self.forcerb = forcerb
1104 if dist <= self.cutoff:
1108 self.rs.add_restraint(pr)
1110 def set_label(self, label):
1113 def add_to_model(self):
1114 self.m.add_restraint(self.rs)
1116 def get_cutoff(self):
1119 def get_kappa(self):
1122 def get_restraint(self):
1125 def get_output(self):
1128 score = self.rs.unprotected_evaluate(
None)
1129 output[
"_TotalScore"] = str(score)
1130 output[
"TemplateRestraint_" + self.label] = str(score)
1136 class MarginalChi3Restraint(object):
1138 def __init__(self, part1, part2):
1139 global impisd2, tools
1140 import IMP.isd2
as impisd2
1143 self.m = part1.get_model()
1146 self.sigmamaxtrans = 0.1
1150 self.sigma = tools.SetupNuisance(
1158 for i
in range(len(self.ps1)):
1159 mc = impisd2.MarginalChi3Restraint(
1163 self.rs.add_restraint(mc)
1165 def set_label(self, label):
1168 def add_to_model(self):
1169 self.m.add_restraint(self.rs)
1171 def get_restraint(self):
1174 def get_particles_to_sample(self):
1176 ps[
"Nuisances_MarginalChi3Restraint_Sigma_" +
1177 self.label] = ([self.sigma], self.sigmamaxtrans)
1180 def get_output(self):
1183 score = self.rs.unprotected_evaluate(
None)
1184 output[
"_TotalScore"] = str(score)
1185 output[
"MarginalChi3Restraint_" + self.label] = str(score)
1186 output[
"MarginalChi3Restraint_Sigma_" +
1187 self.label] = str(self.sigma.get_scale())
1196 this class initializes a CrossLinkMS restraint and contains
1197 all useful information, such as the cross-link database, contained in self.pairs
1198 If restraint_file=None, it will proceed creating simulated data
1201 def __init__(self, prots,
1202 listofxlinkertypes=[
"BS3",
"BS2G",
"EGS"], map_between_protein_names_and_chains=
None,
1203 sigmamin=1.0, sigmamax=1.0, sigmagrid=1, sigmaissampled=
False, typeofprofile=
"gofr"):
1204 global impisd2, tools
1205 import IMP.isd2
as impisd2
1208 if map_between_protein_names_and_chains
is None:
1209 map_between_protein_names_and_chains = {}
1216 self.m = self.prots[0].get_model()
1217 self.sigmamin = sigmamin
1218 self.sigmamax = sigmamax
1219 self.sigmagrid = sigmagrid
1220 self.sigmaissampled = sigmaissampled
1222 self.sigmatrans = 0.1
1223 self.sigmaglobal = tools.SetupNuisance(self.m, self.sigmamin,
1224 self.sigmamin, self.sigmamax, self.sigmaissampled).get_particle()
1225 self.outputlevel =
"low"
1226 self.listofxlinkertypes = listofxlinkertypes
1228 self.reaction_rates =
None
1229 self.allpairs_database =
None
1230 self.residue_list =
None
1232 self.crosslinker_dict = self.get_crosslinker_dict(typeofprofile)
1235 self.mbpnc = map_between_protein_names_and_chains
1239 def get_crosslinker_dict(self, typeofprofile="gofr"):
1244 disttuple = (0.0, 200.0, 1000)
1245 omegatuple = (1.0, 1000.0, 30)
1246 sigmatuple = (self.sigmamin, self.sigmamax, self.sigmagrid)
1248 crosslinker_dict = {}
1249 if "BS3" in self.listofxlinkertypes:
1250 crosslinker_dict[
"BS3"] = tools.get_cross_link_data(
"bs3l",
1251 "pmf_bs3l_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1252 don=
None, doff=
None, prior=0, type_of_profile=typeofprofile)
1253 if "BS2G" in self.listofxlinkertypes:
1254 crosslinker_dict[
"BS2G"] = tools.get_cross_link_data(
"bs2gl",
1255 "pmf_bs2gl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1256 don=
None, doff=
None, prior=0, type_of_profile=typeofprofile)
1257 if "EGS" in self.listofxlinkertypes:
1258 crosslinker_dict[
"EGS"] = tools.get_cross_link_data(
"egl",
1259 "pmf_egl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1260 don=
None, doff=
None, prior=0, type_of_profile=typeofprofile)
1261 if "Short" in self.listofxlinkertypes:
1263 crosslinker_dict[
"Short"] = tools.get_cross_link_data_from_length(
1268 return crosslinker_dict
1274 if restraint_file
is None:
1276 restraint_list = self.allpairs_database
1280 f = open(restraint_file)
1281 restraint_list = f.readlines()
1285 self.added_pairs_list = []
1286 self.missing_residues = []
1287 for line
in restraint_list:
1290 force_restraint =
False
1292 if restraint_file
is None:
1293 if line[
"Is_Detected"]:
1294 crosslinker = line[
"Crosslinker"]
1295 (r1, c1) = line[
"Identified_Pair1"]
1296 (r2, c2) = line[
"Identified_Pair2"]
1302 tokens = line.split()
1304 if (tokens[0] ==
"#"):
1310 crosslinker = tokens[4]
1313 self.index = int(tokens[5])
1317 if (tokens[len(tokens) - 1] ==
"F"):
1318 force_restraint =
True
1322 totallist = eval(line)
1323 self.add_crosslink_according_to_new_file(totallist)
1327 print(
'''CrossLinkMS: attempting to add restraint between
1328 residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
1338 atom_type=IMP.atom.AT_CA)
1339 p1 = (s1.get_selected_particles()[0])
1341 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1348 atom_type=IMP.atom.AT_CA)
1349 p2 = (s2.get_selected_particles()[0])
1351 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1354 for copy
in self.prots:
1359 atom_type=IMP.atom.AT_CA)
1360 p1s.append(s1.get_selected_particles()[0])
1365 atom_type=IMP.atom.AT_CA)
1366 p2s.append(s2.get_selected_particles()[0])
1373 print(
'''CrossLinkMS: WARNING> residue %d of chain %s and
1374 residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
1379 if (p1s[0], p2s[0], crosslinker)
in self.added_pairs_list:
1380 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1382 if (p2s[0], p1s[0], crosslinker)
in self.added_pairs_list:
1383 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1386 print(
"CrossLinkMS: added pair %d %s %d %s" % (r1, c1, r2, c2))
1387 self.added_pairs_list.append((p1s[0], p2s[0], crosslinker))
1389 rs_name =
'restraint_' + str(index)
1391 ln = impisd2.CrossLinkMSRestraint(
1393 self.crosslinker_dict[crosslinker])
1394 for i
in range(len(p1s)):
1395 ln.add_contribution(p1s[i], p2s[i])
1397 for i
in range(len(self.prots)):
1398 self.pairs.append((p1s[i], p2s[i],
1399 crosslinker, rs_name,
1400 100, 100, (r1, c1, i), (r2, c2, i),
1401 crosslinker, i, ln))
1403 self.rs.add_restraint(ln)
1406 self.rs2.add_restraint(
1407 impisd2.UniformPrior(self.sigmaglobal,
1409 self.sigmaglobal.get_upper() - 1.0,
1410 self.sigmaglobal.get_lower() + 0.1))
1411 print(
"CrossLinkMS: missing residues")
1412 for ms
in self.missing_residues:
1413 print(
"CrossLinkMS:missing " + str(ms))
1416 def add_crosslink_according_to_new_file(self, totallist):
1417 force_restraint =
False
1418 ambiguous_list = totallist[0]
1419 crosslinker = totallist[1]
1420 if (totallist[2] ==
"F"):
1421 force_restraint =
True
1430 for pair
in ambiguous_list:
1434 c1 = self.mbpnc[pair[0][0]]
1436 "CrossLinkMS: WARNING> protein name " + \
1437 pair[0][0] +
" was not defined"
1440 c2 = self.mbpnc[pair[1][0]]
1442 "CrossLinkMS: WARNING> protein name " + \
1443 pair[1][0] +
" was not defined"
1445 r1 = int(pair[0][1])
1446 r2 = int(pair[1][1])
1448 print(
'''CrossLinkMS: attempting to add restraint between
1449 residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
1456 atom_type=IMP.atom.AT_CA)
1457 p1 = (s1.get_selected_particles()[0])
1459 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1461 self.missing_residues.append((r1, c1))
1467 atom_type=IMP.atom.AT_CA)
1468 p2 = (s2.get_selected_particles()[0])
1470 print(
"CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1472 self.missing_residues.append((r2, c2))
1480 atom_type=IMP.atom.AT_CA)
1481 p1 = s1.get_selected_particles()[0]
1486 atom_type=IMP.atom.AT_CA)
1487 p2 = s2.get_selected_particles()[0]
1490 if (p1, p2, crosslinker)
in self.added_pairs_list:
1491 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1493 if (p2, p1, crosslinker)
in self.added_pairs_list:
1494 print(
"CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1502 print(
'''CrossLinkMS: WARNING> residue %d of chain %s and
1503 residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
1513 print(
"CrossLinkMS: added pair %d %s %d %s" % (r1, c1, r2, c2))
1515 self.added_pairs_list.append((p1, p2, crosslinker))
1518 rs_name =
'{:05}'.format(self.index % 100000)
1520 ln = impisd2.CrossLinkMSRestraint(
1522 self.crosslinker_dict[crosslinker])
1523 for i
in range(len(p1s)):
1525 ln.add_contribution(p1s[i], p2s[i])
1526 self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name,
1527 100, 100, (r1s[i], c1s[i], i), (r2s[i], c2s[i], i), crosslinker, i, ln))
1528 self.rs.add_restraint(ln)
1532 def simulate_data(self, crosslinker, weights, sensitivity_threshold=0.1,
1533 false_positive_half=0.02, elapsed_time=0.01,
1534 ratemin=100, ratemax=100):
1536 from random
import choice
1537 from random
import randrange
1538 from itertools
import combinations
1539 from numpy.random
import binomial
1541 self.weights = weights
1542 self.sensitivity_threshold = sensitivity_threshold
1543 self.false_positive_half = false_positive_half
1544 self.elapsed_time = elapsed_time
1545 self.ratemin = ratemin
1546 self.ratemax = ratemax
1550 if self.reaction_rates
is None:
1551 self.reaction_rates = {}
1556 atom_type=IMP.atom.AT_CA)
1557 self.residue_list = []
1559 for p1
in s0.get_selected_particles():
1560 (r1, c1) = tools.get_residue_index_and_chain_from_particle(p1)
1561 self.residue_list.append((r1, c1))
1562 if self.ratemin != self.ratemax:
1563 self.reaction_rates[(
1565 c1)] = randrange(self.ratemin,
1569 self.reaction_rates[(r1, c1)] = self.ratemax
1571 if self.allpairs_database
is None:
1572 self.allpairs_database = []
1575 allcomb = list(combinations(self.residue_list, 2))
1576 for ((r1, c1), (r2, c2))
in allcomb:
1581 for copy
in self.prots:
1583 chains=c1, atom_type=IMP.atom.AT_CA)
1584 p1s.append(s1.get_selected_particles()[0])
1586 chains=c2, atom_type=IMP.atom.AT_CA)
1587 p2s.append(s2.get_selected_particles()[0])
1589 ln = impisd2.CrossLinkMSRestraint(
1591 self.crosslinker_dict[crosslinker])
1592 for i
in range(len(p1s)):
1593 ln.add_contribution(p1s[i], p2s[i])
1597 reactionrate1 = self.reaction_rates[(r1, c1)]
1598 reactionrate2 = self.reaction_rates[(r2, c2)]
1599 prob = ln.get_marginal_probabilities()[i]
1600 effrate = float(reactionrate1 * reactionrate2) / \
1601 (reactionrate1 + reactionrate2)
1602 probt = self.weights[i] * \
1603 (1 - exp(-effrate * prob * elapsed_time))
1604 falsepositiveprob = exp(-probt / false_positive_half)
1605 falsepositivebool =
False
1606 falsepositive = binomial(n=1, p=falsepositiveprob)
1607 if (falsepositive == 1):
1608 falsepositivebool =
True
1609 randompair = choice(allcomb)
1610 randpair1 = randompair[0]
1611 randpair2 = randompair[1]
1613 randpair1 = (r1, c1)
1614 randpair2 = (r2, c2)
1615 if (probt > sensitivity_threshold):
1618 detectedbool =
False
1620 self.allpairs_database.append({})
1621 self.allpairs_database[-1][
"Particle1"] = p1s[i]
1622 self.allpairs_database[-1][
"Particle2"] = p2s[i]
1623 self.allpairs_database[-1][
"Distance"] = dist
1624 self.allpairs_database[-1][
"Crosslinker"] = crosslinker
1625 self.allpairs_database[-1][
"IMPRestraint"] = ln
1626 self.allpairs_database[-1][
"IMPRestraint_Probability"] = prob
1627 self.allpairs_database[-1][
"Reaction_Rate1"] = reactionrate1
1628 self.allpairs_database[-1][
"Reaction_Rate2"] = reactionrate2
1629 self.allpairs_database[-1][
"Effective_Rate"] = effrate
1630 self.allpairs_database[-1][
"CrossLink_Fraction"] = probt
1631 self.allpairs_database[
1632 -1][
"Resid1_Chainid1_Copy1"] = (r1, c1, i)
1633 self.allpairs_database[
1634 -1][
"Resid2_Chainid2_Copy2"] = (r2, c2, i)
1635 self.allpairs_database[
1636 -1][
"Is_False_Positive"] = falsepositivebool
1637 self.allpairs_database[-1][
"Identified_Pair1"] = randpair1
1638 self.allpairs_database[-1][
"Identified_Pair2"] = randpair2
1639 self.allpairs_database[-1][
"Is_Detected"] = detectedbool
1641 def set_hierarchy(self, prots):
1645 def initialize_simulated_database(self):
1647 self.allpairs_database =
None
1649 def get_number_detected_inter(self, xl_type):
1653 for el
in self.allpairs_database:
1654 if el[
"Is_Detected"]
and \
1655 ( el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1] )
and \
1656 el[
"Crosslinker"] == xl_type:
1660 def get_number_detected_inter_false_positive(self, xl_type):
1664 for el
in self.allpairs_database:
1665 if el[
"Is_Detected"]
and \
1666 ( el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1] )
and \
1667 el[
"Crosslinker"] == xl_type
and el[
"Is_False_Positive"]:
1671 def show_simulated_data(self, what="Inter"):
1673 if not self.allpairs_database
is None:
1675 for el
in self.allpairs_database:
1677 if el[
"Is_Detected"]:
1678 p1 = el[
"Identified_Pair1"]
1679 p2 = el[
"Identified_Pair2"]
1680 isf = el[
"Is_False_Positive"]
1682 el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1])
1683 cl = el[
"Crosslinker"]
1684 detectedlist.append((p1, p2, isf, cl, isinter))
1686 if el[
"Is_Detected"]
and what ==
"Detected":
1688 if el[
"Is_Detected"]
and el[
"Is_False_Positive"]
and what ==
"FalsePositive":
1690 if el[
"Is_Detected"]
and el[
"Is_False_Positive"] ==
False and what ==
"TruePositive":
1692 if el[
"Is_Detected"]
and what ==
"Intra" and \
1693 (el[
"Identified_Pair1"][1] == el[
"Identified_Pair2"][1]):
1695 if el[
"Is_Detected"]
and what ==
"Inter" and \
1696 (el[
"Identified_Pair1"][1] != el[
"Identified_Pair2"][1]):
1702 print(
"Residue1: %6s, chainid1: %6s, copy1: %6d" % el[
"Resid1_Chainid1_Copy1"])
1703 print(
"Residue2: %6s, chainid2: %6s, copy2: %6d" % el[
"Resid2_Chainid2_Copy2"])
1704 keylist = list(el.keys())
1707 print(
"----", k, el[k])
1708 print(
"r1 c1 r2 c2 FP XL Inter")
1709 for d
in detectedlist:
1710 print(d[0][0], d[0][1], d[1][0], d[1][1], d[2], d[3], d[4])
1712 print(
"CrossLinkMS: Simulated data not initialized")
1715 def dump_simulated_data(self, filename="simulated_cross_link.dat"):
1717 sclf = open(filename,
"w")
1718 for el
in self.allpairs_database:
1723 def write_simulated_data(self, filename="detected_cross_link.dat"):
1725 sclf = open(filename,
"w")
1727 for el
in self.allpairs_database:
1728 if el[
"Is_Detected"]:
1730 p1 = el[
"Identified_Pair1"]
1731 p2 = el[
"Identified_Pair2"]
1732 isf = el[
"Is_False_Positive"]
1733 cl = el[
"Crosslinker"]
1747 def set_label(self, label):
1750 def add_to_model(self):
1751 self.m.add_restraint(self.rs)
1752 self.m.add_restraint(self.rs2)
1754 def get_hierarchies(self):
1758 return self.sigmaglobal
1760 def get_restraint_sets(self):
1761 return self.rs, self.rs2
1763 def get_restraint(self):
1765 tmprs.add_restraint(self.rs)
1766 tmprs.add_restraint(self.rs2)
1769 def set_output_level(self, level="low"):
1771 self.outputlevel = level
1773 def print_chimera_pseudobonds(self, filesuffix, model=0):
1774 f = open(filesuffix +
".chimera",
"w")
1776 for p
in self.pairs:
1777 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
1778 " #" + str(model) +
":" + \
1779 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
1783 def get_particles_to_sample(self):
1785 if self.sigmaissampled:
1786 ps[
"Nuisances_CrossLinkMS_Sigma_" +
1787 self.label] = ([self.sigmaglobal], self.sigmatrans)
1790 def get_output(self):
1795 score = self.rs.unprotected_evaluate(
None)
1796 score2 = self.rs2.unprotected_evaluate(
None)
1797 output[
"_TotalScore"] = str(score + score2)
1799 output[
"CrossLinkMS_Likelihood_" + self.label] = str(score)
1800 output[
"CrossLinkMS_Prior_" + self.label] = str(score2)
1801 output[
"CrossLinkMS_Sigma"] = str(self.sigmaglobal.get_scale())
1803 if self.outputlevel ==
"high":
1804 for i
in range(len(self.pairs)):
1806 p0 = self.pairs[i][0]
1807 p1 = self.pairs[i][1]
1808 crosslinker = self.pairs[i][2]
1809 ln = self.pairs[i][10]
1810 index = self.pairs[i][9]
1811 rsname = self.pairs[i][3]
1812 resid1 = self.pairs[i][6][0]
1813 chain1 = self.pairs[i][6][1]
1814 copy1 = self.pairs[i][6][2]
1815 resid2 = self.pairs[i][7][0]
1816 chain2 = self.pairs[i][7][1]
1817 copy2 = self.pairs[i][7][2]
1818 label_copy = str(rsname) +
":" + str(index) +
"-" + str(resid1) + \
1819 ":" + chain1 +
"_" +
"-" + \
1820 str(resid2) +
":" + chain2 +
"_" + crosslinker
1821 output[
"CrossLinkMS_Partial_Probability_" +
1822 label_copy] = str(ln.get_marginal_probabilities()[index])
1825 label = str(resid1) +
":" + chain1 + \
1826 "_" + str(resid2) +
":" + chain2
1828 output[
"CrossLinkMS_Score_" +
1829 str(rsname) +
"_" + crosslinker +
"_" + label] = str(ln.unprotected_evaluate(
None))
1833 output[
"CrossLinkMS_Distance_" +
1841 class BinomialXLMSRestraint(object):
1843 def __init__(self, m, prots,
1844 listofxlinkertypes=[
"BS3",
"BS2G",
"EGS"], map_between_protein_names_and_chains=
None, typeofprofile=
'pfes'):
1846 if map_between_protein_names_and_chains
is None:
1847 map_between_protein_names_and_chains = {}
1849 global impisd2, tools, exp
1850 import IMP.isd2
as impisd2
1862 self.weightmaxtrans = 0.05
1863 self.weightissampled =
False
1865 self.sigmainit = 5.0
1867 self.sigmaminnuis = 0.0
1868 self.sigmamax = 10.0
1869 self.sigmamaxnuis = 11.0
1871 self.sigmaissampled =
True
1872 self.sigmamaxtrans = 0.1
1878 self.betaissampled =
True
1879 print(
"BinomialXLMSRestraint: beta is sampled")
1881 self.betaissampled =
False
1882 print(
"BinomialXLMSRestraint: beta is NOT sampled")
1883 self.betamaxtrans = 0.01
1886 self.deltainit=0.001
1889 self.deltaissampled=False
1890 self.deltamaxtrans=0.001
1894 self.lamminnuis=0.00001
1896 self.lammaxnuis=100.0
1897 self.lamissampled=False
1898 self.lammaxtrans=0.1
1902 self.psi_dictionary = {}
1904 self.sigma = tools.SetupNuisance(self.m, self.sigmainit,
1905 self.sigmaminnuis, self.sigmamaxnuis, self.sigmaissampled).get_particle()
1907 self.beta = tools.SetupNuisance(self.m, self.betainit,
1908 self.betamin, self.betamax, self.betaissampled).get_particle()
1911 self.delta=tools.SetupNuisance(self.m,self.deltainit,
1912 self.deltamin,self.deltamax,self.deltaissampled).get_particle()
1914 self.lam=tools.SetupNuisance(self.m,self.laminit,
1915 self.lamminnuis,self.lammaxnuis,self.lamissampled).get_particle()
1917 self.weight=tools.SetupWeight(m,False).get_particle()
1919 for n in range(len(self.prots)):
1920 self.weight.add_weight()
1922 self.outputlevel =
"low"
1923 self.listofxlinkertypes = listofxlinkertypes
1925 self.reaction_rates =
None
1926 self.allpairs_database =
None
1927 self.residue_list =
None
1929 self.crosslinker_dict = self.get_crosslinker_dict(typeofprofile)
1932 self.mbpnc = map_between_protein_names_and_chains
1935 def create_psi(self, index, value):
1938 self.psiissampled =
True
1939 print(
"BinomialXLMSRestraint: psi " + str(index) +
" is sampled")
1941 self.psiinit = value
1942 self.psiissampled =
False
1943 print(
"BinomialXLMSRestraint: psi " + str(index) +
" is NOT sampled")
1944 self.psiminnuis = 0.0000001
1945 self.psimaxnuis = 0.4999999
1948 self.psitrans = 0.01
1949 self.psi = tools.SetupNuisance(self.m, self.psiinit,
1950 self.psiminnuis, self.psimaxnuis, self.psiissampled).get_particle()
1951 self.psi_dictionary[index] = (
1956 def get_psi(self, index, value):
1957 if not index
in self.psi_dictionary:
1958 self.create_psi(index, value)
1959 return self.psi_dictionary[index]
1961 def get_crosslinker_dict(self, typeofprofile):
1966 disttuple = (0.0, 200.0, 500)
1967 omegatuple = (0.01, 1000.0, 30)
1968 sigmatuple = (self.sigmamin, self.sigmamax, self.nsigma)
1970 crosslinker_dict = {}
1971 if "BS3" in self.listofxlinkertypes:
1972 crosslinker_dict[
"BS3"] = tools.get_cross_link_data(
"bs3l",
1973 "pmf_bs3l_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=
None, doff=
None, prior=1, type_of_profile=typeofprofile)
1974 if "BS2G" in self.listofxlinkertypes:
1975 crosslinker_dict[
"BS2G"] = tools.get_cross_link_data(
"bs2gl",
1976 "pmf_bs2gl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=
None, doff=
None, prior=1, type_of_profile=typeofprofile)
1977 if "EGS" in self.listofxlinkertypes:
1978 crosslinker_dict[
"EGS"] = tools.get_cross_link_data(
"egl",
1979 "pmf_egl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=
None, doff=
None, prior=1, type_of_profile=typeofprofile)
1980 if "Short" in self.listofxlinkertypes:
1982 crosslinker_dict[
"Short"] = tools.get_cross_link_data_from_length(
1987 return crosslinker_dict
1992 f = open(restraint_file)
1993 restraint_list = f.readlines()
1997 self.added_pairs_list = []
1998 self.missing_residues = []
1999 for line
in restraint_list:
2002 force_restraint =
False
2005 totallist = eval(line)
2006 self.add_crosslink_according_to_new_file(totallist)
2008 self.rs2.add_restraint(
2009 impisd2.UniformPrior(
2016 for psiindex
in self.psi_dictionary:
2017 if self.psi_dictionary[psiindex][2]:
2018 psip = self.psi_dictionary[psiindex][0]
2021 print(
"BinomialXLMSRestraint: setup 0, adding BinomialJeffreysPrior to psi particle " + str(psiindex))
2022 self.rs2.add_restraint(impisd2.BinomialJeffreysPrior(psip))
2025 self.rs2.add_restraint(
2026 impisd2.UniformPrior(
2032 def add_crosslink_according_to_new_file(self, totallist, constructor=0):
2036 force_restraint =
False
2037 ambiguous_list = totallist[0]
2038 crosslinker = totallist[1]
2039 if (totallist[2] ==
"F"):
2040 force_restraint =
True
2051 for pair
in ambiguous_list:
2055 c1 = self.mbpnc[pair[0][0]]
2057 print(
"BinomialXLMSRestraint: WARNING> protein name " + pair[0][0] +
" was not defined")
2060 c2 = self.mbpnc[pair[1][0]]
2062 print(
"BinomialXLMSRestraint: WARNING> protein name " + pair[1][0] +
" was not defined")
2065 r1 = int(pair[0][1])
2066 r2 = int(pair[1][1])
2067 psi = float(pair[2])
2069 psivalue = float(pair[3])
2073 print(
'''CrossLinkMS: attempting to add restraint between
2074 residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
2081 atom_type=IMP.atom.AT_CA)
2082 p1 = (s1.get_selected_particles()[0])
2084 print(
"BinomialXLMSRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
2086 self.missing_residues.append((r1, c1))
2092 atom_type=IMP.atom.AT_CA)
2093 p2 = (s2.get_selected_particles()[0])
2095 print(
"BinomialXLMSRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
2097 self.missing_residues.append((r2, c2))
2105 atom_type=IMP.atom.AT_CA)
2106 p1 = s1.get_selected_particles()[0]
2111 atom_type=IMP.atom.AT_CA)
2112 p2 = s2.get_selected_particles()[0]
2115 if (p1, p2, crosslinker)
in self.added_pairs_list:
2116 print(
"BinomialXLMSRestraint: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2118 if (p2, p1, crosslinker)
in self.added_pairs_list:
2119 print(
"BinomialXLMSRestraint: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2127 print(
'''BinomialXLMSRestraint: WARNING> residue %d of chain %s and
2128 residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
2138 psivalues.append(psivalue)
2140 print(
"BinomialXLMSRestraint: added pair %d %s %d %s" % (r1, c1, r2, c2))
2142 self.added_pairs_list.append((p1, p2, crosslinker))
2145 rs_name =
'{:05}'.format(self.index % 100000)
2148 print(
"BinomialXLMSRestraint: constructor 0")
2149 ln = impisd2.BinomialCrossLinkMSRestraint(
2153 self.crosslinker_dict[crosslinker])
2156 print(
"BinomialXLMSRestraint: constructor 1")
2157 ln = impisd2.BinomialCrossLinkMSRestraint(
2162 self.crosslinker_dict[crosslinker])
2164 for i
in range(len(p1s)):
2165 ln.add_contribution()
2167 psi = self.get_psi(psis[i], psivalues[i])
2169 ln.add_particle_pair(
2174 self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name,
2175 self.index, 100, (r1s[i], c1s[i], i), (r2s[i], c2s[i], i), crosslinker, i, ln))
2181 self.rs2.add_restraint(pr)
2183 self.rs.add_restraint(ln)
2185 print(
"BinomialXLMSRestraint: missing residues")
2186 for ms
in self.missing_residues:
2187 print(
"BinomialXLMSRestraint:missing " + str(ms))
2193 def set_label(self, label):
2196 def add_to_model(self):
2197 self.m.add_restraint(self.rs)
2198 self.m.add_restraint(self.rs2)
2200 def get_hierarchies(self):
2204 return self.sigmaglobal
2206 def get_restraint_sets(self):
2207 return self.rs, self.rs2
2209 def get_restraint(self):
2211 tmprs.add_restraint(self.rs)
2212 tmprs.add_restraint(self.rs2)
2215 def set_output_level(self, level="low"):
2217 self.outputlevel = level
2219 def print_chimera_pseudobonds(self, filesuffix, model=0):
2220 f = open(filesuffix +
".chimera",
"w")
2222 for p
in self.pairs:
2223 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
2224 " #" + str(model) +
":" + \
2225 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
2229 def print_chimera_pseudobonds_with_psiindexes(self, filesuffix, model=0):
2231 f = open(filesuffix +
".chimera",
"w")
2233 for p
in self.pairs:
2234 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
2235 " #" + str(model) +
":" + \
2236 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
2240 def get_output(self):
2245 score = self.rs.unprotected_evaluate(
None)
2246 score2 = self.rs2.unprotected_evaluate(
None)
2247 output[
"_TotalScore"] = str(score + score2)
2248 output[
"CrossLinkMS_Likelihood_" + self.label] = str(score)
2249 output[
"CrossLinkMS_Prior_" + self.label] = str(score2)
2251 if self.outputlevel ==
"high":
2252 for i
in range(len(self.pairs)):
2254 p0 = self.pairs[i][0]
2255 p1 = self.pairs[i][1]
2256 crosslinker = self.pairs[i][2]
2257 psiindex = self.pairs[i][4]
2258 ln = self.pairs[i][10]
2259 index = self.pairs[i][9]
2260 rsname = self.pairs[i][3]
2261 resid1 = self.pairs[i][6][0]
2262 chain1 = self.pairs[i][6][1]
2263 copy1 = self.pairs[i][6][2]
2264 resid2 = self.pairs[i][7][0]
2265 chain2 = self.pairs[i][7][1]
2266 copy2 = self.pairs[i][7][2]
2267 label_copy = str(rsname) +
":" + str(index) +
"-" + str(resid1) + \
2268 ":" + chain1 +
"_" +
"-" + \
2269 str(resid2) +
":" + chain2 +
"_" + crosslinker
2272 label = str(resid1) +
":" + chain1 + \
2273 "_" + str(resid2) +
":" + chain2
2275 output[
"CrossLinkMS_Score_" + str(rsname) +
"_" + str(index) +
"_" +
2276 crosslinker +
"_" + label] = str(ln.unprotected_evaluate(
None))
2280 output[
"CrossLinkMS_Distance_" +
2283 output[
"CrossLinkMS_Sigma_" + self.label] = str(self.sigma.get_scale())
2285 output["CrossLinkMS_Delta_"+self.label]=str(self.delta.get_scale())
2286 output["CrossLinkMS_Lambda_"+self.label]=str(self.lam.get_scale())
2288 output[
"CrossLinkMS_Beta_" + self.label] = str(self.beta.get_scale())
2289 for psiindex
in self.psi_dictionary:
2290 output[
"CrossLinkMS_Psi_" +
2291 str(psiindex) +
"_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2293 for n in range(self.weight.get_number_of_weights()):
2294 output["CrossLinkMS_weights_"+str(n)+"_"+self.label]=str(self.weight.get_weight(n))
2298 def get_particles_to_sample(self):
2300 if self.sigmaissampled:
2301 ps[
"Nuisances_CrossLinkMS_Sigma_" +
2302 self.label] = ([self.sigma], self.sigmamaxtrans)
2304 if self.deltaissampled:
2305 ps["Nuisances_CrossLinkMS_Delta_"+self.label]=([self.delta],self.deltamaxtrans)
2306 if self.lamissampled:
2307 ps["Nuisances_CrossLinkMS_Lambda_"+self.label]=([self.lam],self.lammaxtrans)
2309 if self.betaissampled:
2310 ps[
"Nuisances_CrossLinkMS_Beta_" +
2311 self.label] = ([self.beta], self.betamaxtrans)
2313 for psiindex
in self.psi_dictionary:
2314 if self.psi_dictionary[psiindex][2]:
2315 ps[
"Nuisances_CrossLinkMS_Psi_" +
2316 str(psiindex) +
"_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2319 if self.weightissampled:
2320 ps["Weights_CrossLinkMS_"+self.label]=([self.weight],self.weightmaxtrans)
2327 class CrossLinkMSSimple(object):
2329 def __init__(self, prot, restraints_file, TruncatedHarmonic=True):
2330 """read crosslink restraints between two residue
2331 of different chains from an external text file
2332 syntax: part_name_1 part_name_2 distance error
2333 example: 0 1 1.0 0.1"""
2337 self.m = self.prot.get_model()
2343 if TruncatedHarmonic:
2356 addedd_pairs_list = []
2357 for line
in open(restraints_file):
2361 force_restraint =
True
2362 tokens = line.split()
2365 if (tokens[0] ==
"#"):
2371 crosslinker = tokens[4]
2374 index = int(tokens[5])
2378 if (tokens[len(tokens) - 1] ==
"F"):
2379 force_restraint =
True
2381 print(
"attempting to add restraint between residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
2392 atom_type=IMP.atom.AT_CA)
2393 p1 = (s1.get_selected_particles()[0])
2395 print(
"WARNING> residue %d of chain %s is not there" % (r1, c1))
2402 atom_type=IMP.atom.AT_CA)
2403 p2 = (s2.get_selected_particles()[0])
2405 print(
"WARNING> residue %d of chain %s is not there" % (r2, c2))
2417 print(
"attempting to add restraint between residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
2422 print(
"WARNING> residue %d of chain %s and residue %d of chain %s belong to the same rigid body" % (r1, c1, r2, c2))
2427 if (p1, p2, crosslinker)
in addedd_pairs_list:
2428 print(
"WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2430 if (p2, p1, crosslinker)
in addedd_pairs_list:
2431 print(
"WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2434 print(
"added pair %d %s %d %s" % (r1, c1, r2, c2))
2436 addedd_pairs_list.append((p1, p2, crosslinker))
2438 rs_name =
'restraint_' + str(index)
2441 ln.set_name(
"CrossLinkMSSimple_" + str(r1)
2442 +
":" + str(c1) +
"-" + str(r2) +
":" + str(c2))
2445 self.rs.add_restraint(ln)
2449 self.rs.add_restraint(pr)
2466 def set_label(self, label):
2469 def add_to_model(self):
2470 self.m.add_restraint(self.rs)
2472 def get_restraint(self):
2475 def print_chimera_pseudobonds(self, filesuffix, model=0):
2476 f = open(filesuffix +
".chimera",
"w")
2478 for p
in self.pairs:
2479 s =
"#" + str(model) +
":" + str(p[6][0]) +
"." + p[6][1] +
"@" + atype + \
2480 " #" + str(model) +
":" + \
2481 str(p[7][0]) +
"." + p[7][1] +
"@" + atype
2485 def get_output(self):
2490 score = self.rs.unprotected_evaluate(
None)
2491 output[
"_TotalScore"] = str(score)
2492 output[
"CrossLinkMSSimple_Score_" + self.label] = str(score)
2493 for i
in range(len(self.pairs)):
2495 p0 = self.pairs[i][0]
2496 p1 = self.pairs[i][1]
2497 crosslinker = self.pairs[i][2]
2498 ln = self.pairs[i][9]
2499 pr = self.pairs[i][10]
2500 resid1 = self.pairs[i][6][0]
2501 chain1 = self.pairs[i][6][1]
2502 resid2 = self.pairs[i][7][0]
2503 chain2 = self.pairs[i][7][1]
2505 label = str(resid1) +
":" + chain1 +
"_" + \
2506 str(resid2) +
":" + chain2
2507 output[
"CrossLinkMSSimple_Score_" + crosslinker +
2508 "_" + label] = str(ln.unprotected_evaluate(
None))
2509 output[
"CrossLinkMSSimple_Score_Linear_" + crosslinker +
2510 "_" + label] = str(pr.unprotected_evaluate(
None))
2513 output[
"CrossLinkMSSimple_Distance_" +
2519 def get_residue_index_and_chain_from_particle(p):
2526 def select_calpha_or_residue(
2531 SelectResidue=
False):
2535 residue_index=resid, atom_type=IMP.atom.AT_CA)
2537 ps = s.get_selected_particles()
2543 print(ObjectName +
" multiple residues selected for selection residue %s chain %s " % (resid, chain))
2547 residue_index=resid)
2548 ps = s.get_selected_particles()
2554 print(ObjectName +
" multiple residues selected for selection residue %s chain %s " % (resid, chain))
2557 print(ObjectName +
" residue %s chain %s does not exist" % (resid, chain))
2563 class SetupMembranePoreRestraint(object):
2569 selection_tuples_outside=
None,
2570 selection_tuples_membrane=
None,
2571 selection_tuples_inside=
None,
2572 center=(0.0,0.0,0.0),
2575 membrane_tickness=40.0,
2581 self.m = representation.prot.get_model()
2583 self.representation=representation
2586 self.z_tickness=z_tickness
2587 self.membrane_tickness=membrane_tickness
2593 def create_representation(self):
2596 h.set_name(
"Membrane_Pore"+self.label)
2597 self.representation.prot.add_child(h)
2599 inner=self.z_tickness/2-self.membrane_tickness/2
2600 outer=self.z_tickness/2+self.membrane_tickness/2
2604 particles_surface_1=[]
2605 particles_surface_2=[]
2609 v=self.math.pi/10.0*float(i)+self.math.pi/2.0
2610 u=2.0*self.math.pi/10.0*float(j)
2611 x=(self.radius+inner*self.math.cos(v))*self.math.cos(u)
2612 y=(self.radius+inner*self.math.cos(v))*self.math.sin(u)
2613 z=inner*self.math.sin(v)
2617 d.set_coordinates((x,y,z))
2621 surface_1.append((x,y,z))
2622 particles_surface_1.append(p2)
2624 x=(self.radius+outer*self.math.cos(v))*self.math.cos(u)
2625 y=(self.radius+outer*self.math.cos(v))*self.math.sin(u)
2626 z=outer*self.math.sin(v)
2630 d.set_coordinates((x,y,z))
2634 surface_2.append((x,y,z))
2635 particles_surface_2.append(p2)
2639 x=3.0*self.radius/10.0*float(i)-self.radius*1.5
2640 y=3.0*self.radius/10.0*float(j)-self.radius*1.5
2641 if(self.math.sqrt(x**2+y**2)<self.radius):
2644 for n,z
in enumerate([self.z_tickness/2+self.membrane_tickness/2,
2645 self.z_tickness/2-self.membrane_tickness/2,
2646 -self.z_tickness/2+self.membrane_tickness/2,
2647 -self.z_tickness/2-self.membrane_tickness/2]):
2651 d.set_coordinates((x,y,z))
2655 if n == 0
or n == 3:
2656 surface_1.append((x,y,z))
2657 particles_surface_1.append(p2)
2658 if n == 1
or n == 2:
2659 surface_2.append((x,y,z))
2660 particles_surface_2.append(p2)
2662 from scipy.spatial
import Delaunay
2663 tri = Delaunay(surface_1)
2670 """Add a line between the i-th and j-th points, if not in the list already"""
2671 if (i, j) in edges or (j, i) in edges:
2675 edge_points.append(points[ [i, j] ])
2677 # loop over triangles:
2678 # ia, ib, ic = indices of corner points of the triangle
2679 for ia, ib, ic in tri.vertices:
2688 for ia, ib, ic, id in tri.simplices:
2694 p1=particles_surface_1[e[0]]
2695 p2=particles_surface_1[e[1]]
2696 print(p1,p2,e[0],e[1])
2697 IMP.atom.Bonded.setup_particle(p1)
2698 IMP.atom.Bonded.setup_particle(p2)
2699 IMP.atom.create_bond(IMP.atom.Bonded(p1),IMP.atom.Bonded(p2),1)
2702 for i
in range(len(particles_surface_1)-1):
2703 p1=particles_surface_1[i]
2704 p2=particles_surface_1[i+1]
A filter which returns true if a container contains 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 bounded 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 a fixed array of same-typed values.
def deprecated_module
Mark a Python module as deprecated.
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)
Get the particles from a list of indexes.
Classes to handle different kinds of restraints.
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 spatially-proximal pairs of particles (a,b) from the two SingletonContainers A and B...
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)
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
Score a pair of particles based on the distance between their centers.
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 particles of a Model object.
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.