1 """@namespace IMP.pmi.restraints.crosslinking
2 Restraints for handling crosslinking data.
5 from __future__
import print_function
13 class _NuisancesBase(object):
14 ''' This base class is used to provide nuisance setup and interface
15 for the ISD cross-link restraints '''
17 def create_length(self):
18 ''' a nuisance on the length of the cross-link '''
19 self.lengthinit = 10.0
20 self.lengthissampled =
True
21 self.lengthminnuis = 0.0000001
22 self.lengthmaxnuis = 1000.0
25 self.lengthtrans = 0.2
26 self.length = IMP.pmi.tools.SetupNuisance(self.m, self.lengthinit,
27 self.lengthminnuis, self.lengthmaxnuis, self.lengthissampled).get_particle()
28 self.rslen.add_restraint(
36 def create_sigma(self, resolution):
37 ''' a nuisance on the structural uncertainty '''
38 self.sigmainit = resolution + 2.0
39 self.sigmaissampled =
True
40 self.sigmaminnuis = 0.0000001
41 self.sigmamaxnuis = 1000.0
45 self.sigma = IMP.pmi.tools.SetupNuisance(self.m, self.sigmainit,
46 self.sigmaminnuis, self.sigmamaxnuis, self.sigmaissampled).get_particle()
47 self.sigma_dictionary[resolution] = (
51 self.rssig.add_restraint(
60 def get_sigma(self, resolution):
61 if not resolution
in self.sigma_dictionary:
62 self.create_sigma(resolution)
63 return self.sigma_dictionary[resolution]
65 def create_psi(self, value):
66 ''' a nuisance on the inconsistency '''
68 self.psiissampled =
True
69 self.psiminnuis = 0.0000001
70 self.psimaxnuis = 0.4999999
74 self.psi = IMP.pmi.tools.SetupNuisance(self.m, self.psiinit,
75 self.psiminnuis, self.psimaxnuis,
76 self.psiissampled).get_particle()
77 self.psi_dictionary[value] = (
81 self.rspsi.add_restraint(
90 def get_psi(self, value):
91 if not value
in self.psi_dictionary:
92 self.create_psi(value)
93 return self.psi_dictionary[value]
101 this restraint allows ambiguous crosslinking between multiple copies
102 it is a variant of the SimplifiedCrossLinkMS
113 self.m = representation.prot.get_model()
119 self.outputlevel =
"low"
120 self.expdistance = expdistance
121 self.strength = strength
123 fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
127 tokens = line.split()
129 if (tokens[0] ==
"#"):
138 resolution=resolution,
140 name_is_ambiguous=
True,
142 hrc1 = [representation.hier_db.particle_to_name[p]
for p
in ps1]
144 print(
"ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
149 resolution=resolution,
151 name_is_ambiguous=
True,
153 hrc2 = [representation.hier_db.particle_to_name[p]
for p
in ps2]
155 print(
"ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
162 if self.strength
is None:
171 rad1 = rad1 / len(ps1)
172 rad2 = rad2 / len(ps2)
174 self.strength = 1 / (rad1 ** 2 + rad2 ** 2)
182 self.rs.add_restraint(cr)
183 self.pairs.append((ps1, hrc1, c1, r1, ps2, hrc2, c2, r2, cr))
197 d1.set_radius(uncertainty1)
198 d2.set_radius(uncertainty2)
202 strength = 1 / (uncertainty1 ** 2 + uncertainty2 ** 2)
209 for i
in range(npoints):
213 scores.append(cr.unprotected_evaluate(
None))
214 output.plot_xy_data(dists, scores)
216 def set_label(self, label):
218 self.rs.set_name(label)
219 for r
in self.rs.get_restraints():
222 def add_to_model(self):
223 self.m.add_restraint(self.rs)
225 def get_restraint_sets(self):
228 def get_restraint(self):
231 def set_output_level(self, level="low"):
233 self.outputlevel = level
235 def set_weight(self, weight):
237 self.rs.set_weight(weight)
239 def get_output(self):
245 score = self.weight * self.rs.unprotected_evaluate(
None)
246 output[
"_TotalScore"] = str(score)
247 output[
"ConnectivityCrossLinkMS_Score_" + self.label] = str(score)
248 for n, p
in enumerate(self.pairs):
259 for n1, p1
in enumerate(ps1):
262 for n2, p2
in enumerate(ps2):
266 label = str(r1) +
":" + name1 +
"_" + str(r2) +
":" + name2
267 output[
"ConnectivityCrossLinkMS_Distance_" +
270 label = str(r1) +
":" + c1 +
"_" + str(r2) +
":" + c2
271 output[
"ConnectivityCrossLinkMS_Score_" +
272 label] = str(self.weight * cr.unprotected_evaluate(
None))
277 class SimplifiedCrossLinkMS(object):
288 spheredistancepairscore=
True):
293 if columnmapping
is None:
295 columnmapping[
"Protein1"] = 0
296 columnmapping[
"Protein2"] = 1
297 columnmapping[
"Residue1"] = 2
298 columnmapping[
"Residue2"] = 3
300 self.m = representation.prot.get_model()
305 self.already_added_pairs = {}
307 self.outputlevel =
"low"
308 self.expdistance = expdistance
309 self.strength = strength
310 self.truncated = truncated
311 self.spheredistancepairscore = spheredistancepairscore
316 protein1 = columnmapping[
"Protein1"]
317 protein2 = columnmapping[
"Protein2"]
318 residue1 = columnmapping[
"Residue1"]
319 residue2 = columnmapping[
"Residue2"]
321 fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
325 tokens = line.split()
327 if (tokens[0] ==
"#"):
329 r1 = int(tokens[residue1])
330 c1 = tokens[protein1]
331 r2 = int(tokens[residue2])
332 c2 = tokens[protein2]
336 resolution=resolution,
338 name_is_ambiguous=
False,
342 resolution=resolution,
344 name_is_ambiguous=
False,
349 "residue %d of chain %s selects multiple particles; "
351 % (r1, c1,
"".join(p.get_name()
for p
in ps1)))
353 print(
"SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
358 "residue %d of chain %s selects multiple particles; "
360 % (r2, c2,
"".join(p.get_name()
for p
in ps2)))
362 print(
"SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
368 if (p1, p2)
in self.already_added_pairs:
369 dr = self.already_added_pairs[(p1, p2)]
370 weight = dr.get_weight()
371 dr.set_weight(weight + 1.0)
372 print(
"SimplifiedCrossLinkMS> crosslink %d %s %d %s was already found, adding 1.0 to the weight, weight is now %d" % (r1, c1, r2, c2, weight + 1.0))
378 limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
389 if self.spheredistancepairscore:
394 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2))
396 self.rs.add_restraint(dr)
397 self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
398 self.already_added_pairs[(p1, p2)] = dr
399 self.already_added_pairs[(p2, p1)] = dr
401 def set_label(self, label):
404 def add_to_model(self):
405 self.m.add_restraint(self.rs)
407 def get_restraint_sets(self):
410 def get_restraint(self):
415 for r
in self.rs.get_restraints():
416 rlist.append(IMP.core.PairRestraint.get_from(r))
419 def get_particle_pairs(self):
421 for i
in range(len(self.pairs)):
422 p0 = self.pairs[i][0]
423 p1 = self.pairs[i][1]
424 ppairs.append((p0, p1))
427 def set_output_level(self, level="low"):
429 self.outputlevel = level
431 def set_weight(self, weight):
433 self.rs.set_weight(weight)
435 def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
442 d1.set_radius(radius1)
443 d2.set_radius(radius2)
445 limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
460 for i
in range(npoints):
464 scores.append(dr.unprotected_evaluate(
None))
465 output.plot_xy_data(dists, scores)
467 def get_output(self):
473 score = self.weight * self.rs.unprotected_evaluate(
None)
474 output[
"_TotalScore"] = str(score)
475 output[
"SimplifiedCrossLinkMS_Score_" + self.label] = str(score)
476 for i
in range(len(self.pairs)):
478 p0 = self.pairs[i][0]
479 p1 = self.pairs[i][1]
480 crosslinker =
'standard'
481 ln = self.pairs[i][2]
482 resid1 = self.pairs[i][3]
483 chain1 = self.pairs[i][4]
484 resid2 = self.pairs[i][5]
485 chain2 = self.pairs[i][6]
487 label = str(resid1) +
":" + chain1 +
"_" + \
488 str(resid2) +
":" + chain2
489 output[
"SimplifiedCrossLinkMS_Score_" + crosslinker +
"_" +
490 label] = str(self.weight * ln.unprotected_evaluate(
None))
492 if self.spheredistancepairscore:
498 output[
"SimplifiedCrossLinkMS_Distance_" +
506 class SigmoidalCrossLinkMS(object):
509 self, representation, restraints_file, inflection, slope, amplitude,
510 linear_slope, resolution=
None, columnmapping=
None, csvfile=
False,
511 filters=
None, filelabel=
"None"):
518 if columnmapping
is None:
520 columnmapping[
"Protein1"] = 0
521 columnmapping[
"Protein2"] = 1
522 columnmapping[
"Residue1"] = 2
523 columnmapping[
"Residue2"] = 3
529 db = tools.get_db_from_csv(restraints_file)
531 db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
533 self.m = representation.prot.get_model()
539 self.already_added_pairs = {}
540 self.inflection = inflection
542 self.amplitude = amplitude
543 self.linear_slope = linear_slope
545 self.outputlevel =
"low"
550 protein1 = columnmapping[
"Protein1"]
551 protein2 = columnmapping[
"Protein2"]
552 residue1 = columnmapping[
"Residue1"]
553 residue2 = columnmapping[
"Residue2"]
555 indb = open(
"included." + filelabel +
".xl.db",
"w")
556 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
557 midb = open(
"missing." + filelabel +
".xl.db",
"w")
561 tokens = entry.split()
563 if (tokens[0] ==
"#"):
565 r1 = int(tokens[residue1])
566 c1 = tokens[protein1]
567 r2 = int(tokens[residue2])
568 c2 = tokens[protein2]
570 if filters
is not None:
571 if eval(tools.cross_link_db_filter_parser(filters)) ==
False:
572 exdb.write(str(entry) +
"\n")
574 indb.write(str(entry) +
"\n")
575 r1 = int(entry[residue1])
577 r2 = int(entry[residue2])
582 resolution=resolution,
584 name_is_ambiguous=
False,
588 resolution=resolution,
590 name_is_ambiguous=
False,
594 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
596 print(
"SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
597 midb.write(str(entry) +
"\n")
601 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
603 print(
"SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
604 midb.write(str(entry) +
"\n")
610 if (p1, p2)
in self.already_added_pairs:
611 dr = self.already_added_pairs[(p1, p2)]
612 weight = dr.get_weight()
613 dr.increment_amplitude(amplitude)
614 print(
"SigmoidalCrossLinkMS> crosslink %d %s %d %s was already found, adding %d to the amplitude, amplitude is now %d" % (r1, c1, r2, c2, amplitude, dr.get_amplitude()))
615 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2)
616 +
"-ampl:" + str(dr.get_amplitude()))
629 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2)
630 +
"-ampl:" + str(dr.get_amplitude()))
632 self.rs.add_restraint(dr)
634 self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
635 self.already_added_pairs[(p1, p2)] = dr
636 self.already_added_pairs[(p2, p1)] = dr
638 def set_label(self, label):
641 def add_to_model(self):
642 self.m.add_restraint(self.rs)
644 def get_restraint_sets(self):
647 def get_restraint(self):
652 for r
in self.rs.get_restraints():
653 rlist.append(IMP.core.PairRestraint.get_from(r))
656 def get_particle_pairs(self):
658 for i
in range(len(self.pairs)):
659 p0 = self.pairs[i][0]
660 p1 = self.pairs[i][1]
661 ppairs.append((p0, p1))
664 def set_output_level(self, level="low"):
666 self.outputlevel = level
668 def set_weight(self, weight):
670 self.rs.set_weight(weight)
672 def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
679 d1.set_radius(radius1)
680 d2.set_radius(radius2)
691 for i
in range(npoints):
695 scores.append(dr.unprotected_evaluate(
None))
696 output.plot_xy_data(dists, scores)
698 def get_output(self):
704 score = self.weight * self.rs.unprotected_evaluate(
None)
705 output[
"_TotalScore"] = str(score)
706 output[
"SigmoidalCrossLinkMS_Score_" + self.label] = str(score)
707 for i
in range(len(self.pairs)):
709 p0 = self.pairs[i][0]
710 p1 = self.pairs[i][1]
711 crosslinker =
'standard'
712 ln = self.pairs[i][2]
713 resid1 = self.pairs[i][3]
714 chain1 = self.pairs[i][4]
715 resid2 = self.pairs[i][5]
716 chain2 = self.pairs[i][6]
718 label = str(resid1) +
":" + chain1 +
"_" + \
719 str(resid2) +
":" + chain2
720 output[
"SigmoidalCrossLinkMS_Score_" + crosslinker +
"_" +
721 label +
"_" + self.label] = str(self.weight * ln.unprotected_evaluate(
None))
724 output[
"SigmoidalCrossLinkMS_Distance_" +
730 class ISDCrossLinkMS(_NuisancesBase):
740 def __init__(self, representation,
762 automatic_sigma_classification=
False,
763 attributes_for_label=
None):
773 if type(representation) != list:
774 representations = [representation]
776 representations = representation
778 if columnmapping
is None:
780 columnmapping[
"Protein1"] = 0
781 columnmapping[
"Protein2"] = 1
782 columnmapping[
"Residue1"] = 2
783 columnmapping[
"Residue2"] = 3
784 columnmapping[
"IDScore"] = 4
785 columnmapping[
"XLUniqueID"] = 5
791 db = IMP.pmi.tools.get_db_from_csv(restraints_file)
793 db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
795 indb = open(
"included." + filelabel +
".xl.db",
"w")
796 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
797 midb = open(
"missing." + filelabel +
".xl.db",
"w")
799 self.m = representations[0].prot.get_model()
800 self.marginal = marginal
802 if not self.marginal:
810 self.sigma_dictionary = {}
811 if not self.marginal:
812 self.psi_dictionary = {}
813 self.samplelength = samplelength
814 self.psi_is_sampled =
True
815 self.sigma_is_sampled =
True
820 self.ids_map = IMP.pmi.tools.map()
821 self.ids_map.set_map_element(20.0, 0.05)
822 self.ids_map.set_map_element(65.0, 0.01)
824 self.ids_map = ids_map
826 if radius_map
is None:
827 self.radius_map = IMP.pmi.tools.map()
828 if automatic_sigma_classification:
829 self.radius_map.set_map_element(10, 10)
830 self.radius_map.set_map_element(1, 1)
832 self.radius_map = radius_map
834 self.outputlevel =
"low"
838 self.linear.set_slope(slope)
841 protein1 = columnmapping[
"Protein1"]
842 protein2 = columnmapping[
"Protein2"]
843 residue1 = columnmapping[
"Residue1"]
844 residue2 = columnmapping[
"Residue2"]
845 idscore = columnmapping[
"IDScore"]
847 xluniqueid = columnmapping[
"XLUniqueID"]
857 uniqueid_restraints_map = {}
859 for nxl, entry
in enumerate(db):
861 if not jackknifing
is None:
868 raise NotImplementedError(
"jackknifing not yet implemented")
871 tokens = entry.split()
876 if (tokens[0] ==
"#"):
879 r1 = int(tokens[residue1])
880 c1 = tokens[protein1]
881 r2 = int(tokens[residue2])
882 c2 = tokens[protein2]
884 if offset_dict
is not None:
885 if c1
in offset_dict: r1+=offset_dict[c1]
886 if c2
in offset_dict: r2+=offset_dict[c2]
888 if rename_dict
is not None:
889 if c1
in rename_dict: c1=rename_dict[c1]
890 if c2
in rename_dict: c2=rename_dict[c2]
895 ids = float(tokens[idscore])
896 if xluniqueid
is None:
899 xlid = tokens[xluniqueid]
901 print(
"this line was not accessible " + str(entry))
902 if residue1
not in entry: print(str(residue1)+
" keyword not in database")
903 if residue2
not in entry: print(str(residue2)+
" keyword not in database")
904 if protein1
not in entry: print(str(protein1)+
" keyword not in database")
905 if protein2
not in entry: print(str(protein2)+
" keyword not in database")
906 if idscore
not in entry: print(str(idscore)+
" keyword not in database")
907 if xluniqueid
not in entry: print(str(xluniqueid)+
" keyword not in database")
911 if filters
is not None:
913 exdb.write(str(entry) +
"\n")
917 r1 = int(entry[residue1])
919 r2 = int(entry[residue2])
922 if offset_dict
is not None:
923 if c1
in offset_dict: r1+=offset_dict[c1]
924 if c2
in offset_dict: r2+=offset_dict[c2]
926 if rename_dict
is not None:
927 if c1
in rename_dict: c1=rename_dict[c1]
928 if c2
in rename_dict: c2=rename_dict[c2]
934 ids = float(entry[idscore])
937 if xluniqueid
is None:
940 xlid = entry[xluniqueid]
943 print(
"this line was not accessible " + str(entry))
944 if residue1
not in entry: print(str(residue1)+
" keyword not in database")
945 if residue2
not in entry: print(str(residue2)+
" keyword not in database")
946 if protein1
not in entry: print(str(protein1)+
" keyword not in database")
947 if protein2
not in entry: print(str(protein2)+
" keyword not in database")
948 if idscore
not in entry: print(str(idscore)+
" keyword not in database")
949 if xluniqueid
not in entry: print(str(xluniqueid)+
" keyword not in database")
952 for nstate, r
in enumerate(representations):
956 resolution=resolution,
958 name_is_ambiguous=
False,
962 resolution=resolution,
964 name_is_ambiguous=
False,
968 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
970 print(
"ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
971 midb.write(str(entry) +
"\n")
975 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
977 print(
"ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
978 midb.write(str(entry) +
"\n")
988 if xlid
in uniqueid_restraints_map:
989 print(
"getting a crosslink restraint from id %s" % str(xlid))
990 dr = uniqueid_restraints_map[xlid]
993 if not self.marginal:
994 if not self.samplelength:
995 print(
"generating a new crosslink restraint")
1003 self.create_length()
1009 if not self.samplelength:
1010 dr = IMP.isd_emxl.CrossLinkMSMarginalRestraint(
1016 self.create_length()
1017 dr = IMP.isd_emxl.CrossLinkMSMarginalRestraint(
1021 uniqueid_restraints_map[xlid] = dr
1023 mappedr1 = self.radius_map.get_map_element(
1025 sigma1 = self.get_sigma(mappedr1)[0]
1026 mappedr2 = self.radius_map.get_map_element(
1028 sigma2 = self.get_sigma(mappedr2)[0]
1030 if not self.marginal:
1031 psival = self.ids_map.get_map_element(ids)
1032 psi = self.get_psi(psival)[0]
1035 p1i = p1.get_particle_index()
1036 p2i = p2.get_particle_index()
1042 if not self.marginal:
1044 dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
1045 print(
"--------------")
1046 print(
"ISDCrossLinkMS: generating cross-link restraint between")
1047 print(
"ISDCrossLinkMS: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
1048 print(
"ISDCrossLinkMS: with sigma1 %f sigma2 %f psi %s" % (mappedr1, mappedr2, psival))
1049 print(
"ISDCrossLinkMS: between particles %s and %s" % (p1.get_name(), p2.get_name()))
1050 print(
"==========================================\n")
1053 dr.add_contribution((p1i, p2i), (s1i, s2i))
1054 print(
"--------------")
1055 print(
"ISDCrossLinkMS: generating marginal cross-link restraint between")
1056 print(
"ISDCrossLinkMS: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
1057 print(
"ISDCrossLinkMS: with sigma1 %f sigma2 %f" % (mappedr1, mappedr2))
1058 print(
"ISDCrossLinkMS: between particles %s and %s" % (p1.get_name(), p2.get_name()))
1059 print(
"==========================================\n")
1061 indb.write(str(entry) +
"\n")
1063 restraints.append(dr)
1070 xlattribute =
"intrarb"
1072 xlattribute =
"interrb"
1075 if not attributes_for_label
is None:
1076 for a
in attributes_for_label:
1077 xlattribute = xlattribute +
"_" + str(entry[a])
1079 xlattribute = xlattribute +
"-State:" + str(nstate)
1082 xlattribute +
"-" + c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2) +
"_" + self.label)
1086 xlattribute +
"-" + c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2) +
"_" + self.label)
1087 self.rslin.add_restraint(pr)
1105 self.rs.add_restraint(lw)
1108 def set_slope_linear_term(self, slope):
1109 self.linear.set_slope(slope)
1111 def set_label(self, label):
1114 def add_to_model(self):
1115 self.m.add_restraint(self.rs)
1116 if not self.marginal:
1117 self.m.add_restraint(self.rspsi)
1118 self.m.add_restraint(self.rssig)
1119 self.m.add_restraint(self.rslen)
1120 self.m.add_restraint(self.rslin)
1122 def get_hierarchies(self):
1125 def get_restraint_sets(self):
1128 def get_restraint(self):
1131 def get_restraint_for_rmf(self):
1136 for r
in self.rs.get_restraints():
1137 rlist.append(IMP.core.PairRestraint.get_from(r))
1140 def get_particle_pairs(self):
1142 for i
in range(len(self.pairs)):
1143 p0 = self.pairs[i][0]
1144 p1 = self.pairs[i][1]
1145 ppairs.append((p0, p1))
1148 def set_output_level(self, level="low"):
1150 self.outputlevel = level
1152 def set_psi_is_sampled(self, is_sampled=True):
1153 self.psi_is_sampled = is_sampled
1155 def set_sigma_is_sampled(self, is_sampled=True):
1156 self.sigma_is_sampled = is_sampled
1158 def get_output(self):
1164 score = self.rs.unprotected_evaluate(
None)
1165 output[
"_TotalScore"] = str(score)
1166 output[
"ISDCrossLinkMS_Data_Score_" + self.label] = str(score)
1167 output[
"ISDCrossLinkMS_PriorSig_Score_" +
1168 self.label] = self.rssig.unprotected_evaluate(
None)
1169 if not self.marginal:
1170 output[
"ISDCrossLinkMS_PriorPsi_Score_" +
1171 self.label] = self.rspsi.unprotected_evaluate(
None)
1172 output[
"ISDCrossLinkMS_Linear_Score_" +
1173 self.label] = self.rslin.unprotected_evaluate(
None)
1174 for i
in range(len(self.pairs)):
1176 p0 = self.pairs[i][0]
1177 p1 = self.pairs[i][1]
1178 ln = self.pairs[i][2]
1179 resid1 = self.pairs[i][3]
1180 chain1 = self.pairs[i][4]
1181 resid2 = self.pairs[i][5]
1182 chain2 = self.pairs[i][6]
1183 attribute = self.pairs[i][7]
1184 rad1 = self.pairs[i][8]
1185 rad2 = self.pairs[i][9]
1186 psi = self.pairs[i][10]
1187 xlid= self.pairs[i][11]
1189 label = attribute +
"-" + \
1190 str(resid1) +
":" + chain1 +
"_" + str(resid2) +
":" + \
1191 chain2 +
"-" + str(rad1) +
"-" + str(rad2) +
"-" + str(psi)
1192 output[
"ISDCrossLinkMS_Score_" +
1193 label +
"_" + self.label] = str(-self.log(ln.unprotected_evaluate(
None)))
1196 output[
"ISDCrossLinkMS_Distance_" +
1199 if not self.marginal:
1200 for psiindex
in self.psi_dictionary:
1201 output[
"ISDCrossLinkMS_Psi_" +
1202 str(psiindex) +
"_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
1204 for resolution
in self.sigma_dictionary:
1205 output[
"ISDCrossLinkMS_Sigma_" +
1206 str(resolution) +
"_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
1208 if self.samplelength:
1209 output[
"ISDCrossLinkMS_Length_" +
1210 str(resolution) +
"_" + self.label] = str(self.length.get_scale())
1214 def get_particles_to_sample(self):
1217 for resolution
in self.sigma_dictionary:
1218 if self.sigma_dictionary[resolution][2]
and self.sigma_is_sampled:
1219 ps[
"Nuisances_ISDCrossLinkMS_Sigma_" + str(resolution) +
"_" + self.label] =\
1220 ([self.sigma_dictionary[resolution][0]],
1221 self.sigma_dictionary[resolution][1])
1223 if (
not self.marginal)
and self.psi_is_sampled:
1224 for psiindex
in self.psi_dictionary:
1225 if self.psi_dictionary[psiindex][2]:
1226 ps[
"Nuisances_ISDCrossLinkMS_Psi_" +
1227 str(psiindex) +
"_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
1229 if self.samplelength:
1230 if self.lengthissampled:
1231 ps[
"Nuisances_ISDCrossLinkMS_Length_" +
1232 self.label] = ([self.length], self.lengthtrans)
1237 class CysteineCrossLinkRestraint(object):
1242 def __init__(self, representations, filename, cbeta=False,
1243 betatuple=(0.03, 0.1),
1244 disttuple=(0.0, 25.0, 1000),
1245 omegatuple=(1.0, 1000.0, 50),
1246 sigmatuple=(0.3, 0.3, 1),
1248 sigmaissampled=
False,
1249 weightissampled=
True,
1250 epsilonissampled=
True
1257 self.representations = representations
1258 self.m = self.representations[0].prot.get_model()
1261 self.epsilonmaxtrans = 0.01
1262 self.sigmamaxtrans = 0.1
1263 self.betamaxtrans = 0.01
1264 self.weightmaxtrans = 0.1
1266 self.outputlevel =
"low"
1267 self.betaissampled = betaissampled
1268 self.sigmaissampled = sigmaissampled
1269 self.weightissampled = weightissampled
1270 self.epsilonissampled = epsilonissampled
1272 betalower = betatuple[0]
1273 betaupper = betatuple[1]
1280 self.beta = IMP.pmi.tools.SetupNuisance(
1285 betaissampled).get_particle(
1288 self.sigma = IMP.pmi.tools.SetupNuisance(
1293 sigmaissampled).get_particle()
1295 self.weight = IMP.pmi.tools.SetupWeight(
1297 weightissampled).get_particle(
1301 fl = IMP.pmi.tools.open_file_or_inline_text(filename)
1310 if t[5]
in self.epsilons:
1311 if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
1312 self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
1314 self.epsilons[t[5]] = IMP.pmi.tools.SetupNuisance(self.m,
1315 0.01, 0.01, 1.0 - float(t[4]), epsilonissampled).get_particle()
1316 up = self.epsilons[t[5]].get_upper()
1317 low = self.epsilons[t[5]].get_lower()
1319 self.epsilons[t[5]].set_upper(low)
1321 data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
1325 crossdata = IMP.pmi.tools.get_cross_link_data(
1326 "cysteine",
"cysteine_CA_FES.txt.standard",
1327 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
1329 crossdata = IMP.pmi.tools.get_cross_link_data(
1330 "cysteine",
"cysteine_CB_FES.txt.standard",
1331 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
1334 fmod_grid = IMP.pmi.tools.get_grid(0.0, 1.0, 300,
True)
1335 omega2_grid = IMP.pmi.tools.get_log_grid(0.001, 10000.0, 100)
1336 beta_grid = IMP.pmi.tools.get_log_grid(betalower, betaupper, betangrid)
1339 print(
"--------------")
1340 print(
"CysteineCrossLink: attempting to create a restraint " + str(d))
1359 self.epsilons[epslabel],
1365 for i, representation
in enumerate(self.representations):
1372 resolution=1, name=chain1,
1373 name_is_ambiguous=
False, residue=resid1)[0]
1379 resolution=1, name=chain2,
1380 name_is_ambiguous=
False, residue=resid2)[0]
1389 for t
in range(-1, 2):
1391 resolution=1, name=chain1,
1392 name_is_ambiguous=
False, residue=resid1 + t)
1398 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
1401 resolution=1, name=chain2,
1402 name_is_ambiguous=
False, residue=resid2 + t)
1408 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
1411 if (p1
is not None and p2
is not None):
1412 ccl.add_contribution(p1, p2)
1416 print(
"Distance_" + str(resid1) +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2,
IMP.core.get_distance(d1, d2))
1419 if (len(p1) == 3
and len(p2) == 3):
1420 p11n = p1[0].get_name()
1421 p12n = p1[1].get_name()
1422 p13n = p1[2].get_name()
1423 p21n = p2[0].get_name()
1424 p22n = p2[1].get_name()
1425 p23n = p2[2].get_name()
1427 print(
"CysteineCrossLink: generating CB cysteine cross-link restraint between")
1428 print(
"CysteineCrossLink: residue %d of chain %s and residue %d of chain %s" % (resid1, chain1, resid2, chain2))
1429 print(
"CysteineCrossLink: between particles %s %s %s and %s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
1431 ccl.add_contribution(p1, p2)
1434 self.rs.add_restraint(ccl)
1435 ccl.set_name(
"CysteineCrossLink_" + str(resid1)
1436 +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2)
1438 def set_label(self, label):
1441 def add_to_model(self):
1442 self.m.add_restraint(self.rs)
1444 def get_particles_to_sample(self):
1446 if self.epsilonissampled:
1447 for eps
in self.epsilons.keys():
1448 ps[
"Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
1449 self.label] = ([self.epsilons[eps]], self.epsilonmaxtrans)
1450 if self.betaissampled:
1451 ps[
"Nuisances_CysteineCrossLinkRestraint_beta_" +
1452 self.label] = ([self.beta], self.betamaxtrans)
1453 if self.weightissampled:
1454 ps[
"Weights_CysteineCrossLinkRestraint_" +
1455 self.label] = ([self.weight], self.weightmaxtrans)
1456 if self.sigmaissampled:
1457 ps[
"Nuisances_CysteineCrossLinkRestraint_" +
1458 self.label] = ([self.sigma], self.sigmamaxtrans)
1461 def set_output_level(self, level="low"):
1463 self.outputlevel = level
1465 def get_restraint(self):
1468 def get_sigma(self):
1471 def get_output(self):
1474 score = self.rs.unprotected_evaluate(
None)
1475 output[
"_TotalScore"] = str(score)
1476 output[
"CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
1477 output[
"CysteineCrossLinkRestraint_sigma_" +
1478 self.label] = str(self.sigma.get_scale())
1479 for eps
in self.epsilons.keys():
1480 output[
"CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
1481 self.label] = str(self.epsilons[eps].get_scale())
1482 output[
"CysteineCrossLinkRestraint_beta_" +
1483 self.label] = str(self.beta.get_scale())
1484 for n
in range(self.weight.get_number_of_states()):
1485 output[
"CysteineCrossLinkRestraint_weights_" +
1486 str(n) +
"_" + self.label] = str(self.weight.get_weight(n))
1488 if self.outputlevel ==
"high":
1489 for rst
in self.rs.get_restraints():
1490 output[
"CysteineCrossLinkRestraint_Total_Frequency_" +
1491 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
1492 "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_model_frequency()
1493 output[
"CysteineCrossLinkRestraint_Standard_Error_" +
1494 IMP.isd.CysteineCrossLinkRestraint.get_from(
1497 + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_standard_error()
1498 if len(self.representations) > 1:
1499 for i
in range(len(self.prots)):
1500 output[
"CysteineCrossLinkRestraint_Frequency_Contribution_" +
1501 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
1502 "_State_" + str(i) +
"_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_frequencies()[i]
A function that is harmonic over an interval.
A restraint for ambiguous cross-linking MS data and multiple state approach.
Ints get_index(const kernel::ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
Object used to hold a set of restraints.
Low level functionality (logging, error handling, profiling, command line flags etc) that is used by ...
kernel::RestraintsTemp get_restraints(const Subset &s, const ParticleStatesTable *pst, const DependencyGraph &dg, kernel::RestraintSet *rs)
Add uncertainty to a particle.
A score on the distance between the surfaces of two spheres.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
A class to store an fixed array of same-typed values.
static XYZR setup_particle(kernel::Model *m, ParticleIndex pi)
Simple sigmoidal score calculated between sphere surfaces.
kernel::Restraint * create_connectivity_restraint(const Selections &s, double x0, double k, std::string name="Connectivity%1%")
Create a restraint connecting the selections.
A decorator for a particle with x,y,z coordinates.
Class to handle individual model particles.
this restraint allows ambiguous crosslinking between multiple copies it is a variant of the Simplifie...
Classes for writing output files and processing them.
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
static bool get_is_setup(const IMP::kernel::ParticleAdaptor &p)
Calculate the -Log of a list of restraints.
Applies a PairScore to a Pair.
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A restraint for cysteine cross-linking data.
A decorator for a particle with x,y,z coordinates and a radius.