IMP  2.4.0
The Integrative Modeling Platform
crosslinking.py
1 """@namespace IMP.pmi.restraints.crosslinking
2 Restraints for handling crosslinking data.
3 """
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.core
8 import IMP.base
9 import IMP.algebra
10 import IMP.atom
11 import IMP.container
12 
13 class _NuisancesBase(object):
14  ''' This base class is used to provide nuisance setup and interface
15  for the ISD cross-link restraints '''
16 
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
23  self.lengthmin = 6.0
24  self.lengthmax = 30.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(
30  self.m,
31  self.length,
32  1000000000.0,
33  self.lengthmax,
34  self.lengthmin))
35 
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
42  self.sigmamin = 0.01
43  self.sigmamax = 100.0
44  self.sigmatrans = 0.5
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] = (
48  self.sigma,
49  self.sigmatrans,
50  self.sigmaissampled)
51  self.rssig.add_restraint(
53  self.m,
54  self.sigma,
55  1000000000.0,
56  self.sigmamax,
57  self.sigmamin))
58  # self.rssig.add_restraint(IMP.isd.JeffreysRestraint(self.sigma))
59 
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]
64 
65  def create_psi(self, value):
66  ''' a nuisance on the inconsistency '''
67  self.psiinit = value
68  self.psiissampled = True
69  self.psiminnuis = 0.0000001
70  self.psimaxnuis = 0.4999999
71  self.psimin = 0.01
72  self.psimax = 0.49
73  self.psitrans = 0.1
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] = (
78  self.psi,
79  self.psitrans,
80  self.psiissampled)
81  self.rspsi.add_restraint(
83  self.m,
84  self.psi,
85  1000000000.0,
86  self.psimax,
87  self.psimin))
88  self.rspsi.add_restraint(IMP.isd.JeffreysRestraint(self.m, self.psi))
89 
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]
94 
95 
96 
97 
99 
100  '''
101  this restraint allows ambiguous crosslinking between multiple copies
102  it is a variant of the SimplifiedCrossLinkMS
103  '''
104 
105  def __init__(
106  self,
107  representation,
108  restraints_file,
109  expdistance,
110  strength=None,
111  resolution=None):
112 
113  self.m = representation.prot.get_model()
114  self.rs = IMP.RestraintSet(self.m, 'data')
115  self.weight = 1.0
116  self.label = "None"
117  self.pairs = []
118 
119  self.outputlevel = "low"
120  self.expdistance = expdistance
121  self.strength = strength
122 
123  fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
124 
125  for line in fl:
126 
127  tokens = line.split()
128  # skip character
129  if (tokens[0] == "#"):
130  continue
131  r1 = int(tokens[2])
132  c1 = tokens[0]
133  r2 = int(tokens[3])
134  c2 = tokens[1]
135 
136  ps1 = IMP.pmi.tools.select(
137  representation,
138  resolution=resolution,
139  name=c1,
140  name_is_ambiguous=True,
141  residue=r1)
142  hrc1 = [representation.hier_db.particle_to_name[p] for p in ps1]
143  if len(ps1) == 0:
144  print("ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
145  continue
146 
147  ps2 = IMP.pmi.tools.select(
148  representation,
149  resolution=resolution,
150  name=c2,
151  name_is_ambiguous=True,
152  residue=r2)
153  hrc2 = [representation.hier_db.particle_to_name[p] for p in ps2]
154  if len(ps2) == 0:
155  print("ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
156  continue
157 
158  s1 = IMP.atom.Selection(ps1)
159  s2 = IMP.atom.Selection(ps2)
160 
161  # calculate the radii to estimate the slope of the restraint
162  if self.strength is None:
163  rad1 = 0
164  rad2 = 0
165  for p in ps1:
166  rad1 += IMP.pmi.Uncertainty(p).get_uncertainty()
167 
168  for p in ps2:
169  rad2 += IMP.pmi.Uncertainty(p).get_uncertainty()
170 
171  rad1 = rad1 / len(ps1)
172  rad2 = rad2 / len(ps2)
173 
174  self.strength = 1 / (rad1 ** 2 + rad2 ** 2)
175 
176  sels = [s1, s2]
178  sels,
179  self.expdistance,
180  self.strength)
181 
182  self.rs.add_restraint(cr)
183  self.pairs.append((ps1, hrc1, c1, r1, ps2, hrc2, c2, r2, cr))
184 
185  def plot_restraint(
186  self,
187  uncertainty1,
188  uncertainty2,
189  maxdist=50,
190  npoints=10):
191  import IMP.pmi.output as output
192 
193  p1 = IMP.Particle(self.m)
194  p2 = IMP.Particle(self.m)
197  d1.set_radius(uncertainty1)
198  d2.set_radius(uncertainty2)
199  s1 = IMP.atom.Selection(p1)
200  s2 = IMP.atom.Selection(p2)
201  sels = [s1, s2]
202  strength = 1 / (uncertainty1 ** 2 + uncertainty2 ** 2)
204  sels,
205  self.expdistance,
206  strength)
207  dists = []
208  scores = []
209  for i in range(npoints):
210  d2.set_coordinates(
211  IMP.algebra.Vector3D(maxdist / npoints * float(i), 0, 0))
212  dists.append(IMP.core.get_distance(d1, d2))
213  scores.append(cr.unprotected_evaluate(None))
214  output.plot_xy_data(dists, scores)
215 
216  def set_label(self, label):
217  self.label = label
218  self.rs.set_name(label)
219  for r in self.rs.get_restraints():
220  r.set_name(label)
221 
222  def add_to_model(self):
223  self.m.add_restraint(self.rs)
224 
225  def get_restraint_sets(self):
226  return self.rs
227 
228  def get_restraint(self):
229  return self.rs
230 
231  def set_output_level(self, level="low"):
232  # this might be "low" or "high"
233  self.outputlevel = level
234 
235  def set_weight(self, weight):
236  self.weight = weight
237  self.rs.set_weight(weight)
238 
239  def get_output(self):
240  # content of the crosslink database pairs
241  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
242  self.m.update()
243 
244  output = {}
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):
249 
250  ps1 = p[0]
251  hrc1 = p[1]
252  c1 = p[2]
253  r1 = p[3]
254  ps2 = p[4]
255  hrc2 = p[5]
256  c2 = p[6]
257  r2 = p[7]
258  cr = p[8]
259  for n1, p1 in enumerate(ps1):
260  name1 = hrc1[n1]
261 
262  for n2, p2 in enumerate(ps2):
263  name2 = hrc2[n2]
264  d1 = IMP.core.XYZR(p1)
265  d2 = IMP.core.XYZR(p2)
266  label = str(r1) + ":" + name1 + "_" + str(r2) + ":" + name2
267  output["ConnectivityCrossLinkMS_Distance_" +
268  label] = str(IMP.core.get_distance(d1, d2))
269 
270  label = str(r1) + ":" + c1 + "_" + str(r2) + ":" + c2
271  output["ConnectivityCrossLinkMS_Score_" +
272  label] = str(self.weight * cr.unprotected_evaluate(None))
273 
274  return output
275 
276 
277 class SimplifiedCrossLinkMS(object):
278 
279  def __init__(
280  self,
281  representation,
282  restraints_file,
283  expdistance,
284  strength,
285  resolution=None,
286  columnmapping=None,
287  truncated=True,
288  spheredistancepairscore=True):
289  # columnindexes is a list of column indexes for protein1, protein2, residue1, residue2
290  # by default column 0 = protein1; column 1 = protein2; column 2 =
291  # residue1; column 3 = residue2
292 
293  if columnmapping is None:
294  columnmapping = {}
295  columnmapping["Protein1"] = 0
296  columnmapping["Protein2"] = 1
297  columnmapping["Residue1"] = 2
298  columnmapping["Residue2"] = 3
299 
300  self.m = representation.prot.get_model()
301  self.rs = IMP.RestraintSet(self.m, 'data')
302  self.weight = 1.0
303  self.label = "None"
304  self.pairs = []
305  self.already_added_pairs = {}
306 
307  self.outputlevel = "low"
308  self.expdistance = expdistance
309  self.strength = strength
310  self.truncated = truncated
311  self.spheredistancepairscore = spheredistancepairscore
312 
313  # fill the cross-linker pmfs
314  # to accelerate the init the list listofxlinkertypes might contain only
315  # yht needed crosslinks
316  protein1 = columnmapping["Protein1"]
317  protein2 = columnmapping["Protein2"]
318  residue1 = columnmapping["Residue1"]
319  residue2 = columnmapping["Residue2"]
320 
321  fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
322 
323  for line in fl:
324 
325  tokens = line.split()
326  # skip character
327  if (tokens[0] == "#"):
328  continue
329  r1 = int(tokens[residue1])
330  c1 = tokens[protein1]
331  r2 = int(tokens[residue2])
332  c2 = tokens[protein2]
333 
334  ps1 = IMP.pmi.tools.select(
335  representation,
336  resolution=resolution,
337  name=c1,
338  name_is_ambiguous=False,
339  residue=r1)
340  ps2 = IMP.pmi.tools.select(
341  representation,
342  resolution=resolution,
343  name=c2,
344  name_is_ambiguous=False,
345  residue=r2)
346 
347  if len(ps1) > 1:
348  raise ValueError(
349  "residue %d of chain %s selects multiple particles; "
350  "particles are: %s"
351  % (r1, c1, "".join(p.get_name() for p in ps1)))
352  elif len(ps1) == 0:
353  print("SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
354  continue
355 
356  if len(ps2) > 1:
357  raise ValueError(
358  "residue %d of chain %s selects multiple particles; "
359  "particles are: %s"
360  % (r2, c2, "".join(p.get_name() for p in ps2)))
361  elif len(ps2) == 0:
362  print("SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
363  continue
364 
365  p1 = ps1[0]
366  p2 = ps2[0]
367 
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))
373  continue
374 
375  else:
376 
377  if self.truncated:
378  limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
380  self.expdistance,
381  self.strength,
382  self.expdistance +
383  15.,
384  limit)
385  else:
387  self.expdistance,
388  self.strength)
389  if self.spheredistancepairscore:
391  else:
393  dr = IMP.core.PairRestraint(df, (p1, p2))
394  dr.set_name(c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2))
395 
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
400 
401  def set_label(self, label):
402  self.label = label
403 
404  def add_to_model(self):
405  self.m.add_restraint(self.rs)
406 
407  def get_restraint_sets(self):
408  return self.rs
409 
410  def get_restraint(self):
411  return self.rs
412 
413  def get_restraints(self):
414  rlist = []
415  for r in self.rs.get_restraints():
416  rlist.append(IMP.core.PairRestraint.get_from(r))
417  return rlist
418 
419  def get_particle_pairs(self):
420  ppairs = []
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))
425  return ppairs
426 
427  def set_output_level(self, level="low"):
428  # this might be "low" or "high"
429  self.outputlevel = level
430 
431  def set_weight(self, weight):
432  self.weight = weight
433  self.rs.set_weight(weight)
434 
435  def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
436  import IMP.pmi.output as output
437 
438  p1 = IMP.Particle(self.m)
439  p2 = IMP.Particle(self.m)
442  d1.set_radius(radius1)
443  d2.set_radius(radius2)
444  if self.truncated:
445  limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
447  self.expdistance,
448  self.strength,
449  self.expdistance +
450  15.,
451  limit)
452  else:
454  self.expdistance,
455  self.strength)
457  dr = IMP.core.PairRestraint(df, (p1, p2))
458  dists = []
459  scores = []
460  for i in range(npoints):
461  d2.set_coordinates(
462  IMP.algebra.Vector3D(maxdist / npoints * float(i), 0, 0))
463  dists.append(IMP.core.get_distance(d1, d2))
464  scores.append(dr.unprotected_evaluate(None))
465  output.plot_xy_data(dists, scores)
466 
467  def get_output(self):
468  # content of the crosslink database pairs
469  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
470  self.m.update()
471 
472  output = {}
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)):
477 
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]
486 
487  label = str(resid1) + ":" + chain1 + "_" + \
488  str(resid2) + ":" + chain2
489  output["SimplifiedCrossLinkMS_Score_" + crosslinker + "_" +
490  label] = str(self.weight * ln.unprotected_evaluate(None))
491 
492  if self.spheredistancepairscore:
493  d0 = IMP.core.XYZR(p0)
494  d1 = IMP.core.XYZR(p1)
495  else:
496  d0 = IMP.core.XYZ(p0)
497  d1 = IMP.core.XYZ(p1)
498  output["SimplifiedCrossLinkMS_Distance_" +
499  label] = str(IMP.core.get_distance(d0, d1))
500 
501  return output
502 
503 #
504 
505 
506 class SigmoidalCrossLinkMS(object):
507 
508  def __init__(
509  self, representation, restraints_file, inflection, slope, amplitude,
510  linear_slope, resolution=None, columnmapping=None, csvfile=False,
511  filters=None, filelabel="None"):
512  # columnindexes is a list of column indexes for protein1, protein2, residue1, residue2
513  # by default column 0 = protein1; column 1 = protein2; column 2 = residue1; column 3 = residue2
514  # the filters applies to the csvfile, the format is
515  # filters=[("Field1",">|<|=|>=|<=",value),("Field2","is","String"),("Field2","in","String")]
516  import IMP.pmi.tools as tools
517 
518  if columnmapping is None:
519  columnmapping = {}
520  columnmapping["Protein1"] = 0
521  columnmapping["Protein2"] = 1
522  columnmapping["Residue1"] = 2
523  columnmapping["Residue2"] = 3
524 
525  if csvfile:
526  # in case the file is a cvs file
527  # columnmapping will contain the field names
528  # that compare in the first line of the cvs file
529  db = tools.get_db_from_csv(restraints_file)
530  else:
531  db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
532 
533  self.m = representation.prot.get_model()
534  self.rs = IMP.RestraintSet(self.m, 'data')
535  self.weight = 1.0
536 
537  self.label = "None"
538  self.pairs = []
539  self.already_added_pairs = {}
540  self.inflection = inflection
541  self.slope = slope
542  self.amplitude = amplitude
543  self.linear_slope = linear_slope
544 
545  self.outputlevel = "low"
546 
547  # fill the cross-linker pmfs
548  # to accelerate the init the list listofxlinkertypes might contain only
549  # yht needed crosslinks
550  protein1 = columnmapping["Protein1"]
551  protein2 = columnmapping["Protein2"]
552  residue1 = columnmapping["Residue1"]
553  residue2 = columnmapping["Residue2"]
554 
555  indb = open("included." + filelabel + ".xl.db", "w")
556  exdb = open("excluded." + filelabel + ".xl.db", "w")
557  midb = open("missing." + filelabel + ".xl.db", "w")
558 
559  for entry in db:
560  if not csvfile:
561  tokens = entry.split()
562  # skip character
563  if (tokens[0] == "#"):
564  continue
565  r1 = int(tokens[residue1])
566  c1 = tokens[protein1]
567  r2 = int(tokens[residue2])
568  c2 = tokens[protein2]
569  else:
570  if filters is not None:
571  if eval(tools.cross_link_db_filter_parser(filters)) == False:
572  exdb.write(str(entry) + "\n")
573  continue
574  indb.write(str(entry) + "\n")
575  r1 = int(entry[residue1])
576  c1 = entry[protein1]
577  r2 = int(entry[residue2])
578  c2 = entry[protein2]
579 
580  ps1 = IMP.pmi.tools.select(
581  representation,
582  resolution=resolution,
583  name=c1,
584  name_is_ambiguous=False,
585  residue=r1)
586  ps2 = IMP.pmi.tools.select(
587  representation,
588  resolution=resolution,
589  name=c2,
590  name_is_ambiguous=False,
591  residue=r2)
592 
593  if len(ps1) > 1:
594  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
595  elif len(ps1) == 0:
596  print("SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
597  midb.write(str(entry) + "\n")
598  continue
599 
600  if len(ps2) > 1:
601  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
602  elif len(ps2) == 0:
603  print("SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
604  midb.write(str(entry) + "\n")
605  continue
606 
607  p1 = ps1[0]
608  p2 = ps2[0]
609 
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()))
617  continue
618 
619  else:
620 
622  self.m,
623  p1,
624  p2,
625  self.inflection,
626  self.slope,
627  self.amplitude,
628  self.linear_slope)
629  dr.set_name(c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2)
630  + "-ampl:" + str(dr.get_amplitude()))
631 
632  self.rs.add_restraint(dr)
633 
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
637 
638  def set_label(self, label):
639  self.label = label
640 
641  def add_to_model(self):
642  self.m.add_restraint(self.rs)
643 
644  def get_restraint_sets(self):
645  return self.rs
646 
647  def get_restraint(self):
648  return self.rs
649 
650  def get_restraints(self):
651  rlist = []
652  for r in self.rs.get_restraints():
653  rlist.append(IMP.core.PairRestraint.get_from(r))
654  return rlist
655 
656  def get_particle_pairs(self):
657  ppairs = []
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))
662  return ppairs
663 
664  def set_output_level(self, level="low"):
665  # this might be "low" or "high"
666  self.outputlevel = level
667 
668  def set_weight(self, weight):
669  self.weight = weight
670  self.rs.set_weight(weight)
671 
672  def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
673  import IMP.pmi.output as output
674 
675  p1 = IMP.Particle(self.m)
676  p2 = IMP.Particle(self.m)
679  d1.set_radius(radius1)
680  d2.set_radius(radius2)
682  self.m,
683  p1,
684  p2,
685  self.inflection,
686  self.slope,
687  self.amplitude,
688  self.linear_slope)
689  dists = []
690  scores = []
691  for i in range(npoints):
692  d2.set_coordinates(
693  IMP.algebra.Vector3D(maxdist / npoints * float(i), 0, 0))
694  dists.append(IMP.core.get_distance(d1, d2))
695  scores.append(dr.unprotected_evaluate(None))
696  output.plot_xy_data(dists, scores)
697 
698  def get_output(self):
699  # content of the crosslink database pairs
700  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
701  self.m.update()
702 
703  output = {}
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)):
708 
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]
717 
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))
722  d0 = IMP.core.XYZR(p0)
723  d1 = IMP.core.XYZR(p1)
724  output["SigmoidalCrossLinkMS_Distance_" +
725  label + "_" + self.label] = str(IMP.core.get_distance(d0, d1))
726 
727  return output
728 
729 
730 class ISDCrossLinkMS(_NuisancesBase):
731  import IMP.isd
732  try:
733  import IMP.isd_emxl
734  no_isd_emxl = False
735  except:
736  no_isd_emxl = True
737  import IMP.pmi.tools
738  from math import log
739 
740  def __init__(self, representation,
741  restraints_file,
742  length,
743  jackknifing=None,
744  # jackknifing (Float [0,1]) default=None; percentage of data to be
745  # removed randomly
746  resolution=None,
747  # resolution (Non-negative Integer) default=None; percentage of data
748  # to be removed randomly
749  slope=0.0,
750  inner_slope=0.0,
751  columnmapping=None,
752  rename_dict=None,
753  offset_dict=None,
754  csvfile=False,
755  samplelength=False,
756  ids_map=None,
757  radius_map=None,
758  filters=None,
759  label="None",
760  filelabel="None",
761  marginal=False,
762  automatic_sigma_classification=False,
763  attributes_for_label=None):
764 
765  # columnindexes is a list of column indexes for protein1, protein2, residue1, residue2,idscore, XL unique id
766  # by default column 0 = protein1; column 1 = protein2; column 2 = residue1; column 3 = residue2;
767  # column 4 = idscores
768  # attributes_for_label: anithing in the csv database that must be added to the label
769  # slope is the slope defined on the linear function
770  # inner_slope is the slope defined on the restraint dierectly
771  # suggestion: do not use both!
772 
773  if type(representation) != list:
774  representations = [representation]
775  else:
776  representations = representation
777 
778  if columnmapping is None:
779  columnmapping = {}
780  columnmapping["Protein1"] = 0
781  columnmapping["Protein2"] = 1
782  columnmapping["Residue1"] = 2
783  columnmapping["Residue2"] = 3
784  columnmapping["IDScore"] = 4
785  columnmapping["XLUniqueID"] = 5
786 
787  if csvfile:
788  # in case the file is a cvs file
789  # columnmapping will contain the field names
790  # that compare in the first line of the cvs file
791  db = IMP.pmi.tools.get_db_from_csv(restraints_file)
792  else:
793  db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
794 
795  indb = open("included." + filelabel + ".xl.db", "w")
796  exdb = open("excluded." + filelabel + ".xl.db", "w")
797  midb = open("missing." + filelabel + ".xl.db", "w")
798 
799  self.m = representations[0].prot.get_model()
800  self.marginal = marginal
801  self.rs = IMP.RestraintSet(self.m, 'data')
802  if not self.marginal:
803  self.rspsi = IMP.RestraintSet(self.m, 'prior_psi')
804  self.rssig = IMP.RestraintSet(self.m, 'prior_sigmas')
805  self.rslin = IMP.RestraintSet(self.m, 'prior_linear')
806  self.rslen = IMP.RestraintSet(self.m, 'prior_length')
807 
808  self.label = label
809  self.pairs = []
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
816 
817  # isd_map is a dictionary/map that is used to determine the psi
818  # parameter from identity scores (such as ID-Score, or FDR)
819  if ids_map is None:
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)
823  else:
824  self.ids_map = ids_map
825 
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)
831  else:
832  self.radius_map = radius_map
833 
834  self.outputlevel = "low"
835 
836  # small linear contribution for long range
837  self.linear = IMP.core.Linear(0, 0.0)
838  self.linear.set_slope(slope)
839  dps2 = IMP.core.DistancePairScore(self.linear)
840 
841  protein1 = columnmapping["Protein1"]
842  protein2 = columnmapping["Protein2"]
843  residue1 = columnmapping["Residue1"]
844  residue2 = columnmapping["Residue2"]
845  idscore = columnmapping["IDScore"]
846  try:
847  xluniqueid = columnmapping["XLUniqueID"]
848  except:
849  xluniqueid = None
850 
851  restraints = []
852 
853  # we need this dictionary to create ambiguity (i.e., multistate)
854  # if one id is already present in the dictionary, add the term to the
855  # corresponding already geenrated restraint
856 
857  uniqueid_restraints_map = {}
858 
859  for nxl, entry in enumerate(db):
860 
861  if not jackknifing is None:
862 
863  # to be implemented
864  # the problem is that in the replica exchange
865  # you have to broadcast the same restraints to every
866  # replica
867 
868  raise NotImplementedError("jackknifing not yet implemented")
869 
870  if not csvfile:
871  tokens = entry.split()
872  if len(tokens)==0:
873  continue
874 
875  # skip character
876  if (tokens[0] == "#"):
877  continue
878  try:
879  r1 = int(tokens[residue1])
880  c1 = tokens[protein1]
881  r2 = int(tokens[residue2])
882  c2 = tokens[protein2]
883 
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]
887 
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]
891 
892  if idscore is None:
893  ids = 1.0
894  else:
895  ids = float(tokens[idscore])
896  if xluniqueid is None:
897  xlid = str(nxl)
898  else:
899  xlid = tokens[xluniqueid]
900  except:
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")
908  continue
909 
910  else:
911  if filters is not None:
912  if eval(IMP.pmi.tools.cross_link_db_filter_parser(filters)) == False:
913  exdb.write(str(entry) + "\n")
914  continue
915 
916  try:
917  r1 = int(entry[residue1])
918  c1 = entry[protein1]
919  r2 = int(entry[residue2])
920  c2 = entry[protein2]
921 
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]
925 
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]
929 
930  if idscore is None:
931  ids = 1.0
932  else:
933  try:
934  ids = float(entry[idscore])
935  except ValueError:
936  ids = entry[idscore]
937  if xluniqueid is None:
938  xlid = str(nxl)
939  else:
940  xlid = entry[xluniqueid]
941 
942  except:
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")
950  continue
951 
952  for nstate, r in enumerate(representations):
953  # loop over every state
954  ps1 = IMP.pmi.tools.select(
955  r,
956  resolution=resolution,
957  name=c1,
958  name_is_ambiguous=False,
959  residue=r1)
960  ps2 = IMP.pmi.tools.select(
961  r,
962  resolution=resolution,
963  name=c2,
964  name_is_ambiguous=False,
965  residue=r2)
966 
967  if len(ps1) > 1:
968  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
969  elif len(ps1) == 0:
970  print("ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
971  midb.write(str(entry) + "\n")
972  continue
973 
974  if len(ps2) > 1:
975  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
976  elif len(ps2) == 0:
977  print("ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
978  midb.write(str(entry) + "\n")
979  continue
980 
981  p1 = ps1[0]
982  p2 = ps2[0]
983 
984  # remove in the future!!!
985  if p1 == p2:
986  continue
987 
988  if xlid in uniqueid_restraints_map:
989  print("getting a crosslink restraint from id %s" % str(xlid))
990  dr = uniqueid_restraints_map[xlid]
991 
992  else:
993  if not self.marginal:
994  if not self.samplelength:
995  print("generating a new crosslink restraint")
997  self.m,
998  length,
999  inner_slope)
1000  else:
1001  # this will create a xl length particle that will
1002  # be sampled
1003  self.create_length()
1005  self.m,
1006  self.length)
1007 
1008  else:
1009  if not self.samplelength:
1010  dr = IMP.isd_emxl.CrossLinkMSMarginalRestraint(
1011  self.m,
1012  length)
1013  else:
1014  # this will create a xl length particle that will
1015  # be sampled
1016  self.create_length()
1017  dr = IMP.isd_emxl.CrossLinkMSMarginalRestraint(
1018  self.m,
1019  self.length)
1020 
1021  uniqueid_restraints_map[xlid] = dr
1022 
1023  mappedr1 = self.radius_map.get_map_element(
1024  IMP.pmi.Uncertainty(p1).get_uncertainty())
1025  sigma1 = self.get_sigma(mappedr1)[0]
1026  mappedr2 = self.radius_map.get_map_element(
1027  IMP.pmi.Uncertainty(p2).get_uncertainty())
1028  sigma2 = self.get_sigma(mappedr2)[0]
1029 
1030  if not self.marginal:
1031  psival = self.ids_map.get_map_element(ids)
1032  psi = self.get_psi(psival)[0]
1033 
1034 
1035  p1i = p1.get_particle_index()
1036  p2i = p2.get_particle_index()
1037  s1i = sigma1.get_particle().get_index()
1038  s2i = sigma2.get_particle().get_index()
1039 
1040  #print nstate, p1i, p2i, p1.get_name(), p2.get_name()
1041 
1042  if not self.marginal:
1043  psii = psi.get_particle().get_index()
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")
1051  else:
1052  psival = None
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")
1060 
1061  indb.write(str(entry) + "\n")
1062 
1063  restraints.append(dr)
1064 
1065  # check if the two residues belong to the same rigid body
1068  IMP.core.RigidMember(p1).get_rigid_body() ==
1069  IMP.core.RigidMember(p2).get_rigid_body()):
1070  xlattribute = "intrarb"
1071  else:
1072  xlattribute = "interrb"
1073 
1074  if csvfile:
1075  if not attributes_for_label is None:
1076  for a in attributes_for_label:
1077  xlattribute = xlattribute + "_" + str(entry[a])
1078 
1079  xlattribute = xlattribute + "-State:" + str(nstate)
1080 
1081  dr.set_name(
1082  xlattribute + "-" + c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2) + "_" + self.label)
1083 
1084  pr = IMP.core.PairRestraint(dps2, IMP.ParticlePair(p1, p2))
1085  pr.set_name(
1086  xlattribute + "-" + c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2) + "_" + self.label)
1087  self.rslin.add_restraint(pr)
1088 
1089 
1090  self.pairs.append(
1091  (p1,
1092  p2,
1093  dr,
1094  r1,
1095  c1,
1096  r2,
1097  c2,
1098  xlattribute,
1099  mappedr1,
1100  mappedr2,
1101  psival,
1102  xlid))
1103 
1104  lw = IMP.isd.LogWrapper(restraints,1.0)
1105  self.rs.add_restraint(lw)
1106 
1107 
1108  def set_slope_linear_term(self, slope):
1109  self.linear.set_slope(slope)
1110 
1111  def set_label(self, label):
1112  self.label = label
1113 
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)
1121 
1122  def get_hierarchies(self):
1123  return self.prot
1124 
1125  def get_restraint_sets(self):
1126  return self.rs
1127 
1128  def get_restraint(self):
1129  return self.rs
1130 
1131  def get_restraint_for_rmf(self):
1132  return self.rslin
1133 
1134  def get_restraints(self):
1135  rlist = []
1136  for r in self.rs.get_restraints():
1137  rlist.append(IMP.core.PairRestraint.get_from(r))
1138  return rlist
1139 
1140  def get_particle_pairs(self):
1141  ppairs = []
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))
1146  return ppairs
1147 
1148  def set_output_level(self, level="low"):
1149  # this might be "low" or "high"
1150  self.outputlevel = level
1151 
1152  def set_psi_is_sampled(self, is_sampled=True):
1153  self.psi_is_sampled = is_sampled
1154 
1155  def set_sigma_is_sampled(self, is_sampled=True):
1156  self.sigma_is_sampled = is_sampled
1157 
1158  def get_output(self):
1159  # content of the crosslink database pairs
1160  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
1161  self.m.update()
1162 
1163  output = {}
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)):
1175 
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]
1188 
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)))
1194  d0 = IMP.core.XYZ(p0)
1195  d1 = IMP.core.XYZ(p1)
1196  output["ISDCrossLinkMS_Distance_" +
1197  label + "_" + self.label] = str(IMP.core.get_distance(d0, d1))
1198 
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())
1203 
1204  for resolution in self.sigma_dictionary:
1205  output["ISDCrossLinkMS_Sigma_" +
1206  str(resolution) + "_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
1207 
1208  if self.samplelength:
1209  output["ISDCrossLinkMS_Length_" +
1210  str(resolution) + "_" + self.label] = str(self.length.get_scale())
1211 
1212  return output
1213 
1214  def get_particles_to_sample(self):
1215  ps = {}
1216 
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])
1222 
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])
1228 
1229  if self.samplelength:
1230  if self.lengthissampled:
1231  ps["Nuisances_ISDCrossLinkMS_Length_" +
1232  self.label] = ([self.length], self.lengthtrans)
1233 
1234  return ps
1235 
1236 #
1237 class CysteineCrossLinkRestraint(object):
1238 
1239  import IMP.isd
1240  import IMP.pmi.tools
1241 
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),
1247  betaissampled=True,
1248  sigmaissampled=False,
1249  weightissampled=True,
1250  epsilonissampled=True
1251  ):
1252  # the file must have residue1 chain1 residue2 chain2 fractionvalue epsilonname
1253  # epsilonname is a name for the epsilon particle that must be used for that particular
1254  # residue pair, eg, "Epsilon-Intra-Solvent", or
1255  # "Epsilon-Solvent-Membrane", etc.
1256 
1257  self.representations = representations
1258  self.m = self.representations[0].prot.get_model()
1259  self.rs = IMP.RestraintSet(self.m, 'Cysteine_Crosslink')
1260  self.cbeta = cbeta
1261  self.epsilonmaxtrans = 0.01
1262  self.sigmamaxtrans = 0.1
1263  self.betamaxtrans = 0.01
1264  self.weightmaxtrans = 0.1
1265  self.label = "None"
1266  self.outputlevel = "low"
1267  self.betaissampled = betaissampled
1268  self.sigmaissampled = sigmaissampled
1269  self.weightissampled = weightissampled
1270  self.epsilonissampled = epsilonissampled
1271 
1272  betalower = betatuple[0]
1273  betaupper = betatuple[1]
1274  beta = 0.04
1275  sigma = 10.0
1276  betangrid = 30
1277  crossdataprior = 1
1278 
1279  # beta
1280  self.beta = IMP.pmi.tools.SetupNuisance(
1281  self.m,
1282  beta,
1283  betalower,
1284  betaupper,
1285  betaissampled).get_particle(
1286  )
1287  # sigma
1288  self.sigma = IMP.pmi.tools.SetupNuisance(
1289  self.m,
1290  sigma,
1291  sigmatuple[0],
1292  sigmatuple[1],
1293  sigmaissampled).get_particle()
1294  # population particle
1295  self.weight = IMP.pmi.tools.SetupWeight(
1296  self.m,
1297  weightissampled).get_particle(
1298  )
1299 
1300  # read the file
1301  fl = IMP.pmi.tools.open_file_or_inline_text(filename)
1302  # determine the upperlimit for epsilon
1303  # and also how many epsilon are needed
1304  self.epsilons = {}
1305  data = []
1306  for line in fl:
1307  t = line.split()
1308  if t[0][0] == "#":
1309  continue
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]))
1313  else:
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()
1318  if up < low:
1319  self.epsilons[t[5]].set_upper(low)
1320 
1321  data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
1322 
1323  # create CrossLinkData
1324  if not self.cbeta:
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)
1328  else:
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)
1332 
1333  # create grids needed by CysteineCrossLinkData
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)
1337 
1338  for d in data:
1339  print("--------------")
1340  print("CysteineCrossLink: attempting to create a restraint " + str(d))
1341  resid1 = d[0]
1342  chain1 = d[1]
1343  resid2 = d[2]
1344  chain2 = d[3]
1345  fexp = d[4]
1346  epslabel = d[5]
1347 
1348  # CysteineCrossLinkData
1349 
1351  fexp,
1352  fmod_grid,
1353  omega2_grid,
1354  beta_grid)
1355 
1357  self.beta,
1358  self.sigma,
1359  self.epsilons[epslabel],
1360  self.weight,
1361  crossdata,
1362  ccldata)
1363 
1364  failed = False
1365  for i, representation in enumerate(self.representations):
1366 
1367  if not self.cbeta:
1368  p1 = None
1369  p2 = None
1370 
1371  p1 = IMP.pmi.tools.select(representation,
1372  resolution=1, name=chain1,
1373  name_is_ambiguous=False, residue=resid1)[0]
1374 
1375  if p1 is None:
1376  failed = True
1377 
1378  p2 = IMP.pmi.tools.select(representation,
1379  resolution=1, name=chain2,
1380  name_is_ambiguous=False, residue=resid2)[0]
1381 
1382  if p2 is None:
1383  failed = True
1384 
1385  else:
1386  # use cbetas
1387  p1 = []
1388  p2 = []
1389  for t in range(-1, 2):
1390  p = IMP.pmi.tools.select(representation,
1391  resolution=1, name=chain1,
1392  name_is_ambiguous=False, residue=resid1 + t)
1393 
1394  if len(p) == 1:
1395  p1 += p
1396  else:
1397  failed = True
1398  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
1399 
1400  p = IMP.pmi.tools.select(representation,
1401  resolution=1, name=chain2,
1402  name_is_ambiguous=False, residue=resid2 + t)
1403 
1404  if len(p) == 1:
1405  p2 += p
1406  else:
1407  failed = True
1408  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
1409 
1410  if not self.cbeta:
1411  if (p1 is not None and p2 is not None):
1412  ccl.add_contribution(p1, p2)
1413  d1 = IMP.core.XYZ(p1)
1414  d2 = IMP.core.XYZ(p2)
1415 
1416  print("Distance_" + str(resid1) + "_" + chain1 + ":" + str(resid2) + "_" + chain2, IMP.core.get_distance(d1, d2))
1417 
1418  else:
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()
1426 
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))
1430 
1431  ccl.add_contribution(p1, p2)
1432 
1433  if not failed:
1434  self.rs.add_restraint(ccl)
1435  ccl.set_name("CysteineCrossLink_" + str(resid1)
1436  + "_" + chain1 + ":" + str(resid2) + "_" + chain2)
1437 
1438  def set_label(self, label):
1439  self.label = label
1440 
1441  def add_to_model(self):
1442  self.m.add_restraint(self.rs)
1443 
1444  def get_particles_to_sample(self):
1445  ps = {}
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)
1459  return ps
1460 
1461  def set_output_level(self, level="low"):
1462  # this might be "low" or "high"
1463  self.outputlevel = level
1464 
1465  def get_restraint(self):
1466  return self.rs
1467 
1468  def get_sigma(self):
1469  return self.sigma
1470 
1471  def get_output(self):
1472  self.m.update()
1473  output = {}
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))
1487 
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(
1495  rst).get_name(
1496  ) + "_"
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]
1503 
1504  return output
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)
Miscellaneous utilities.
Definition: tools.py:1
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.
Definition: Uncertainty.h:24
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.
Definition: XYZR.h:87
Uniform distribution with harmonic boundaries.
Definition: UniformPrior.h:20
A class to store an fixed array of same-typed values.
Definition: Array.h:33
static XYZR setup_particle(kernel::Model *m, ParticleIndex pi)
Definition: XYZR.h:48
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.
Definition: XYZ.h:30
Class to handle individual model particles.
this restraint allows ambiguous crosslinking between multiple copies it is a variant of the Simplifie...
Definition: crosslinking.py:98
Linear function
Definition: Linear.h:19
Classes for writing output files and processing them.
Definition: output.py:1
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)
Definition: rigid_bodies.h:479
VectorD< 3 > Vector3D
Definition: VectorD.h:395
Calculate the -Log of a list of restraints.
Definition: LogWrapper.h:19
Applies a PairScore to a Pair.
Definition: PairRestraint.h:29
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:633
def cross_link_db_filter_parser
example '"{ID_Score}" > 28 AND "{Sample}" == "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample...
Definition: tools.py:425
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:62
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.
Definition: XYZR.h:27