1 """@namespace IMP.pmi1.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.pmi1.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.pmi1.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.pmi1.io.crosslink.CrossLinkDataBase):
80 raise TypeError(
"CrossLinkingMassSpectrometryRestraint: CrossLinkDataBase should be an IMP.pmi1.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"
110 xl_groups = [p.get_cross_link_group(self)
111 for p, state
in representations[0]._protocol_output]
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]
156 ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
159 zip(representations[0]._protocol_output,
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")
284 for p, ex_xl
in zip(representations[0]._protocol_output,
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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1659 class ISDCrossLinkMS(IMP.pmi1.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.pmi1.tools.get_db_from_csv(restraints_file)
1711 db = IMP.pmi1.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.pmi1.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.pmi1.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)
2009 def set_weight(self, weight):
2010 self.weight = weight
2011 self.rs.set_weight(weight)
2013 def set_slope_linear_term(self, slope):
2014 self.linear.set_slope(slope)
2016 def set_label(self, label):
2019 def add_to_model(self):
2026 def get_hierarchies(self):
2029 def get_restraint_sets(self):
2032 def get_restraint(self):
2035 def get_restraint_for_rmf(self):
2038 def get_restraints(self):
2040 for r
in self.rs.get_restraints():
2041 rlist.append(IMP.core.PairRestraint.get_from(r))
2044 def get_particle_pairs(self):
2046 for i
in range(len(self.pairs)):
2047 p0 = self.pairs[i][0]
2048 p1 = self.pairs[i][1]
2049 ppairs.append((p0, p1))
2052 def set_output_level(self, level="low"):
2054 self.outputlevel = level
2056 def set_psi_is_sampled(self, is_sampled=True):
2057 self.psi_is_sampled = is_sampled
2059 def set_sigma_is_sampled(self, is_sampled=True):
2060 self.sigma_is_sampled = is_sampled
2062 def get_label(self,pairs_index):
2063 resid1 = self.pairs[pairs_index][3]
2064 chain1 = self.pairs[pairs_index][4]
2065 resid2 = self.pairs[pairs_index][5]
2066 chain2 = self.pairs[pairs_index][6]
2067 attribute = self.pairs[pairs_index][7]
2068 rad1 = self.pairs[pairs_index][8]
2069 rad2 = self.pairs[pairs_index][9]
2070 psi = self.pairs[pairs_index][10]
2071 xlid= self.pairs[pairs_index][11]
2072 label = attribute +
"-" + \
2073 str(resid1) +
":" + chain1 +
"_" + str(resid2) +
":" + \
2074 chain2 +
"-" + str(rad1) +
"-" + str(rad2) +
"-" + str(psi)
2077 def write_db(self,filename):
2078 cldb=IMP.pmi1.output.CrossLinkIdentifierDatabase()
2080 for pairs_index
in range(len(self.pairs)):
2082 resid1 = self.pairs[pairs_index][3]
2083 chain1 = self.pairs[pairs_index][4]
2084 resid2 = self.pairs[pairs_index][5]
2085 chain2 = self.pairs[pairs_index][6]
2086 attribute = self.pairs[pairs_index][7]
2087 rad1 = self.pairs[pairs_index][8]
2088 rad2 = self.pairs[pairs_index][9]
2089 psi = self.pairs[pairs_index][10]
2090 xlid= self.pairs[pairs_index][11]
2091 nstate=self.pairs[pairs_index][12]
2092 ids=self.pairs[pairs_index][13]
2094 label=self.get_label(pairs_index)
2095 cldb.set_unique_id(label,xlid)
2096 cldb.set_protein1(label,chain1)
2097 cldb.set_protein2(label,chain2)
2098 cldb.set_residue1(label,resid1)
2099 cldb.set_residue2(label,resid2)
2100 cldb.set_idscore(label,ids)
2101 cldb.set_state(label,nstate)
2102 cldb.set_sigma1(label,rad1)
2103 cldb.set_sigma2(label,rad2)
2104 cldb.set_psi(label,psi)
2105 cldb.write(filename)
2108 def get_output(self):
2112 score = self.weight * self.rs.unprotected_evaluate(
None)
2113 output[
"_TotalScore"] = str(score)
2114 output[
"ISDCrossLinkMS_Data_Score_" + self.label] = str(score)
2115 output[
"ISDCrossLinkMS_PriorSig_Score_" +
2116 self.label] = self.rssig.unprotected_evaluate(
None)
2117 output[
"ISDCrossLinkMS_PriorPsi_Score_" +
2118 self.label] = self.rspsi.unprotected_evaluate(
None)
2119 output[
"ISDCrossLinkMS_Linear_Score_" +
2120 self.label] = self.rslin.unprotected_evaluate(
None)
2121 for i
in range(len(self.pairs)):
2123 label=self.get_label(i)
2124 ln = self.pairs[i][2]
2125 p0 = self.pairs[i][0]
2126 p1 = self.pairs[i][1]
2127 output[
"ISDCrossLinkMS_Score_" +
2128 label +
"_" + self.label] = str(self.weight * -log(ln.unprotected_evaluate(
None)))
2132 output[
"ISDCrossLinkMS_Distance_" +
2136 for psiindex
in self.psi_dictionary:
2137 output[
"ISDCrossLinkMS_Psi_" +
2138 str(psiindex) +
"_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2140 for resolution
in self.sigma_dictionary:
2141 output[
"ISDCrossLinkMS_Sigma_" +
2142 str(resolution) +
"_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
2147 def get_particles_to_sample(self):
2149 for resolution
in self.sigma_dictionary:
2150 if self.sigma_dictionary[resolution][2]
and self.sigma_is_sampled:
2151 ps[
"Nuisances_ISDCrossLinkMS_Sigma_" + str(resolution) +
"_" + self.label] =\
2152 ([self.sigma_dictionary[resolution][0]],
2153 self.sigma_dictionary[resolution][1])
2154 if self.psi_is_sampled:
2155 for psiindex
in self.psi_dictionary:
2156 if self.psi_dictionary[psiindex][2]:
2157 ps[
"Nuisances_ISDCrossLinkMS_Psi_" +
2158 str(psiindex) +
"_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2162 class CysteineCrossLinkRestraint(object):
2163 def __init__(self, representations, filename, cbeta=False,
2164 betatuple=(0.03, 0.1),
2165 disttuple=(0.0, 25.0, 1000),
2166 omegatuple=(1.0, 1000.0, 50),
2167 sigmatuple=(0.3, 0.3, 1),
2169 sigmaissampled=
False,
2170 weightissampled=
True,
2171 epsilonissampled=
True
2178 self.representations = representations
2179 self.m = self.representations[0].prot.get_model()
2182 self.epsilonmaxtrans = 0.01
2183 self.sigmamaxtrans = 0.1
2184 self.betamaxtrans = 0.01
2185 self.weightmaxtrans = 0.1
2187 self.outputlevel =
"low"
2188 self.betaissampled = betaissampled
2189 self.sigmaissampled = sigmaissampled
2190 self.weightissampled = weightissampled
2191 self.epsilonissampled = epsilonissampled
2193 betalower = betatuple[0]
2194 betaupper = betatuple[1]
2201 self.beta = IMP.pmi1.tools.SetupNuisance(
2206 betaissampled).get_particle(
2209 self.sigma = IMP.pmi1.tools.SetupNuisance(
2214 sigmaissampled).get_particle()
2216 self.weight = IMP.pmi1.tools.SetupWeight(
2218 isoptimized=
False).get_particle(
2222 fl = IMP.pmi1.tools.open_file_or_inline_text(filename)
2231 if t[5]
in self.epsilons:
2232 if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
2233 self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
2235 self.epsilons[t[5]] = IMP.pmi1.tools.SetupNuisance(self.m,
2236 0.01, 0.01, 1.0 - float(t[4]), epsilonissampled).get_particle()
2237 up = self.epsilons[t[5]].get_upper()
2238 low = self.epsilons[t[5]].get_lower()
2240 self.epsilons[t[5]].set_upper(low)
2242 data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
2246 crossdata = IMP.pmi1.tools.get_cross_link_data(
2247 "cysteine",
"cysteine_CA_FES.txt.standard",
2248 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2250 crossdata = IMP.pmi1.tools.get_cross_link_data(
2251 "cysteine",
"cysteine_CB_FES.txt.standard",
2252 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2255 fmod_grid = IMP.pmi1.tools.get_grid(0.0, 1.0, 300,
True)
2256 omega2_grid = IMP.pmi1.tools.get_log_grid(0.001, 10000.0, 100)
2257 beta_grid = IMP.pmi1.tools.get_log_grid(betalower, betaupper, betangrid)
2260 print(
"--------------")
2261 print(
"CysteineCrossLink: attempting to create a restraint " + str(d))
2281 self.epsilons[epslabel],
2287 for i, representation
in enumerate(self.representations):
2294 resolution=1, name=chain1,
2295 name_is_ambiguous=
False, residue=resid1)[0]
2301 resolution=1, name=chain2,
2302 name_is_ambiguous=
False, residue=resid2)[0]
2311 for t
in range(-1, 2):
2313 resolution=1, name=chain1,
2314 name_is_ambiguous=
False, residue=resid1 + t)
2320 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
2323 resolution=1, name=chain2,
2324 name_is_ambiguous=
False, residue=resid2 + t)
2330 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
2333 if (p1
is not None and p2
is not None):
2334 ccl.add_contribution(p1, p2)
2338 print(
"Distance_" + str(resid1) +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2,
IMP.core.get_distance(d1, d2))
2341 if (len(p1) == 3
and len(p2) == 3):
2342 p11n = p1[0].get_name()
2343 p12n = p1[1].get_name()
2344 p13n = p1[2].get_name()
2345 p21n = p2[0].get_name()
2346 p22n = p2[1].get_name()
2347 p23n = p2[2].get_name()
2349 print(
"CysteineCrossLink: generating CB cysteine cross-link restraint between")
2350 print(
"CysteineCrossLink: residue %d of chain %s and residue %d of chain %s" % (resid1, chain1, resid2, chain2))
2351 print(
"CysteineCrossLink: between particles %s %s %s and %s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
2353 ccl.add_contribution(p1, p2)
2356 self.rs.add_restraint(ccl)
2357 ccl.set_name(
"CysteineCrossLink_" + str(resid1)
2358 +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2)
2361 self.weight.get_particle()
2362 ).set_weights_are_optimized(weightissampled)
2364 def set_label(self, label):
2367 def add_to_model(self):
2370 def get_particles_to_sample(self):
2372 if self.epsilonissampled:
2373 for eps
in self.epsilons.keys():
2374 ps[
"Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
2375 self.label] = ([self.epsilons[eps]], self.epsilonmaxtrans)
2376 if self.betaissampled:
2377 ps[
"Nuisances_CysteineCrossLinkRestraint_beta_" +
2378 self.label] = ([self.beta], self.betamaxtrans)
2379 if self.weightissampled:
2380 ps[
"Weights_CysteineCrossLinkRestraint_" +
2381 self.label] = ([self.weight], self.weightmaxtrans)
2382 if self.sigmaissampled:
2383 ps[
"Nuisances_CysteineCrossLinkRestraint_" +
2384 self.label] = ([self.sigma], self.sigmamaxtrans)
2387 def set_output_level(self, level="low"):
2389 self.outputlevel = level
2391 def get_restraint(self):
2394 def get_sigma(self):
2397 def get_output(self):
2399 score = self.rs.unprotected_evaluate(
None)
2400 output[
"_TotalScore"] = str(score)
2401 output[
"CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
2402 output[
"CysteineCrossLinkRestraint_sigma_" +
2403 self.label] = str(self.sigma.get_scale())
2404 for eps
in self.epsilons.keys():
2405 output[
"CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
2406 self.label] = str(self.epsilons[eps].get_scale())
2407 output[
"CysteineCrossLinkRestraint_beta_" +
2408 self.label] = str(self.beta.get_scale())
2409 for n
in range(self.weight.get_number_of_weights()):
2410 output[
"CysteineCrossLinkRestraint_weights_" +
2411 str(n) +
"_" + self.label] = str(self.weight.get_weight(n))
2413 if self.outputlevel ==
"high":
2414 for rst
in self.rs.get_restraints():
2415 output[
"CysteineCrossLinkRestraint_Total_Frequency_" +
2416 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2417 "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_model_frequency()
2418 output[
"CysteineCrossLinkRestraint_Standard_Error_" +
2419 IMP.isd.CysteineCrossLinkRestraint.get_from(
2422 + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_standard_error()
2423 if len(self.representations) > 1:
2424 for i
in range(len(self.prots)):
2425 output[
"CysteineCrossLinkRestraint_Frequency_Contribution_" +
2426 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2427 "_State_" + str(i) +
"_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_frequencies()[i]
Setup cross-link distance restraints from mass spectrometry data.
Add weights to a particle.
A function that is harmonic over an bounded interval.
A restraint for ambiguous cross-linking MS data and multiple state approach.
def get_restraints
get the restraints in a list
def get_restraint_sets
get the restraint set
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
Handles cross-link data sets.
static XYZR setup_particle(Model *m, ParticleIndex pi)
Restrain atom pairs based on a set of crosslinks.
def get_best_stats
return the probability, best distance, two coords, and possibly the psi for each xl ...
def set_output_level
Set the output level of the output.
Classes to handle different kinds of restraints.
def get_hierarchies
get the hierarchy
Add scale parameter to particle.
A score on the distance between the surfaces of two spheres.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Score distance between two particle centers using a harmonic function.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
def get_restraint_for_rmf
get the dummy restraints to be displayed in the rmf file
Object used to hold a set of restraints.
Class for storing model, its restraints, constraints, and particles.
Setup cross-link distance restraints at atomic level The "atomic" aspect is that it models the part...
def get_output
Get outputs to write to stat files.
Add uncertainty to a particle.
def set_psi_is_sampled
Switch on/off the sampling of psi particles.
def load_nuisances_from_stat_file
Read a stat file and load all the sigmas.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def create_sigma
This is called internally.
Simple sigmoidal score calculated between sphere surfaces.
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_particles_to_sample
Get the particles to be sampled by the IMP.pmi1.sampler object.
def set_sigma_is_sampled
Switch on/off the sampling of sigma particles.
A decorator for a particle with x,y,z coordinates.
int get_state_index(Hierarchy h)
Walk up the hierarchy to find the current state.
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi1.sampler object.
Base class for PMI restraints, which wrap IMP.Restraint(s).
Modify a set of continuous variables using a normal distribution.
Score a pair of particles based on the distance between their centers.
def deprecated_object
Python decorator to mark a class as deprecated.
def get_movers
Get all need data to construct a mover in IMP.pmi1.dof class.
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 plot_violations
Create CMM files, one for each state, of all crosslinks.
def get_restraint_for_rmf
Get the restraint for visualization in an RMF file.
The general base class for IMP exceptions.
Residue get_residue(Atom d, bool nothrow=false)
Return the Residue containing this atom.
Calculate the -Log of a list of restraints.
This restraint allows ambiguous crosslinking between multiple copies it is a variant of the Simplifie...
def get_particle_pairs
Get a list of tuples containing the particle pairs.
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_output
Get the output of the restraint to be used by the IMP.pmi1.output object.
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.
Classes for writing output files and processing them.
def create_psi
This is called internally.
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.
Select hierarchy particles identified by the biological name.
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...
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.