1 """@namespace IMP.pmi.restraints.crosslinking
2 Restraints for handling crosslinking data.
5 from __future__
import print_function
17 from collections
import defaultdict
25 """Setup cross-link distance restraints from mass spectrometry data.
26 The noise in the data and the structural uncertainty of cross-linked amino-acids
27 is inferred using Bayes theory of probability
28 @note Wraps an IMP::isd::CrossLinkMSRestraint
30 def __init__(self, representation=None,
32 CrossLinkDataBase=
None,
38 attributes_for_label=
None,
41 @param representation DEPRECATED The IMP.pmi.representation.Representation
42 object that contain the molecular system
43 @param root_hier The canonical hierarchy containing all the states
44 @param CrossLinkDataBase The IMP.pmi.io.crosslink.CrossLinkDataBase
45 object that contains the cross-link dataset
46 @param length maximal cross-linker length (including the residue sidechains)
47 @param resolution what representation resolution should the cross-link
48 restraint be applied to.
49 @param slope The slope of a distance-linear scoring function that
50 funnels the score when the particles are
51 too far away. Suggested value 0.02.
52 @param label the extra text to label the restraint so that it is
53 searchable in the output
54 @param filelabel automatically generated file containing missing/included/excluded
55 cross-links will be labeled using this text
56 @param attributes_for_label
57 @param weight Weight of restraint
61 if representation
is not None:
63 if type(representation) != list:
64 representations = [representation]
66 representations = representation
67 model = representations[0].prot.get_model()
68 elif root_hier
is not None:
70 model = root_hier.get_model()
72 raise Exception(
"You must pass either representation or root_hier")
74 super(CrossLinkingMassSpectrometryRestraint, self).
__init__(
75 model, weight=weight, label=label)
77 if CrossLinkDataBase
is None:
78 raise Exception(
"You must pass a CrossLinkDataBase")
79 if not isinstance(CrossLinkDataBase,IMP.pmi.io.crosslink.CrossLinkDataBase):
80 raise TypeError(
"CrossLinkingMassSpectrometryRestraint: CrossLinkDataBase should be an IMP.pmi.io.crosslink.CrossLinkDataBase object")
81 self.CrossLinkDataBase = CrossLinkDataBase
83 if resolution==0
or resolution
is None:
84 raise Exception(
"You must pass a resolution and it can't be zero")
86 indb = open(
"included." + filelabel +
".xl.db",
"w")
87 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
88 midb = open(
"missing." + filelabel +
".xl.db",
"w")
90 self.rs.set_name(self.rs.get_name() +
"_Data")
91 self.rspsi = self._create_restraint_set(
"PriorPsi")
92 self.rssig = self._create_restraint_set(
"PriorSig")
93 self.rslin = self._create_restraint_set(
"Linear")
97 self.linear.set_slope(0.0)
100 self.psi_is_sampled =
True
101 self.sigma_is_sampled =
True
102 self.psi_dictionary={}
103 self.sigma_dictionary={}
105 self.outputlevel =
"low"
109 xl_groups = [p.get_cross_link_group(self)
110 for p, state
in IMP.pmi.tools._all_protocol_outputs(
111 representations, root_hier)]
115 copies_to_add = defaultdict(int)
116 print(
'gathering copies')
117 for xlid
in self.CrossLinkDataBase.xlid_iterator():
118 for xl
in self.CrossLinkDataBase[xlid]:
119 r1 = xl[self.CrossLinkDataBase.residue1_key]
120 c1 = xl[self.CrossLinkDataBase.protein1_key]
121 r2 = xl[self.CrossLinkDataBase.residue2_key]
122 c2 = xl[self.CrossLinkDataBase.protein2_key]
123 for c,r
in ((c1,r1),(c2,r2)):
124 if c
in copies_to_add:
130 resolution=resolution).get_selected_particles()
132 copies_to_add[c] = len(sel)-1
134 for molname
in copies_to_add:
135 if copies_to_add[molname]==0:
138 self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein1_key,molname+
'.0',fo1)
140 self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein2_key,molname+
'.0',fo2)
141 for ncopy
in range(copies_to_add[molname]):
142 self.CrossLinkDataBase.clone_protein(
'%s.0'%molname,
'%s.%i'%(molname,ncopy+1))
143 print(
'done pmi2 prelims')
145 for xlid
in self.CrossLinkDataBase.xlid_iterator():
146 new_contribution=
True
147 for xl
in self.CrossLinkDataBase[xlid]:
149 r1 = xl[self.CrossLinkDataBase.residue1_key]
150 c1 = xl[self.CrossLinkDataBase.protein1_key]
151 r2 = xl[self.CrossLinkDataBase.residue2_key]
152 c2 = xl[self.CrossLinkDataBase.protein2_key]
155 ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
158 zip(IMP.pmi.tools._all_protocol_outputs(
159 representations, root_hier),
163 iterlist = range(len(IMP.atom.get_by_type(root_hier,IMP.atom.STATE_TYPE)))
165 iterlist = representations
166 for nstate, r
in enumerate(iterlist):
168 xl[self.CrossLinkDataBase.state_key]=nstate
169 xl[self.CrossLinkDataBase.data_set_name_key]=self.label
177 name1,copy1 = c1.split(
'.')
179 name2,copy2 = c2.split(
'.')
183 copy_index=int(copy1),
185 resolution=resolution).get_selected_particles()
189 copy_index=int(copy2),
191 resolution=resolution).get_selected_particles()
198 resolution=resolution,
200 name_is_ambiguous=
False,
204 resolution=resolution,
206 name_is_ambiguous=
False,
210 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
212 print(
"CrossLinkingMassSpectrometryRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
213 midb.write(str(xl) +
"\n")
217 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
219 print(
"CrossLinkingMassSpectrometryRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
220 midb.write(str(xl) +
"\n")
226 if p1 == p2
and r1 == r2:
227 print(
"CrossLinkingMassSpectrometryRestraint: WARNING> same particle and same residue, skipping cross-link")
231 print(
"generating a new crosslink restraint")
232 new_contribution=
False
237 restraints.append(dr)
240 if self.CrossLinkDataBase.sigma1_key
not in xl.keys():
242 xl[self.CrossLinkDataBase.sigma1_key]=sigma1name
244 sigma1name=xl[self.CrossLinkDataBase.sigma1_key]
247 if self.CrossLinkDataBase.sigma2_key
not in xl.keys():
249 xl[self.CrossLinkDataBase.sigma2_key]=sigma2name
251 sigma2name=xl[self.CrossLinkDataBase.sigma2_key]
254 if self.CrossLinkDataBase.psi_key
not in xl.keys():
256 xl[self.CrossLinkDataBase.psi_key]=psiname
258 psiname=xl[self.CrossLinkDataBase.psi_key]
263 p1i = p1.get_particle_index()
265 p2i = p2.get_particle_index()
268 xl[
"Particle_sigma1"]=sigma1
270 xl[
"Particle_sigma2"]=sigma2
272 xl[
"Particle_psi"]=psi
274 dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
277 print(
"--------------")
278 print(
"CrossLinkingMassSpectrometryRestraint: generating cross-link restraint between")
279 print(
"CrossLinkingMassSpectrometryRestraint: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
280 print(
"CrossLinkingMassSpectrometryRestraint: with sigma1 %s sigma2 %s psi %s" % (sigma1name, sigma2name, psiname))
281 print(
"CrossLinkingMassSpectrometryRestraint: between particles %s and %s" % (p1.get_name(), p2.get_name()))
282 print(
"==========================================\n")
283 for p, ex_xl
in zip(IMP.pmi.tools._all_protocol_outputs(
284 representations, root_hier),
286 p[0].add_cross_link(p[1], ex_xl[0], p1, p2, length,
287 sigma1, sigma2, psi, ex_xl[1])
294 xl[
"IntraRigidBody"]=
True
296 xl[
"IntraRigidBody"]=
False
298 xl_label=self.CrossLinkDataBase.get_short_cross_link_string(xl)
299 xl[
"ShortLabel"]=xl_label
300 dr.set_name(xl_label)
304 pr.set_name(xl_label)
305 self.rslin.add_restraint(pr)
307 self.xl_list.append(xl)
309 indb.write(str(xl) +
"\n")
311 if len(self.xl_list) == 0:
312 raise SystemError(
"CrossLinkingMassSpectrometryRestraint: no crosslink was constructed")
313 self.xl_restraints = restraints
315 self.rs.add_restraint(lw)
317 def __set_dataset(self, ds):
318 self.CrossLinkDataBase.dataset = ds
319 dataset = property(
lambda self: self.CrossLinkDataBase.dataset,
323 """ get the hierarchy """
327 """ get the restraint set """
331 """ get the restraints in a list """
332 return self.xl_restraints
335 """ get the dummy restraints to be displayed in the rmf file """
339 """ Get a list of tuples containing the particle pairs """
341 for i
in range(len(self.pairs)):
342 p0 = self.pairs[i][0]
343 p1 = self.pairs[i][1]
344 ppairs.append((p0, p1))
348 """ Set the output level of the output """
349 self.outputlevel = level
352 """ Switch on/off the sampling of psi particles """
353 self.psi_is_sampled = is_sampled
356 """ Switch on/off the sampling of sigma particles """
357 self.sigma_is_sampled = is_sampled
361 """ This is called internally. Creates a nuisance
362 on the structural uncertainty """
363 if name
in self.sigma_dictionary:
364 return self.sigma_dictionary[name][0]
367 sigmaminnuis = 0.0000001
368 sigmamaxnuis = 1000.0
372 sigma = IMP.pmi.tools.SetupNuisance(self.model, sigmainit,
373 sigmaminnuis, sigmamaxnuis, self.sigma_is_sampled).get_particle()
374 self.sigma_dictionary[name] = (
377 self.sigma_is_sampled)
378 self.rssig.add_restraint(
388 """ This is called internally. Creates a nuisance
389 on the data uncertainty """
390 if name
in self.psi_dictionary:
391 return self.psi_dictionary[name][0]
394 psiminnuis = 0.0000001
395 psimaxnuis = 0.4999999
399 psi = IMP.pmi.tools.SetupNuisance(self.model, psiinit,
400 psiminnuis, psimaxnuis,
401 self.psi_is_sampled).get_particle()
402 self.psi_dictionary[name] = (
407 self.rspsi.add_restraint(
419 """ Get the output of the restraint to be used by the IMP.pmi.output object"""
420 output = super(CrossLinkingMassSpectrometryRestraint, self).
get_output()
422 for xl
in self.xl_list:
424 xl_label=xl[
"ShortLabel"]
428 output[
"CrossLinkingMassSpectrometryRestraint_Score_" +
429 xl_label] = str(-log(ln.unprotected_evaluate(
None)))
433 output[
"CrossLinkingMassSpectrometryRestraint_Distance_" +
437 for psiname
in self.psi_dictionary:
438 output[
"CrossLinkingMassSpectrometryRestraint_Psi_" +
439 str(psiname) + self._label_suffix] = str(
440 self.psi_dictionary[psiname][0].get_scale())
442 for sigmaname
in self.sigma_dictionary:
443 output[
"CrossLinkingMassSpectrometryRestraint_Sigma_" +
444 str(sigmaname) + self._label_suffix] = str(
445 self.sigma_dictionary[sigmaname][0].get_scale())
451 """ Get all need data to construct a mover in IMP.pmi.dof class"""
453 if self.sigma_is_sampled:
454 for sigmaname
in self.sigma_dictionary:
455 mover_name=
"Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" + str(sigmaname) +
"_" + self.label
456 particle=self.sigma_dictionary[sigmaname][0]
457 maxstep=(self.sigma_dictionary[sigmaname][1])
460 mv.set_name(mover_name)
463 if self.psi_is_sampled:
464 for psiname
in self.psi_dictionary:
465 mover_name=
"Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" + str(psiname) +
"_" + self.label
466 particle=self.psi_dictionary[psiname][0]
467 maxstep=(self.psi_dictionary[psiname][1])
470 mv.set_name(mover_name)
477 """ Get the particles to be sampled by the IMP.pmi.sampler object """
479 if self.sigma_is_sampled:
480 for sigmaname
in self.sigma_dictionary:
481 ps[
"Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" +
482 str(sigmaname) + self._label_suffix] =\
483 ([self.sigma_dictionary[sigmaname][0]],
484 self.sigma_dictionary[sigmaname][1])
486 if self.psi_is_sampled:
487 for psiname
in self.psi_dictionary:
488 ps[
"Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
489 str(psiname) + self._label_suffix] =\
490 ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
496 """Setup cross-link distance restraints at atomic level
497 The "atomic" aspect is that it models the particle uncertainty with a Gaussian
498 The noise in the data and the structural uncertainty of cross-linked amino-acids
499 is inferred using Bayes' theory of probability
500 @note Wraps an IMP::isd::AtomicCrossLinkMSRestraint
510 nuisances_are_optimized=
True,
517 Automatically creates one "sigma" per crosslinked residue and one "psis" per pair.
518 Other nuisance options are available.
519 @note Will return an error if the data+extra_sel don't specify two particles per XL pair.
520 @param root_hier The root hierarchy on which you'll do selection
521 @param xldb CrossLinkDataBase object
522 @param atom_type Can either be "NZ" or "CA"
523 @param length The XL linker length
524 @param slope Linear term to add to the restraint function to help when far away
525 @param nstates The number of states to model. Defaults to the number of states in root.
526 @param label The output label for the restraint
527 @param nuisances_are_optimized Whether to optimize nuisances
528 @param sigma_init The initial value for all the sigmas
529 @param psi_init The initial value for all the psis
530 @param one_psi Use a single psi for all restraints (if False, creates one per XL)
531 @param filelabel automatically generated file containing missing/included/excluded
532 cross-links will be labeled using this text
533 @param weight Weight of restraint
538 self.root = root_hier
539 rname =
"AtomicXLRestraint"
540 super(AtomicCrossLinkMSRestraint, self).
__init__(
541 self.root.get_model(), name=
"AtomicXLRestraint", label=label,
545 self.sigma_is_sampled = nuisances_are_optimized
546 self.psi_is_sampled = nuisances_are_optimized
548 if atom_type
not in(
"NZ",
"CA"):
549 raise Exception(
"AtomicXLRestraint: atom_type must be NZ or CA")
550 self.atom_type = atom_type
551 self.nstates = nstates
553 self.nstates = len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE))
554 elif nstates!=len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE)):
555 print(
"Warning: nstates is not the same as the number of states in root")
557 self.rs_psi = self._create_restraint_set(
"psi")
558 self.rs_sig = self._create_restraint_set(
"sigma")
559 self.rs_lin = self._create_restraint_set(
"linear")
561 self.psi_dictionary = {}
562 self.sigma_dictionary = {}
564 self.particles=defaultdict(set)
565 self.one_psi = one_psi
566 self.bonded_pairs = []
568 print(
'creating a single psi for all XLs')
570 print(
'creating one psi for each XL')
573 if filelabel
is not None:
574 indb = open(
"included." + filelabel +
".xl.db",
"w")
575 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
576 midb = open(
"missing." + filelabel +
".xl.db",
"w")
582 # read ahead to get the number of XL's per residue
583 num_xls_per_res=defaultdict(int)
584 for unique_id in data:
585 for nstate in range(self.nstates):
586 for xl in data[unique_id]:
587 num_xls_per_res[str(xl.r1)]+=1
588 num_xls_per_res[str(xl.r2)]+=1
590 # Will setup two sigmas based on promiscuity of the residue
592 self.sig_low = setup_nuisance(self.m,self.rs_nuis,init_val=sigma_init,min_val=1.0,
593 max_val=100.0,is_opt=self.nuis_opt)
594 self.sig_high = setup_nuisance(self.m,self.rs_nuis,init_val=sigma_init,min_val=1.0,
595 max_val=100.0,is_opt=self.nuis_opt)
597 self._create_sigma(
'sigma',sigma_init)
599 self._create_psi(
'psi',psi_init)
601 for xlid
in self.xldb.xlid_iterator():
602 self._create_psi(xlid,psi_init)
606 for xlid
in self.xldb.xlid_iterator():
609 psip = self.psi_dictionary[
'psi'][0].get_particle_index()
611 psip = self.psi_dictionary[unique_id][0].get_particle_index()
620 for nstate
in self.nstates:
621 for xl
in self.xldb[xlid]:
622 r1 = xl[self.xldb.residue1_key]
623 c1 = xl[self.xldb.protein1_key].strip()
624 r2 = xl[self.xldb.residue2_key]
625 c2 = xl[self.xldb.protein2_key].strip()
632 residue_index=r1).get_selected_particles()
637 residue_index=r2).get_selected_particles()
639 print(
"AtomicXLRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
640 if filelabel
is not None:
641 midb.write(str(xl) +
"\n")
645 print(
"AtomicXLRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
646 if filelabel
is not None:
647 midb.write(str(xl) +
"\n")
653 num1=num_xls_per_res[str(xl.r1)]
654 num2=num_xls_per_res[str(xl.r2)]
655 if num1<sig_threshold:
659 if num2<sig_threshold:
664 sig1 = self.sigma_dictionary[
'sigma'][0]
665 sig2 = self.sigma_dictionary[
'sigma'][0]
668 for p1,p2
in itertools.product(ps1,ps2):
671 self.particles[nstate]|=set((p1,p2))
672 r.add_contribution([p1.get_index(),p2.get_index()],
673 [sig1.get_particle_index(),sig2.get_particle_index()])
676 if num_contributions>0:
677 print(
'XL:',xlid,
'num contributions:',num_contributions)
680 raise Exception(
"You didn't create any XL restraints")
681 print(
'created',len(xlrs),
'XL restraints')
682 rname = self.rs.get_name()
684 self.rs.set_name(rname)
685 self.rs.set_weight(self.weight)
686 self.restraint_sets = [self.rs] + self.restraint_sets[1:]
688 def get_hierarchy(self):
691 def _create_sigma(self, name,sigmainit):
692 """ This is called internally. Creates a nuisance
693 on the structural uncertainty """
694 if name
in self.sigma_dictionary:
695 return self.sigma_dictionary[name][0]
697 sigmaminnuis = 0.0000001
698 sigmamaxnuis = 1000.0
702 sigma = IMP.pmi.tools.SetupNuisance(self.model,
706 self.sigma_is_sampled).get_particle()
707 self.sigma_dictionary[name] = (
710 self.sigma_is_sampled)
711 self.rs_sig.add_restraint(
720 def _create_psi(self, name,psiinit):
721 """ This is called internally. Creates a nuisance
722 on the data uncertainty """
723 if name
in self.psi_dictionary:
724 return self.psi_dictionary[name][0]
726 psiminnuis = 0.0000001
727 psimaxnuis = 0.4999999
731 psi = IMP.pmi.tools.SetupNuisance(self.model,
735 self.psi_is_sampled).get_particle()
736 self.psi_dictionary[name] = (
741 self.rs_psi.add_restraint(
753 """ create dummy harmonic restraints for each XL but don't add to model
754 Makes it easy to see each contribution to each XL in RMF
756 class MyGetRestraint(object):
765 for nxl
in range(self.get_number_of_restraints()):
766 xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
768 for ncontr
in range(xl.get_number_of_contributions()):
769 ps=xl.get_contribution(ncontr)
771 'xl%i_contr%i'%(nxl,ncontr))
773 dummy_rs.append(MyGetRestraint(rs))
778 """ Get the particles to be sampled by the IMP.pmi.sampler object """
780 if self.sigma_is_sampled:
781 for sigmaname
in self.sigma_dictionary:
782 ps[
"Nuisances_AtomicCrossLinkingMSRestraint_Sigma_" +
783 str(sigmaname) + self._label_suffix] = \
784 ([self.sigma_dictionary[sigmaname][0]],
785 self.sigma_dictionary[sigmaname][1])
786 if self.psi_is_sampled:
787 for psiname
in self.psi_dictionary:
788 ps[
"Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
789 str(psiname) + self._label_suffix] =\
790 ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
793 def get_bonded_pairs(self):
794 return self.bonded_pairs
796 def get_number_of_restraints(self):
797 return self.rs.get_number_of_restraints()
800 return 'XL restraint with '+str(len(self.rs.get_restraint(0).get_number_of_restraints())) \
804 """Read a stat file and load all the sigmas.
805 This is potentially quite stupid.
806 It's also a hack since the sigmas should be stored in the RMF file.
807 Also, requires one sigma and one psi for ALL XLs.
810 sig_val = float(subprocess.check_output([
"process_output.py",
"-f",in_fn,
811 "-s",
"AtomicXLRestraint_sigma"]).split(
'\n>')[1+nframe])
812 psi_val = float(subprocess.check_output([
"process_output.py",
"-f",in_fn,
813 "-s",
"AtomicXLRestraint_psi"]).split(
'\n>')[1+nframe])
814 for nxl
in range(self.rs.get_number_of_restraints()):
815 xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
818 for contr
in range(xl.get_number_of_contributions()):
819 sig1,sig2=xl.get_contribution_sigmas(contr)
822 print(
'loaded nuisances from file')
825 max_prob_for_violation=0.1,
826 min_dist_for_violation=1e9,
828 limit_to_chains=
None,
830 """Create CMM files, one for each state, of all crosslinks.
831 will draw in GREEN if non-violated in all states (or if only one state)
832 will draw in PURPLE if non-violated only in a subset of states (draws nothing elsewhere)
833 will draw in RED in ALL states if all violated
834 (if only one state, you'll only see green and red)
836 @param out_prefix Output xlink files prefix
837 @param max_prob_for_violation It's a violation if the probability is below this
838 @param min_dist_for_violation It's a violation if the min dist is above this
839 @param coarsen Use CA positions
840 @param limit_to_chains Try to visualize just these chains
841 @param exclude_to_chains Try to NOT visualize these chains
843 print(
'going to calculate violations and plot CMM files')
845 all_dists = [s[
"low_dist"]
for s
in all_stats]
851 cmds = defaultdict(set)
852 for nstate
in range(self.nstates):
853 outf=open(out_prefix+str(nstate)+
'.cmm',
'w')
854 outf.write(
'<marker_set name="xlinks_state%i"> \n' % nstate)
857 print(
'will limit to',limit_to_chains)
858 print(
'will exclude',exclude_chains)
863 for nxl
in range(self.rs.get_number_of_restraints()):
867 for nstate
in range(self.nstates):
868 prob = state_info[nstate][nxl][
"prob"]
869 low_dist = state_info[nstate][nxl][
"low_dist"]
870 if prob<max_prob_for_violation
or low_dist>min_dist_for_violation:
878 if len(npass)==self.nstates:
880 elif len(nviol)==self.nstates:
884 print(nxl,
'state dists:',[state_info[nstate][nxl][
"low_dist"]
for nstate
in range(self.nstates)],
885 'viol states:',nviol,
'all viol?',all_viol)
886 for nstate
in range(self.nstates):
888 r=0.365; g=0.933; b=0.365;
891 r=0.980; g=0.302; b=0.247;
898 r=0.365; g=0.933; b=0.365;
900 pp = state_info[nstate][nxl][
"low_pp"]
909 cmds[nstate].add((ch1,r1))
910 cmds[nstate].add((ch2,r2))
912 outf = out_fns[nstate]
914 outf.write(
'<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
915 'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,c1[0],c1[1],c1[2],r,g,b))
916 outf.write(
'<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
917 'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv+1,c2[0],c2[1],c2[2],r,g,b))
918 outf.write(
'<link id1= "%d" id2="%d" radius="0.8" '
919 'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,nv+1,r,g,b))
922 for nstate
in range(self.nstates):
923 out_fns[nstate].write(
'</marker_set>\n')
924 out_fns[nstate].close()
926 for ch,r
in cmds[nstate]:
927 cmd+=
'#%i:%i.%s '%(nstate,r,ch)
931 def _get_contribution_info(self,xl,ncontr,use_CA=False):
932 """Return the particles at that contribution. If requested will return CA's instead"""
933 idx1=xl.get_contribution(ncontr)[0]
934 idx2=xl.get_contribution(ncontr)[1]
942 return idx1,idx2,dist
944 def get_best_stats(self,limit_to_state=None,limit_to_chains=None,exclude_chains='',use_CA=False):
945 ''' return the probability, best distance, two coords, and possibly the psi for each xl
946 @param limit_to_state Only examine contributions from one state
947 @param limit_to_chains Returns the particles for certain "easy to visualize" chains
948 @param exclude_chains Even if you limit, don't let one end be in this list.
949 Only works if you also limit chains
950 @param use_CA Limit to CA particles
953 for nxl
in range(self.rs.get_number_of_restraints()):
955 xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
962 for contr
in range(xl.get_number_of_contributions()):
963 pp = xl.get_contribution(contr)
970 if limit_to_state
is not None:
972 if nstate!=limit_to_state:
974 state_contrs.append(contr)
977 if limit_to_chains
is not None:
980 if (c1
in limit_to_chains
or c2
in limit_to_chains)
and (
981 c1
not in exclude_chains
and c2
not in exclude_chains):
982 if dist<low_dist_lim:
989 if limit_to_state
is not None:
990 this_info[
"prob"] = xl.evaluate_for_contributions(state_contrs,
None)
992 this_info[
"prob"] = xl.unprotected_evaluate(
None)
993 if limit_to_chains
is not None:
994 this_info[
"low_pp"] = low_pp_lim
996 this_info[
"low_pp"] = low_pp
998 this_info[
"low_dist"] = low_dist
1001 this_info[
"psi"] = pval
1002 ret.append(this_info)
1005 def print_stats(self):
1008 for nxl,s
in enumerate(stats):
1010 print(s[
"low_dist"])
1014 output = super(AtomicCrossLinkMSRestraint, self).
get_output()
1025 for nxl,s
in enumerate(stats):
1026 if s[
'low_dist']>20.0:
1028 output[
"AtomicXLRestraint_%i_%s"%(nxl,
"Prob")]=str(s[
'prob'])
1029 output[
"AtomicXLRestraint_%i_%s"%(nxl,
"BestDist")]=str(s[
'low_dist'])
1030 if not self.one_psi:
1031 output[
"AtomicXLRestraint_%i_%s"%(nxl,
"psi")]=str(s[
'psi'])
1032 output[
"AtomicXLRestraint_NumViol"] = str(bad_count)
1038 """This restraint allows ambiguous crosslinking between multiple copies
1039 it is a variant of the SimplifiedCrossLinkMS
1049 self.m = representation.prot.get_model()
1055 self.outputlevel =
"low"
1056 self.expdistance = expdistance
1057 self.strength = strength
1059 fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1063 tokens = line.split()
1065 if (tokens[0] ==
"#"):
1074 resolution=resolution,
1076 name_is_ambiguous=
True,
1078 hrc1 = [representation.hier_db.particle_to_name[p]
for p
in ps1]
1080 print(
"ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1085 resolution=resolution,
1087 name_is_ambiguous=
True,
1089 hrc2 = [representation.hier_db.particle_to_name[p]
for p
in ps2]
1091 print(
"ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1096 s1_inds = s1.get_selected_particle_indexes()
1097 s2_inds = s2.get_selected_particle_indexes()
1098 if not set(s1_inds).isdisjoint(set(s2_inds)):
1099 print(
"ConnectivityCrossLinkMS: WARNING> %s %s and %s %s "
1100 "select the same bead(s) - ignoring" % (c1, r1, c2, r2))
1104 if self.strength
is None:
1113 rad1 = rad1 / len(ps1)
1114 rad2 = rad2 / len(ps2)
1116 self.strength = 1 / (rad1 ** 2 + rad2 ** 2)
1124 self.rs.add_restraint(cr)
1125 self.pairs.append((ps1, hrc1, c1, r1, ps2, hrc2, c2, r2, cr))
1138 d1.set_radius(uncertainty1)
1139 d2.set_radius(uncertainty2)
1143 strength = 1 / (uncertainty1 ** 2 + uncertainty2 ** 2)
1150 for i
in range(npoints):
1154 scores.append(cr.unprotected_evaluate(
None))
1155 IMP.pmi.output.plot_xy_data(dists, scores)
1157 def set_label(self, label):
1159 self.rs.set_name(label)
1160 for r
in self.rs.get_restraints():
1163 def add_to_model(self):
1166 def get_restraint_sets(self):
1169 def get_restraint(self):
1172 def set_output_level(self, level="low"):
1174 self.outputlevel = level
1176 def set_weight(self, weight):
1177 self.weight = weight
1178 self.rs.set_weight(weight)
1180 def get_output(self):
1184 score = self.weight * self.rs.unprotected_evaluate(
None)
1185 output[
"_TotalScore"] = str(score)
1186 output[
"ConnectivityCrossLinkMS_Score_" + self.label] = str(score)
1187 for n, p
in enumerate(self.pairs):
1198 for n1, p1
in enumerate(ps1):
1201 for n2, p2
in enumerate(ps2):
1205 label = str(r1) +
":" + name1 +
"_" + str(r2) +
":" + name2
1206 output[
"ConnectivityCrossLinkMS_Distance_" +
1209 label = str(r1) +
":" + c1 +
"_" + str(r2) +
":" + c2
1210 output[
"ConnectivityCrossLinkMS_Score_" +
1211 label] = str(self.weight * cr.unprotected_evaluate(
None))
1215 @
IMP.deprecated_object(
"2.5",
"Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1216 class SimplifiedCrossLinkMS(object):
1226 spheredistancepairscore=
True):
1231 if columnmapping
is None:
1233 columnmapping[
"Protein1"] = 0
1234 columnmapping[
"Protein2"] = 1
1235 columnmapping[
"Residue1"] = 2
1236 columnmapping[
"Residue2"] = 3
1238 self.m = representation.prot.get_model()
1243 self.already_added_pairs = {}
1245 self.outputlevel =
"low"
1246 self.expdistance = expdistance
1247 self.strength = strength
1248 self.truncated = truncated
1249 self.spheredistancepairscore = spheredistancepairscore
1254 protein1 = columnmapping[
"Protein1"]
1255 protein2 = columnmapping[
"Protein2"]
1256 residue1 = columnmapping[
"Residue1"]
1257 residue2 = columnmapping[
"Residue2"]
1259 fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1263 tokens = line.split()
1265 if (tokens[0] ==
"#"):
1267 r1 = int(tokens[residue1])
1268 c1 = tokens[protein1]
1269 r2 = int(tokens[residue2])
1270 c2 = tokens[protein2]
1274 resolution=resolution,
1276 name_is_ambiguous=
False,
1280 resolution=resolution,
1282 name_is_ambiguous=
False,
1287 "residue %d of chain %s selects multiple particles; "
1289 % (r1, c1,
"".join(p.get_name()
for p
in ps1)))
1291 print(
"SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1296 "residue %d of chain %s selects multiple particles; "
1298 % (r2, c2,
"".join(p.get_name()
for p
in ps2)))
1300 print(
"SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1306 if (p1, p2)
in self.already_added_pairs:
1307 dr = self.already_added_pairs[(p1, p2)]
1308 weight = dr.get_weight()
1309 dr.set_weight(weight + 1.0)
1310 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))
1316 limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
1327 if self.spheredistancepairscore:
1332 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2))
1334 self.rs.add_restraint(dr)
1335 self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
1336 self.already_added_pairs[(p1, p2)] = dr
1337 self.already_added_pairs[(p2, p1)] = dr
1339 def set_label(self, label):
1342 def add_to_model(self):
1345 def get_restraint_sets(self):
1348 def get_restraint(self):
1351 def get_restraints(self):
1353 for r
in self.rs.get_restraints():
1354 rlist.append(IMP.core.PairRestraint.get_from(r))
1357 def get_particle_pairs(self):
1359 for i
in range(len(self.pairs)):
1360 p0 = self.pairs[i][0]
1361 p1 = self.pairs[i][1]
1362 ppairs.append((p0, p1))
1365 def set_output_level(self, level="low"):
1367 self.outputlevel = level
1369 def set_weight(self, weight):
1370 self.weight = weight
1371 self.rs.set_weight(weight)
1373 def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
1379 d1.set_radius(radius1)
1380 d2.set_radius(radius2)
1382 limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
1397 for i
in range(npoints):
1401 scores.append(dr.unprotected_evaluate(
None))
1402 IMP.pmi.output.plot_xy_data(dists, scores)
1404 def get_output(self):
1408 score = self.weight * self.rs.unprotected_evaluate(
None)
1409 output[
"_TotalScore"] = str(score)
1410 output[
"SimplifiedCrossLinkMS_Score_" + self.label] = str(score)
1411 for i
in range(len(self.pairs)):
1413 p0 = self.pairs[i][0]
1414 p1 = self.pairs[i][1]
1415 crosslinker =
'standard'
1416 ln = self.pairs[i][2]
1417 resid1 = self.pairs[i][3]
1418 chain1 = self.pairs[i][4]
1419 resid2 = self.pairs[i][5]
1420 chain2 = self.pairs[i][6]
1422 label = str(resid1) +
":" + chain1 +
"_" + \
1423 str(resid2) +
":" + chain2
1424 output[
"SimplifiedCrossLinkMS_Score_" + crosslinker +
"_" +
1425 label] = str(self.weight * ln.unprotected_evaluate(
None))
1427 if self.spheredistancepairscore:
1433 output[
"SimplifiedCrossLinkMS_Distance_" +
1440 @
IMP.deprecated_object(
"2.5",
"Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1441 class SigmoidalCrossLinkMS(object):
1443 self, representation, restraints_file, inflection, slope, amplitude,
1444 linear_slope, resolution=
None, columnmapping=
None, csvfile=
False,
1445 filters=
None, filelabel=
"None"):
1451 if columnmapping
is None:
1453 columnmapping[
"Protein1"] = 0
1454 columnmapping[
"Protein2"] = 1
1455 columnmapping[
"Residue1"] = 2
1456 columnmapping[
"Residue2"] = 3
1462 db = tools.get_db_from_csv(restraints_file)
1464 db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1466 self.m = representation.prot.get_model()
1472 self.already_added_pairs = {}
1473 self.inflection = inflection
1475 self.amplitude = amplitude
1476 self.linear_slope = linear_slope
1478 self.outputlevel =
"low"
1483 protein1 = columnmapping[
"Protein1"]
1484 protein2 = columnmapping[
"Protein2"]
1485 residue1 = columnmapping[
"Residue1"]
1486 residue2 = columnmapping[
"Residue2"]
1488 indb = open(
"included." + filelabel +
".xl.db",
"w")
1489 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
1490 midb = open(
"missing." + filelabel +
".xl.db",
"w")
1494 tokens = entry.split()
1496 if (tokens[0] ==
"#"):
1498 r1 = int(tokens[residue1])
1499 c1 = tokens[protein1]
1500 r2 = int(tokens[residue2])
1501 c2 = tokens[protein2]
1503 if filters
is not None:
1504 if eval(tools.cross_link_db_filter_parser(filters)) ==
False:
1505 exdb.write(str(entry) +
"\n")
1507 indb.write(str(entry) +
"\n")
1508 r1 = int(entry[residue1])
1509 c1 = entry[protein1]
1510 r2 = int(entry[residue2])
1511 c2 = entry[protein2]
1515 resolution=resolution,
1517 name_is_ambiguous=
False,
1521 resolution=resolution,
1523 name_is_ambiguous=
False,
1527 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
1529 print(
"SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1530 midb.write(str(entry) +
"\n")
1534 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
1536 print(
"SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1537 midb.write(str(entry) +
"\n")
1543 if (p1, p2)
in self.already_added_pairs:
1544 dr = self.already_added_pairs[(p1, p2)]
1545 weight = dr.get_weight()
1546 dr.increment_amplitude(amplitude)
1547 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()))
1548 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2)
1549 +
"-ampl:" + str(dr.get_amplitude()))
1562 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2)
1563 +
"-ampl:" + str(dr.get_amplitude()))
1565 self.rs.add_restraint(dr)
1567 self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
1568 self.already_added_pairs[(p1, p2)] = dr
1569 self.already_added_pairs[(p2, p1)] = dr
1571 def set_label(self, label):
1574 def add_to_model(self):
1577 def get_restraint_sets(self):
1580 def get_restraint(self):
1583 def get_restraints(self):
1585 for r
in self.rs.get_restraints():
1586 rlist.append(IMP.core.PairRestraint.get_from(r))
1589 def get_particle_pairs(self):
1591 for i
in range(len(self.pairs)):
1592 p0 = self.pairs[i][0]
1593 p1 = self.pairs[i][1]
1594 ppairs.append((p0, p1))
1597 def set_output_level(self, level="low"):
1599 self.outputlevel = level
1601 def set_weight(self, weight):
1602 self.weight = weight
1603 self.rs.set_weight(weight)
1605 def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
1610 d1.set_radius(radius1)
1611 d2.set_radius(radius2)
1622 for i
in range(npoints):
1626 scores.append(dr.unprotected_evaluate(
None))
1627 IMP.pmi.output.plot_xy_data(dists, scores)
1629 def get_output(self):
1633 score = self.weight * self.rs.unprotected_evaluate(
None)
1634 output[
"_TotalScore"] = str(score)
1635 output[
"SigmoidalCrossLinkMS_Score_" + self.label] = str(score)
1636 for i
in range(len(self.pairs)):
1638 p0 = self.pairs[i][0]
1639 p1 = self.pairs[i][1]
1640 crosslinker =
'standard'
1641 ln = self.pairs[i][2]
1642 resid1 = self.pairs[i][3]
1643 chain1 = self.pairs[i][4]
1644 resid2 = self.pairs[i][5]
1645 chain2 = self.pairs[i][6]
1647 label = str(resid1) +
":" + chain1 +
"_" + \
1648 str(resid2) +
":" + chain2
1649 output[
"SigmoidalCrossLinkMS_Score_" + crosslinker +
"_" +
1650 label +
"_" + self.label] = str(self.weight * ln.unprotected_evaluate(
None))
1653 output[
"SigmoidalCrossLinkMS_Distance_" +
1658 @
IMP.deprecated_object(
"2.5",
"Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1659 class ISDCrossLinkMS(IMP.pmi.restraints._NuisancesBase):
1660 def __init__(self, representation,
1680 automatic_sigma_classification=
False,
1681 attributes_for_label=
None):
1691 if type(representation) != list:
1692 representations = [representation]
1694 representations = representation
1696 if columnmapping
is None:
1698 columnmapping[
"Protein1"] = 0
1699 columnmapping[
"Protein2"] = 1
1700 columnmapping[
"Residue1"] = 2
1701 columnmapping[
"Residue2"] = 3
1702 columnmapping[
"IDScore"] = 4
1703 columnmapping[
"XLUniqueID"] = 5
1709 db = IMP.pmi.tools.get_db_from_csv(restraints_file)
1711 db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1713 indb = open(
"included." + filelabel +
".xl.db",
"w")
1714 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
1715 midb = open(
"missing." + filelabel +
".xl.db",
"w")
1717 self.m = representations[0].prot.get_model()
1727 self.sigma_dictionary = {}
1728 self.psi_dictionary = {}
1729 self.psi_is_sampled =
True
1730 self.sigma_is_sampled =
True
1732 self.dataset = representations[0].get_file_dataset(restraints_file)
1733 if not self.dataset
and os.path.exists(restraints_file):
1734 l = ihm.location.InputFileLocation(restraints_file,
1735 details=
"Crosslinks")
1736 self.dataset = ihm.dataset.CXMSDataset(l)
1738 xl_groups = [p.get_cross_link_group(self)
1739 for p, state
in representations[0]._protocol_output]
1744 self.ids_map = IMP.pmi.tools.map()
1745 self.ids_map.set_map_element(20.0, 0.05)
1746 self.ids_map.set_map_element(65.0, 0.01)
1748 self.ids_map = ids_map
1750 if radius_map
is None:
1751 self.radius_map = IMP.pmi.tools.map()
1752 if automatic_sigma_classification:
1753 self.radius_map.set_map_element(10, 10)
1754 self.radius_map.set_map_element(1, 1)
1756 self.radius_map = radius_map
1758 self.outputlevel =
"low"
1762 self.linear.set_slope(slope)
1765 protein1 = columnmapping[
"Protein1"]
1766 protein2 = columnmapping[
"Protein2"]
1767 residue1 = columnmapping[
"Residue1"]
1768 residue2 = columnmapping[
"Residue2"]
1769 idscore = columnmapping[
"IDScore"]
1771 xluniqueid = columnmapping[
"XLUniqueID"]
1781 uniqueid_restraints_map = {}
1783 for nxl, entry
in enumerate(db):
1785 if not jackknifing
is None:
1792 raise NotImplementedError(
"jackknifing not yet implemented")
1795 tokens = entry.split()
1800 if (tokens[0] ==
"#"):
1803 r1 = int(tokens[residue1])
1804 c1 = tokens[protein1]
1805 r2 = int(tokens[residue2])
1806 c2 = tokens[protein2]
1808 if offset_dict
is not None:
1809 if c1
in offset_dict: r1+=offset_dict[c1]
1810 if c2
in offset_dict: r2+=offset_dict[c2]
1812 if rename_dict
is not None:
1813 if c1
in rename_dict: c1=rename_dict[c1]
1814 if c2
in rename_dict: c2=rename_dict[c2]
1819 ids = float(tokens[idscore])
1820 if xluniqueid
is None:
1823 xlid = tokens[xluniqueid]
1825 print(
"this line was not accessible " + str(entry))
1826 if residue1
not in entry: print(str(residue1)+
" keyword not in database")
1827 if residue2
not in entry: print(str(residue2)+
" keyword not in database")
1828 if protein1
not in entry: print(str(protein1)+
" keyword not in database")
1829 if protein2
not in entry: print(str(protein2)+
" keyword not in database")
1830 if idscore
not in entry: print(str(idscore)+
" keyword not in database")
1831 if xluniqueid
not in entry: print(str(xluniqueid)+
" keyword not in database")
1835 if filters
is not None:
1837 exdb.write(str(entry) +
"\n")
1841 r1 = int(entry[residue1])
1842 c1 = entry[protein1]
1843 r2 = int(entry[residue2])
1844 c2 = entry[protein2]
1846 if offset_dict
is not None:
1847 if c1
in offset_dict: r1+=offset_dict[c1]
1848 if c2
in offset_dict: r2+=offset_dict[c2]
1850 if rename_dict
is not None:
1851 if c1
in rename_dict: c1=rename_dict[c1]
1852 if c2
in rename_dict: c2=rename_dict[c2]
1858 ids = float(entry[idscore])
1860 ids = entry[idscore]
1861 if xluniqueid
is None:
1864 xlid = entry[xluniqueid]
1867 print(
"this line was not accessible " + str(entry))
1868 if residue1
not in entry: print(str(residue1)+
" keyword not in database")
1869 if residue2
not in entry: print(str(residue2)+
" keyword not in database")
1870 if protein1
not in entry: print(str(protein1)+
" keyword not in database")
1871 if protein2
not in entry: print(str(protein2)+
" keyword not in database")
1872 if idscore
not in entry: print(str(idscore)+
" keyword not in database")
1873 if xluniqueid
not in entry: print(str(xluniqueid)+
" keyword not in database")
1877 ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
1879 for p, group
in zip(representations[0]._protocol_output,
1882 for nstate, r
in enumerate(representations):
1887 resolution=resolution,
1889 name_is_ambiguous=
False,
1893 resolution=resolution,
1895 name_is_ambiguous=
False,
1899 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
1901 print(
"ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1902 midb.write(str(entry) +
"\n")
1906 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
1908 print(
"ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1909 midb.write(str(entry) +
"\n")
1915 if (p1 == p2)
and (r1 == r2) :
1916 print(
"ISDCrossLinkMS Restraint: WARNING> on the identical bead particles and the identical residues, thus skipping this cross-link.")
1919 if xlid
in uniqueid_restraints_map:
1920 print(
"Appending a crosslink restraint into the uniqueID %s" % str(xlid))
1921 dr = uniqueid_restraints_map[xlid]
1923 print(
"Generating a NEW crosslink restraint with a uniqueID %s" % str(xlid))
1928 restraints.append(dr)
1929 uniqueid_restraints_map[xlid] = dr
1931 mappedr1 = self.radius_map.get_map_element(
1933 sigma1 = self.get_sigma(mappedr1)[0]
1934 mappedr2 = self.radius_map.get_map_element(
1936 sigma2 = self.get_sigma(mappedr2)[0]
1938 psival = self.ids_map.get_map_element(ids)
1939 psi = self.get_psi(psival)[0]
1942 p1i = p1.get_particle_index()
1943 p2i = p2.get_particle_index()
1951 dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
1952 print(
"--------------")
1953 print(
"ISDCrossLinkMS: generating cross-link restraint between")
1954 print(
"ISDCrossLinkMS: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
1955 print(
"ISDCrossLinkMS: with sigma1 %f sigma2 %f psi %s" % (mappedr1, mappedr2, psival))
1956 print(
"ISDCrossLinkMS: between particles %s and %s" % (p1.get_name(), p2.get_name()))
1957 print(
"==========================================\n")
1958 indb.write(str(entry) +
"\n")
1959 for p, ex_xl
in zip(representations[0]._protocol_output,
1961 p[0].add_cross_link(p[1], ex_xl[0], p1, p2, length,
1962 sigma1, sigma2, psi, ex_xl[1])
1969 xlattribute =
"intrarb"
1971 xlattribute =
"interrb"
1974 if not attributes_for_label
is None:
1975 for a
in attributes_for_label:
1976 xlattribute = xlattribute +
"_" + str(entry[a])
1978 xlattribute = xlattribute +
"-State:" + str(nstate)
1981 xlattribute +
"-" + c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2) +
"_" + self.label)
1986 xlattribute +
"-" + c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2) +
"_" + self.label)
1987 self.rslin.add_restraint(pr)
2007 self.rs.add_restraint(lw)
2010 model = property(
lambda s: s.m)
2012 def set_weight(self, weight):
2013 self.weight = weight
2014 self.rs.set_weight(weight)
2016 def set_slope_linear_term(self, slope):
2017 self.linear.set_slope(slope)
2019 def set_label(self, label):
2022 def add_to_model(self):
2029 def get_hierarchies(self):
2032 def get_restraint_sets(self):
2035 def get_restraint(self):
2038 def get_restraint_for_rmf(self):
2041 def get_restraints(self):
2043 for r
in self.rs.get_restraints():
2044 rlist.append(IMP.core.PairRestraint.get_from(r))
2047 def get_particle_pairs(self):
2049 for i
in range(len(self.pairs)):
2050 p0 = self.pairs[i][0]
2051 p1 = self.pairs[i][1]
2052 ppairs.append((p0, p1))
2055 def set_output_level(self, level="low"):
2057 self.outputlevel = level
2059 def set_psi_is_sampled(self, is_sampled=True):
2060 self.psi_is_sampled = is_sampled
2062 def set_sigma_is_sampled(self, is_sampled=True):
2063 self.sigma_is_sampled = is_sampled
2065 def get_label(self,pairs_index):
2066 resid1 = self.pairs[pairs_index][3]
2067 chain1 = self.pairs[pairs_index][4]
2068 resid2 = self.pairs[pairs_index][5]
2069 chain2 = self.pairs[pairs_index][6]
2070 attribute = self.pairs[pairs_index][7]
2071 rad1 = self.pairs[pairs_index][8]
2072 rad2 = self.pairs[pairs_index][9]
2073 psi = self.pairs[pairs_index][10]
2074 xlid= self.pairs[pairs_index][11]
2075 label = attribute +
"-" + \
2076 str(resid1) +
":" + chain1 +
"_" + str(resid2) +
":" + \
2077 chain2 +
"-" + str(rad1) +
"-" + str(rad2) +
"-" + str(psi)
2080 def write_db(self,filename):
2081 cldb=IMP.pmi.output.CrossLinkIdentifierDatabase()
2083 for pairs_index
in range(len(self.pairs)):
2085 resid1 = self.pairs[pairs_index][3]
2086 chain1 = self.pairs[pairs_index][4]
2087 resid2 = self.pairs[pairs_index][5]
2088 chain2 = self.pairs[pairs_index][6]
2089 attribute = self.pairs[pairs_index][7]
2090 rad1 = self.pairs[pairs_index][8]
2091 rad2 = self.pairs[pairs_index][9]
2092 psi = self.pairs[pairs_index][10]
2093 xlid= self.pairs[pairs_index][11]
2094 nstate=self.pairs[pairs_index][12]
2095 ids=self.pairs[pairs_index][13]
2097 label=self.get_label(pairs_index)
2098 cldb.set_unique_id(label,xlid)
2099 cldb.set_protein1(label,chain1)
2100 cldb.set_protein2(label,chain2)
2101 cldb.set_residue1(label,resid1)
2102 cldb.set_residue2(label,resid2)
2103 cldb.set_idscore(label,ids)
2104 cldb.set_state(label,nstate)
2105 cldb.set_sigma1(label,rad1)
2106 cldb.set_sigma2(label,rad2)
2107 cldb.set_psi(label,psi)
2108 cldb.write(filename)
2111 def get_output(self):
2115 score = self.weight * self.rs.unprotected_evaluate(
None)
2116 output[
"_TotalScore"] = str(score)
2117 output[
"ISDCrossLinkMS_Data_Score_" + self.label] = str(score)
2118 output[
"ISDCrossLinkMS_PriorSig_Score_" +
2119 self.label] = self.rssig.unprotected_evaluate(
None)
2120 output[
"ISDCrossLinkMS_PriorPsi_Score_" +
2121 self.label] = self.rspsi.unprotected_evaluate(
None)
2122 output[
"ISDCrossLinkMS_Linear_Score_" +
2123 self.label] = self.rslin.unprotected_evaluate(
None)
2124 for i
in range(len(self.pairs)):
2126 label=self.get_label(i)
2127 ln = self.pairs[i][2]
2128 p0 = self.pairs[i][0]
2129 p1 = self.pairs[i][1]
2130 output[
"ISDCrossLinkMS_Score_" +
2131 label +
"_" + self.label] = str(self.weight * -log(ln.unprotected_evaluate(
None)))
2135 output[
"ISDCrossLinkMS_Distance_" +
2139 for psiindex
in self.psi_dictionary:
2140 output[
"ISDCrossLinkMS_Psi_" +
2141 str(psiindex) +
"_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2143 for resolution
in self.sigma_dictionary:
2144 output[
"ISDCrossLinkMS_Sigma_" +
2145 str(resolution) +
"_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
2150 def get_particles_to_sample(self):
2152 for resolution
in self.sigma_dictionary:
2153 if self.sigma_dictionary[resolution][2]
and self.sigma_is_sampled:
2154 ps[
"Nuisances_ISDCrossLinkMS_Sigma_" + str(resolution) +
"_" + self.label] =\
2155 ([self.sigma_dictionary[resolution][0]],
2156 self.sigma_dictionary[resolution][1])
2157 if self.psi_is_sampled:
2158 for psiindex
in self.psi_dictionary:
2159 if self.psi_dictionary[psiindex][2]:
2160 ps[
"Nuisances_ISDCrossLinkMS_Psi_" +
2161 str(psiindex) +
"_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2165 class CysteineCrossLinkRestraint(object):
2166 def __init__(self, representations, filename, cbeta=False,
2167 betatuple=(0.03, 0.1),
2168 disttuple=(0.0, 25.0, 1000),
2169 omegatuple=(1.0, 1000.0, 50),
2170 sigmatuple=(0.3, 0.3, 1),
2172 sigmaissampled=
False,
2173 weightissampled=
True,
2174 epsilonissampled=
True
2181 self.representations = representations
2182 self.m = self.representations[0].prot.get_model()
2185 self.epsilonmaxtrans = 0.01
2186 self.sigmamaxtrans = 0.1
2187 self.betamaxtrans = 0.01
2188 self.weightmaxtrans = 0.1
2190 self.outputlevel =
"low"
2191 self.betaissampled = betaissampled
2192 self.sigmaissampled = sigmaissampled
2193 self.weightissampled = weightissampled
2194 self.epsilonissampled = epsilonissampled
2196 betalower = betatuple[0]
2197 betaupper = betatuple[1]
2204 self.beta = IMP.pmi.tools.SetupNuisance(
2209 betaissampled).get_particle(
2212 self.sigma = IMP.pmi.tools.SetupNuisance(
2217 sigmaissampled).get_particle()
2219 self.weight = IMP.pmi.tools.SetupWeight(
2221 weightissampled).get_particle(
2225 fl = IMP.pmi.tools.open_file_or_inline_text(filename)
2234 if t[5]
in self.epsilons:
2235 if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
2236 self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
2238 self.epsilons[t[5]] = IMP.pmi.tools.SetupNuisance(self.m,
2239 0.01, 0.01, 1.0 - float(t[4]), epsilonissampled).get_particle()
2240 up = self.epsilons[t[5]].get_upper()
2241 low = self.epsilons[t[5]].get_lower()
2243 self.epsilons[t[5]].set_upper(low)
2245 data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
2249 crossdata = IMP.pmi.tools.get_cross_link_data(
2250 "cysteine",
"cysteine_CA_FES.txt.standard",
2251 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2253 crossdata = IMP.pmi.tools.get_cross_link_data(
2254 "cysteine",
"cysteine_CB_FES.txt.standard",
2255 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2258 fmod_grid = IMP.pmi.tools.get_grid(0.0, 1.0, 300,
True)
2259 omega2_grid = IMP.pmi.tools.get_log_grid(0.001, 10000.0, 100)
2260 beta_grid = IMP.pmi.tools.get_log_grid(betalower, betaupper, betangrid)
2263 print(
"--------------")
2264 print(
"CysteineCrossLink: attempting to create a restraint " + str(d))
2284 self.epsilons[epslabel],
2290 for i, representation
in enumerate(self.representations):
2297 resolution=1, name=chain1,
2298 name_is_ambiguous=
False, residue=resid1)[0]
2304 resolution=1, name=chain2,
2305 name_is_ambiguous=
False, residue=resid2)[0]
2314 for t
in range(-1, 2):
2316 resolution=1, name=chain1,
2317 name_is_ambiguous=
False, residue=resid1 + t)
2323 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
2326 resolution=1, name=chain2,
2327 name_is_ambiguous=
False, residue=resid2 + t)
2333 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
2336 if (p1
is not None and p2
is not None):
2337 ccl.add_contribution(p1, p2)
2341 print(
"Distance_" + str(resid1) +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2,
IMP.core.get_distance(d1, d2))
2344 if (len(p1) == 3
and len(p2) == 3):
2345 p11n = p1[0].get_name()
2346 p12n = p1[1].get_name()
2347 p13n = p1[2].get_name()
2348 p21n = p2[0].get_name()
2349 p22n = p2[1].get_name()
2350 p23n = p2[2].get_name()
2352 print(
"CysteineCrossLink: generating CB cysteine cross-link restraint between")
2353 print(
"CysteineCrossLink: residue %d of chain %s and residue %d of chain %s" % (resid1, chain1, resid2, chain2))
2354 print(
"CysteineCrossLink: between particles %s %s %s and %s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
2356 ccl.add_contribution(p1, p2)
2359 self.rs.add_restraint(ccl)
2360 ccl.set_name(
"CysteineCrossLink_" + str(resid1)
2361 +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2)
2363 def set_label(self, label):
2366 def add_to_model(self):
2369 def get_particles_to_sample(self):
2371 if self.epsilonissampled:
2372 for eps
in self.epsilons.keys():
2373 ps[
"Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
2374 self.label] = ([self.epsilons[eps]], self.epsilonmaxtrans)
2375 if self.betaissampled:
2376 ps[
"Nuisances_CysteineCrossLinkRestraint_beta_" +
2377 self.label] = ([self.beta], self.betamaxtrans)
2378 if self.weightissampled:
2379 ps[
"Weights_CysteineCrossLinkRestraint_" +
2380 self.label] = ([self.weight], self.weightmaxtrans)
2381 if self.sigmaissampled:
2382 ps[
"Nuisances_CysteineCrossLinkRestraint_" +
2383 self.label] = ([self.sigma], self.sigmamaxtrans)
2386 def set_output_level(self, level="low"):
2388 self.outputlevel = level
2390 def get_restraint(self):
2393 def get_sigma(self):
2396 def get_output(self):
2398 score = self.rs.unprotected_evaluate(
None)
2399 output[
"_TotalScore"] = str(score)
2400 output[
"CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
2401 output[
"CysteineCrossLinkRestraint_sigma_" +
2402 self.label] = str(self.sigma.get_scale())
2403 for eps
in self.epsilons.keys():
2404 output[
"CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
2405 self.label] = str(self.epsilons[eps].get_scale())
2406 output[
"CysteineCrossLinkRestraint_beta_" +
2407 self.label] = str(self.beta.get_scale())
2408 for n
in range(self.weight.get_number_of_states()):
2409 output[
"CysteineCrossLinkRestraint_weights_" +
2410 str(n) +
"_" + self.label] = str(self.weight.get_weight(n))
2412 if self.outputlevel ==
"high":
2413 for rst
in self.rs.get_restraints():
2414 output[
"CysteineCrossLinkRestraint_Total_Frequency_" +
2415 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2416 "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_model_frequency()
2417 output[
"CysteineCrossLinkRestraint_Standard_Error_" +
2418 IMP.isd.CysteineCrossLinkRestraint.get_from(
2421 + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_standard_error()
2422 if len(self.representations) > 1:
2423 for i
in range(len(self.prots)):
2424 output[
"CysteineCrossLinkRestraint_Frequency_Contribution_" +
2425 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2426 "_State_" + str(i) +
"_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_frequencies()[i]
def set_sigma_is_sampled
Switch on/off the sampling of sigma particles.
Setup cross-link distance restraints from mass spectrometry data.
A function that is harmonic over an interval.
A restraint for ambiguous cross-linking MS data and multiple state approach.
def get_best_stats
return the probability, best distance, two coords, and possibly the psi for each xl ...
def get_hierarchies
get the hierarchy
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
static XYZR setup_particle(Model *m, ParticleIndex pi)
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi.sampler object.
Restrain atom pairs based on a set of crosslinks.
Handles cross-link data sets.
Add scale parameter to particle.
Add uncertainty to a particle.
A score on the distance between the surfaces of two spheres.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi.sampler object.
def get_restraint_for_rmf
Get the restraint for visualization in an RMF file.
Object used to hold a set of restraints.
Class for storing model, its restraints, constraints, and particles.
Classes to handle different kinds of restraints.
def set_output_level
Set the output level of the output.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def get_restraints
get the restraints in a list
A decorator for a particle representing an atom.
This class allows to create filter functions that can be passed to the CrossLinkDataBase in this way:...
def get_movers
Get all need data to construct a mover in IMP.pmi.dof class.
def get_particle_pairs
Get a list of tuples containing the particle pairs.
def create_psi
This is called internally.
Simple sigmoidal score calculated between sphere surfaces.
A decorator for a particle with x,y,z coordinates.
def create_sigma
This is called internally.
This restraint allows ambiguous crosslinking between multiple copies it is a variant of the Simplifie...
int get_state_index(Hierarchy h)
Walk up the hierarchy to find the current state.
def plot_violations
Create CMM files, one for each state, of all crosslinks.
Modify a set of continuous variables using a normal distribution.
Classes for writing output files and processing them.
def deprecated_object
Python decorator to mark a class as deprecated.
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...
def create_restraints_for_rmf
create dummy harmonic restraints for each XL but don't add to model Makes it easy to see each contrib...
The general base class for IMP exceptions.
Setup cross-link distance restraints at atomic level The "atomic" aspect is that it models the part...
Residue get_residue(Atom d, bool nothrow=false)
Return the Residue containing this atom.
Calculate the -Log of a list of restraints.
Class to handle individual particles of a Model object.
double get_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the distance between two vectors.
def get_restraint_sets
get the restraint set
Restraint * create_connectivity_restraint(const Selections &s, double x0, double k, std::string name="Connectivity%1%")
Create a restraint connecting the selections.
Applies a PairScore to a Pair.
def get_output
Get the output of the restraint to be used by the IMP.pmi.output object.
def get_restraint_for_rmf
get the dummy restraints to be displayed in the rmf file
Functionality for loading, creating, manipulating and scoring atomic structures.
std::string get_chain_id(Hierarchy h)
Walk up the hierarchy to determine the chain id.
def load_nuisances_from_stat_file
Read a stat file and load all the sigmas.
Select hierarchy particles identified by the biological name.
def get_output
Get outputs to write to stat files.
def set_psi_is_sampled
Switch on/off the sampling of psi particles.
Base class for PMI restraints, which wrap IMP.Restraint(s).
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.