1 """@namespace IMP.pmi.restraints.crosslinking
2 Restraints for handling crosslinking data.
5 from __future__
import print_function
16 from collections
import defaultdict
21 """Setup cross-link distance restraints from mass spectrometry data.
22 The noise in the data and the structural uncertainty of cross-linked amino-acids
23 is inferred using Bayes theory of probability
24 \note Wraps an IMP::isd::CrossLinkMSRestraint
26 def __init__(self, representation=None,
28 CrossLinkDataBase=
None,
34 attributes_for_label=
None):
36 @param representation DEPRECATED The IMP.pmi.representation.Representation
37 object that contain the molecular system
38 @param root_hier The canonical hierarchy containing all the states
39 @param CrossLinkDataBase The IMP.pmi.io.crosslink.CrossLinkDataBase
40 object that contains the cross-link dataset
41 @param length maximal cross-linker length (including the residue sidechains)
42 @param resolution what representation resolution should the cross-link
43 restraint be applied to.
44 @param slope The slope of a distance-linear scoring function that
45 funnels the score when the particles are
46 too far away. Suggested value 0.02.
47 @param label the extra text to label the restraint so that it is
48 searchable in the output
49 @param filelabel automatically generated file containing missing/included/excluded
50 cross-links will be labeled using this text
51 @param attributes_for_label
55 if representation
is not None:
57 if type(representation) != list:
58 representations = [representation]
60 representations = representation
61 self.m = representations[0].prot.get_model()
62 elif root_hier
is not None:
63 self.m = root_hier.get_model()
65 raise Exception(
"You must pass either representation or root_hier")
67 if CrossLinkDataBase
is None:
68 raise Exception(
"You must pass a CrossLinkDataBase")
69 if not isinstance(CrossLinkDataBase,IMP.pmi.io.crosslink.CrossLinkDataBase):
70 raise TypeError(
"CrossLinkingMassSpectrometryRestraint: CrossLinkDataBase should be an IMP.pmi.io.crosslink.CrossLinkDataBase object")
71 self.CrossLinkDataBase = CrossLinkDataBase
73 if resolution==0
or resolution
is None:
74 raise Exception(
"You must pass a resolution and it can't be zero")
76 indb = open(
"included." + filelabel +
".xl.db",
"w")
77 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
78 midb = open(
"missing." + filelabel +
".xl.db",
"w")
88 self.linear.set_slope(0.0)
92 self.psi_is_sampled =
True
93 self.sigma_is_sampled =
True
94 self.psi_dictionary={}
95 self.sigma_dictionary={}
97 self.outputlevel =
"low"
103 copies_to_add = defaultdict(int)
104 print(
'gathering copies')
105 for xlid
in self.CrossLinkDataBase.xlid_iterator():
106 for xl
in self.CrossLinkDataBase[xlid]:
107 r1 = xl[self.CrossLinkDataBase.residue1_key]
108 c1 = xl[self.CrossLinkDataBase.protein1_key]
109 r2 = xl[self.CrossLinkDataBase.residue2_key]
110 c2 = xl[self.CrossLinkDataBase.protein2_key]
111 for c,r
in ((c1,r1),(c2,r2)):
112 if c
in copies_to_add:
118 resolution=resolution).get_selected_particles()
120 copies_to_add[c] = len(sel)-1
122 for molname
in copies_to_add:
123 if copies_to_add[molname]==0:
126 self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein1_key,molname+
'.0',fo1)
128 self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein2_key,molname+
'.0',fo2)
129 for ncopy
in range(copies_to_add[molname]):
130 self.CrossLinkDataBase.clone_protein(
'%s.0'%molname,
'%s.%i'%(molname,ncopy+1))
131 print(
'done pmi2 prelims')
133 for xlid
in self.CrossLinkDataBase.xlid_iterator():
134 new_contribution=
True
135 for xl
in self.CrossLinkDataBase[xlid]:
137 r1 = xl[self.CrossLinkDataBase.residue1_key]
138 c1 = xl[self.CrossLinkDataBase.protein1_key]
139 r2 = xl[self.CrossLinkDataBase.residue2_key]
140 c2 = xl[self.CrossLinkDataBase.protein2_key]
143 iterlist = range(len(IMP.atom.get_by_type(root_hier,IMP.atom.STATE_TYPE)))
145 iterlist = representations
146 for nstate, r
in enumerate(iterlist):
148 xl[self.CrossLinkDataBase.state_key]=nstate
149 xl[self.CrossLinkDataBase.data_set_name_key]=self.label
157 name1,copy1 = c1.split(
'.')
159 name2,copy2 = c2.split(
'.')
163 copy_index=int(copy1),
165 resolution=resolution).get_selected_particles()
169 copy_index=int(copy2),
171 resolution=resolution).get_selected_particles()
178 resolution=resolution,
180 name_is_ambiguous=
False,
184 resolution=resolution,
186 name_is_ambiguous=
False,
190 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
192 print(
"CrossLinkingMassSpectrometryRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
193 midb.write(str(xl) +
"\n")
197 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
199 print(
"CrossLinkingMassSpectrometryRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
200 midb.write(str(xl) +
"\n")
206 if p1 == p2
and r1 == r2:
207 print(
"CrossLinkingMassSpectrometryRestraint: WARNING> same particle and same residue, skipping cross-link")
211 print(
"generating a new crosslink restraint")
212 new_contribution=
False
217 restraints.append(dr)
220 if self.CrossLinkDataBase.sigma1_key
not in xl.keys():
222 xl[self.CrossLinkDataBase.sigma1_key]=sigma1name
224 sigma1name=xl[self.CrossLinkDataBase.sigma1_key]
227 if self.CrossLinkDataBase.sigma2_key
not in xl.keys():
229 xl[self.CrossLinkDataBase.sigma2_key]=sigma2name
231 sigma2name=xl[self.CrossLinkDataBase.sigma2_key]
234 if self.CrossLinkDataBase.psi_key
not in xl.keys():
236 xl[self.CrossLinkDataBase.psi_key]=psiname
238 psiname=xl[self.CrossLinkDataBase.psi_key]
243 p1i = p1.get_particle_index()
245 p2i = p2.get_particle_index()
248 xl[
"Particle_sigma1"]=sigma1
250 xl[
"Particle_sigma2"]=sigma2
252 xl[
"Particle_psi"]=psi
254 dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
257 print(
"--------------")
258 print(
"CrossLinkingMassSpectrometryRestraint: generating cross-link restraint between")
259 print(
"CrossLinkingMassSpectrometryRestraint: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
260 print(
"CrossLinkingMassSpectrometryRestraint: with sigma1 %s sigma2 %s psi %s" % (sigma1name, sigma2name, psiname))
261 print(
"CrossLinkingMassSpectrometryRestraint: between particles %s and %s" % (p1.get_name(), p2.get_name()))
262 print(
"==========================================\n")
269 xl[
"IntraRigidBody"]=
True
271 xl[
"IntraRigidBody"]=
False
273 xl_label=self.CrossLinkDataBase.get_short_cross_link_string(xl)
274 xl[
"ShortLabel"]=xl_label
275 dr.set_name(xl_label)
279 pr.set_name(xl_label)
280 self.rslin.add_restraint(pr)
282 self.xl_list.append(xl)
284 indb.write(str(xl) +
"\n")
286 if len(self.xl_list) == 0:
287 raise SystemError(
"CrossLinkingMassSpectrometryRestraint: no crosslink was constructed")
290 self.rs.add_restraint(lw)
293 """ Add the restraint to the model so that it is evaluated """
300 """ get the hierarchy """
304 """ get the restraint set """
308 """ get the restraint set (redundant with get_restraint_sets)"""
312 """ get the dummy restraints to be displayed in the rmf file """
316 """ get the restraints in a list """
318 for r
in self.rs.get_restraints():
319 rlist.append(IMP.core.PairRestraint.get_from(r))
323 """ Get a list of tuples containing the particle pairs """
325 for i
in range(len(self.pairs)):
326 p0 = self.pairs[i][0]
327 p1 = self.pairs[i][1]
328 ppairs.append((p0, p1))
332 """ Set the output level of the output """
333 self.outputlevel = level
336 """ Switch on/off the sampling of psi particles """
337 self.psi_is_sampled = is_sampled
340 """ Switch on/off the sampling of sigma particles """
341 self.sigma_is_sampled = is_sampled
345 """ This is called internally. Creates a nuisance
346 on the structural uncertainty """
347 if name
in self.sigma_dictionary:
348 return self.sigma_dictionary[name][0]
351 sigmaminnuis = 0.0000001
352 sigmamaxnuis = 1000.0
356 sigma = IMP.pmi.tools.SetupNuisance(self.m, sigmainit,
357 sigmaminnuis, sigmamaxnuis, self.sigma_is_sampled).get_particle()
358 self.sigma_dictionary[name] = (
361 self.sigma_is_sampled)
362 self.rssig.add_restraint(
372 """ This is called internally. Creates a nuisance
373 on the data uncertainty """
374 if name
in self.psi_dictionary:
375 return self.psi_dictionary[name][0]
378 psiminnuis = 0.0000001
379 psimaxnuis = 0.4999999
383 psi = IMP.pmi.tools.SetupNuisance(self.m, psiinit,
384 psiminnuis, psimaxnuis,
385 self.psi_is_sampled).get_particle()
386 self.psi_dictionary[name] = (
391 self.rspsi.add_restraint(
403 """ Set the restraint output label """
407 """ Get the output of the restraint to be used by the IMP.pmi.output object"""
411 score = self.rs.unprotected_evaluate(
None)
412 output[
"_TotalScore"] = str(score)
413 output[
"CrossLinkingMassSpectrometryRestraint_Data_Score_" + self.label] = str(score)
414 output[
"CrossLinkingMassSpectrometryRestraint_PriorSig_Score_" +
415 self.label] = self.rssig.unprotected_evaluate(
None)
416 output[
"CrossLinkingMassSpectrometryRestraint_PriorPsi_Score_" +
417 self.label] = self.rspsi.unprotected_evaluate(
None)
418 output[
"CrossLinkingMassSpectrometryRestraint_Linear_Score_" +
419 self.label] = self.rslin.unprotected_evaluate(
None)
420 for xl
in self.xl_list:
422 xl_label=xl[
"ShortLabel"]
426 output[
"CrossLinkingMassSpectrometryRestraint_Score_" +
427 xl_label] = str(-log(ln.unprotected_evaluate(
None)))
431 output[
"CrossLinkingMassSpectrometryRestraint_Distance_" +
435 for psiname
in self.psi_dictionary:
436 output[
"CrossLinkingMassSpectrometryRestraint_Psi_" +
437 str(psiname) +
"_" + self.label] = str(self.psi_dictionary[psiname][0].get_scale())
439 for sigmaname
in self.sigma_dictionary:
440 output[
"CrossLinkingMassSpectrometryRestraint_Sigma_" +
441 str(sigmaname) +
"_" + self.label] = str(self.sigma_dictionary[sigmaname][0].get_scale())
447 """ Get the particles to be sampled by the IMP.pmi.sampler object """
449 if self.sigma_is_sampled:
450 for sigmaname
in self.sigma_dictionary:
451 ps[
"Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" + str(sigmaname) +
"_" + self.label] =\
452 ([self.sigma_dictionary[sigmaname][0]],
453 self.sigma_dictionary[sigmaname][1])
455 if self.psi_is_sampled:
456 for psiname
in self.psi_dictionary:
457 ps[
"Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
458 str(psiname) +
"_" + self.label] =\
459 ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
465 """Setup cross-link distance restraints at atomic level
466 The "atomic" aspect is that it models the particle uncertainty with a Gaussian
467 The noise in the data and the structural uncertainty of cross-linked amino-acids
468 is inferred using Bayes' theory of probability
469 \note Wraps an IMP::isd::AtomicCrossLinkMSRestraint
479 nuisances_are_optimized=
True,
485 Automatically creates one "sigma" per crosslinked residue and one "psis" per pair.
486 Other nuisance options are available.
487 \note Will return an error if the data+extra_sel don't specify two particles per XL pair.
488 @param root_hier The root hierarchy on which you'll do selection
489 @param xldb CrossLinkDataBase object
490 @param atom_type Can either be "NZ" or "CA"
491 @param length The XL linker length
492 @param slope Linear term to add to the restraint function to help when far away
493 @param nstates The number of states to model. Defaults to the number of states in root.
494 @param label The output label for the restraint
495 @param nuisances_are_optimized Whether to optimize nuisances
496 @param sigma_init The initial value for all the sigmas
497 @param psi_init The initial value for all the psis
498 @param one_psi Use a single psi for all restraints (if False, creates one per XL)
499 @param filelabel automatically generated file containing missing/included/excluded
500 cross-links will be labeled using this text
504 self.root = root_hier
505 self.mdl = self.root.get_model()
510 self.sigma_is_sampled = nuisances_are_optimized
511 self.psi_is_sampled = nuisances_are_optimized
513 if atom_type
not in(
"NZ",
"CA"):
514 raise Exception(
"AtomicXLRestraint: atom_type must be NZ or CA")
515 self.atom_type = atom_type
516 self.nstates = nstates
518 self.nstates = len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE))
519 elif nstates!=len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE)):
520 print(
"Warning: nstates is not the same as the number of states in root")
528 self.psi_dictionary = {}
529 self.sigma_dictionary = {}
531 self.particles=defaultdict(set)
532 self.one_psi = one_psi
533 self.bonded_pairs = []
535 print(
'creating a single psi for all XLs')
537 print(
'creating one psi for each XL')
540 if filelabel
is not None:
541 indb = open(
"included." + filelabel +
".xl.db",
"w")
542 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
543 midb = open(
"missing." + filelabel +
".xl.db",
"w")
549 # read ahead to get the number of XL's per residue
550 num_xls_per_res=defaultdict(int)
551 for unique_id in data:
552 for nstate in range(self.nstates):
553 for xl in data[unique_id]:
554 num_xls_per_res[str(xl.r1)]+=1
555 num_xls_per_res[str(xl.r2)]+=1
557 # Will setup two sigmas based on promiscuity of the residue
559 self.sig_low = setup_nuisance(self.mdl,self.rs_nuis,init_val=sigma_init,min_val=1.0,
560 max_val=100.0,is_opt=self.nuis_opt)
561 self.sig_high = setup_nuisance(self.mdl,self.rs_nuis,init_val=sigma_init,min_val=1.0,
562 max_val=100.0,is_opt=self.nuis_opt)
564 self._create_sigma(
'sigma',sigma_init)
566 self._create_psi(
'psi',psi_init)
568 for xlid
in self.xldb.xlid_iterator():
569 self._create_psi(xlid,psi_init)
573 for xlid
in self.xldb.xlid_iterator():
576 psip = self.psi_dictionary[
'psi'][0].get_particle_index()
578 psip = self.psi_dictionary[unique_id][0].get_particle_index()
587 for nstate
in self.nstates:
588 for xl
in self.xldb[xlid]:
589 r1 = xl[self.xldb.residue1_key]
590 c1 = xl[self.xldb.protein1_key].strip()
591 r2 = xl[self.xldb.residue2_key]
592 c2 = xl[self.xldb.protein2_key].strip()
599 residue_index=r1).get_selected_particles()
604 residue_index=r2).get_selected_particles()
606 print(
"AtomicXLRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
607 if filelabel
is not None:
608 midb.write(str(xl) +
"\n")
612 print(
"AtomicXLRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
613 if filelabel
is not None:
614 midb.write(str(xl) +
"\n")
620 num1=num_xls_per_res[str(xl.r1)]
621 num2=num_xls_per_res[str(xl.r2)]
622 if num1<sig_threshold:
626 if num2<sig_threshold:
631 sig1 = self.sigma_dictionary[
'sigma'][0]
632 sig2 = self.sigma_dictionary[
'sigma'][0]
635 for p1,p2
in itertools.product(ps1,ps2):
638 self.particles[nstate]|=set((p1,p2))
639 r.add_contribution([p1.get_index(),p2.get_index()],
640 [sig1.get_particle_index(),sig2.get_particle_index()])
643 if num_contributions>0:
644 print(
'XL:',xlid,
'num contributions:',num_contributions)
647 raise Exception(
"You didn't create any XL restraints")
648 print(
'created',len(xlrs),
'XL restraints')
651 def set_weight(self,weight):
653 self.rs.set_weight(weight)
655 def set_label(self, label):
658 def add_to_model(self):
663 def get_hierarchy(self):
666 def get_restraint_set(self):
669 def _create_sigma(self, name,sigmainit):
670 """ This is called internally. Creates a nuisance
671 on the structural uncertainty """
672 if name
in self.sigma_dictionary:
673 return self.sigma_dictionary[name][0]
675 sigmaminnuis = 0.0000001
676 sigmamaxnuis = 1000.0
680 sigma = IMP.pmi.tools.SetupNuisance(self.mdl,
684 self.sigma_is_sampled).get_particle()
685 self.sigma_dictionary[name] = (
688 self.sigma_is_sampled)
689 self.rs_sig.add_restraint(
698 def _create_psi(self, name,psiinit):
699 """ This is called internally. Creates a nuisance
700 on the data uncertainty """
701 if name
in self.psi_dictionary:
702 return self.psi_dictionary[name][0]
704 psiminnuis = 0.0000001
705 psimaxnuis = 0.4999999
709 psi = IMP.pmi.tools.SetupNuisance(self.mdl,
713 self.psi_is_sampled).get_particle()
714 self.psi_dictionary[name] = (
719 self.rs_psi.add_restraint(
731 """ create dummy harmonic restraints for each XL but don't add to model
732 Makes it easy to see each contribution to each XL in RMF
734 class MyGetRestraint(object):
737 def get_restraint_for_rmf(self):
743 for nxl
in range(self.get_number_of_restraints()):
744 xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
746 for ncontr
in range(xl.get_number_of_contributions()):
747 ps=xl.get_contribution(ncontr)
749 'xl%i_contr%i'%(nxl,ncontr))
751 dummy_rs.append(MyGetRestraint(rs))
756 """ Get the particles to be sampled by the IMP.pmi.sampler object """
758 if self.sigma_is_sampled:
759 for sigmaname
in self.sigma_dictionary:
760 ps[
"Nuisances_AtomicCrossLinkingMSRestraint_Sigma_" + str(sigmaname) +
"_" + self.label] = \
761 ([self.sigma_dictionary[sigmaname][0]],
762 self.sigma_dictionary[sigmaname][1])
763 if self.psi_is_sampled:
764 for psiname
in self.psi_dictionary:
765 ps[
"Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
766 str(psiname) +
"_" + self.label] =\
767 ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
770 def get_bonded_pairs(self):
771 return self.bonded_pairs
773 def get_number_of_restraints(self):
774 return self.rs.get_number_of_restraints()
777 return 'XL restraint with '+str(len(self.rs.get_restraint(0).get_number_of_restraints())) \
781 """Read a stat file and load all the sigmas.
782 This is potentially quite stupid.
783 It's also a hack since the sigmas should be stored in the RMF file.
784 Also, requires one sigma and one psi for ALL XLs.
787 sig_val = float(subprocess.check_output([
"process_output.py",
"-f",in_fn,
788 "-s",
"AtomicXLRestraint_sigma"]).split(
'\n>')[1+nframe])
789 psi_val = float(subprocess.check_output([
"process_output.py",
"-f",in_fn,
790 "-s",
"AtomicXLRestraint_psi"]).split(
'\n>')[1+nframe])
791 for nxl
in range(self.rs.get_number_of_restraints()):
792 xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
795 for contr
in range(xl.get_number_of_contributions()):
796 sig1,sig2=xl.get_contribution_sigmas(contr)
799 print(
'loaded nuisances from file')
802 max_prob_for_violation=0.1,
803 min_dist_for_violation=1e9,
805 limit_to_chains=
None,
807 """Create CMM files, one for each state, of all crosslinks.
808 will draw in GREEN if non-violated in all states (or if only one state)
809 will draw in PURPLE if non-violated only in a subset of states (draws nothing elsewhere)
810 will draw in RED in ALL states if all violated
811 (if only one state, you'll only see green and red)
813 @param out_prefix Output xlink files prefix
814 @param max_prob_for_violation It's a violation if the probability is below this
815 @param min_dist_for_violation It's a violation if the min dist is above this
816 @param coarsen Use CA positions
817 @param limit_to_chains Try to visualize just these chains
818 @param exclude_to_chains Try to NOT visualize these chains
820 print(
'going to calculate violations and plot CMM files')
822 all_dists = [s[
"low_dist"]
for s
in all_stats]
828 cmds = defaultdict(set)
829 for nstate
in range(self.nstates):
830 outf=open(out_prefix+str(nstate)+
'.cmm',
'w')
831 outf.write(
'<marker_set name="xlinks_state%i"> \n' % nstate)
834 print(
'will limit to',limit_to_chains)
835 print(
'will exclude',exclude_chains)
840 for nxl
in range(self.rs.get_number_of_restraints()):
844 for nstate
in range(self.nstates):
845 prob = state_info[nstate][nxl][
"prob"]
846 low_dist = state_info[nstate][nxl][
"low_dist"]
847 if prob<max_prob_for_violation
or low_dist>min_dist_for_violation:
855 if len(npass)==self.nstates:
857 elif len(nviol)==self.nstates:
861 print(nxl,
'state dists:',[state_info[nstate][nxl][
"low_dist"]
for nstate
in range(self.nstates)],
862 'viol states:',nviol,
'all viol?',all_viol)
863 for nstate
in range(self.nstates):
865 r=0.365; g=0.933; b=0.365;
868 r=0.980; g=0.302; b=0.247;
875 r=0.365; g=0.933; b=0.365;
877 pp = state_info[nstate][nxl][
"low_pp"]
886 cmds[nstate].add((ch1,r1))
887 cmds[nstate].add((ch2,r2))
889 outf = out_fns[nstate]
891 outf.write(
'<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
892 'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,c1[0],c1[1],c1[2],r,g,b))
893 outf.write(
'<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
894 'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv+1,c2[0],c2[1],c2[2],r,g,b))
895 outf.write(
'<link id1= "%d" id2="%d" radius="0.8" '
896 'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,nv+1,r,g,b))
899 for nstate
in range(self.nstates):
900 out_fns[nstate].write(
'</marker_set>\n')
901 out_fns[nstate].close()
903 for ch,r
in cmds[nstate]:
904 cmd+=
'#%i:%i.%s '%(nstate,r,ch)
908 def _get_contribution_info(self,xl,ncontr,use_CA=False):
909 """Return the particles at that contribution. If requested will return CA's instead"""
910 idx1=xl.get_contribution(ncontr)[0]
911 idx2=xl.get_contribution(ncontr)[1]
919 return idx1,idx2,dist
921 def get_best_stats(self,limit_to_state=None,limit_to_chains=None,exclude_chains='',use_CA=False):
922 ''' return the probability, best distance, two coords, and possibly the psi for each xl
923 @param limit_to_state Only examine contributions from one state
924 @param limit_to_chains Returns the particles for certain "easy to visualize" chains
925 @param exclude_chains Even if you limit, don't let one end be in this list.
926 Only works if you also limit chains
927 @param use_CA Limit to CA particles
930 for nxl
in range(self.rs.get_number_of_restraints()):
932 xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
939 for contr
in range(xl.get_number_of_contributions()):
940 pp = xl.get_contribution(contr)
947 if limit_to_state
is not None:
949 if nstate!=limit_to_state:
951 state_contrs.append(contr)
954 if limit_to_chains
is not None:
957 if (c1
in limit_to_chains
or c2
in limit_to_chains)
and (
958 c1
not in exclude_chains
and c2
not in exclude_chains):
959 if dist<low_dist_lim:
966 if limit_to_state
is not None:
967 this_info[
"prob"] = xl.evaluate_for_contributions(state_contrs,
None)
969 this_info[
"prob"] = xl.unprotected_evaluate(
None)
970 if limit_to_chains
is not None:
971 this_info[
"low_pp"] = low_pp_lim
973 this_info[
"low_pp"] = low_pp
975 this_info[
"low_dist"] = low_dist
978 this_info[
"psi"] = pval
979 ret.append(this_info)
982 def print_stats(self):
985 for nxl,s
in enumerate(stats):
990 def get_output(self):
993 score = self.weight * self.rs.unprotected_evaluate(
None)
994 output[
"_TotalScore"] = str(score)
995 output[
"AtomicXLRestraint" + self.label] = str(score)
999 output[
"AtomicXLRestraint_sigma"] = self.rs_sig.unprotected_evaluate(
None)
1000 output[
"AtomicXLRestraint_psi"] = self.rs_psi.unprotected_evaluate(
None)
1008 for nxl,s
in enumerate(stats):
1009 if s[
'low_dist']>20.0:
1011 output[
"AtomicXLRestraint_%i_%s"%(nxl,
"Prob")]=str(s[
'prob'])
1012 output[
"AtomicXLRestraint_%i_%s"%(nxl,
"BestDist")]=str(s[
'low_dist'])
1013 if not self.one_psi:
1014 output[
"AtomicXLRestraint_%i_%s"%(nxl,
"psi")]=str(s[
'psi'])
1015 output[
"AtomicXLRestraint_NumViol"] = str(bad_count)
1021 """This restraint allows ambiguous crosslinking between multiple copies
1022 it is a variant of the SimplifiedCrossLinkMS
1032 self.m = representation.prot.get_model()
1038 self.outputlevel =
"low"
1039 self.expdistance = expdistance
1040 self.strength = strength
1042 fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1046 tokens = line.split()
1048 if (tokens[0] ==
"#"):
1057 resolution=resolution,
1059 name_is_ambiguous=
True,
1061 hrc1 = [representation.hier_db.particle_to_name[p]
for p
in ps1]
1063 print(
"ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1068 resolution=resolution,
1070 name_is_ambiguous=
True,
1072 hrc2 = [representation.hier_db.particle_to_name[p]
for p
in ps2]
1074 print(
"ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1081 if self.strength
is None:
1090 rad1 = rad1 / len(ps1)
1091 rad2 = rad2 / len(ps2)
1093 self.strength = 1 / (rad1 ** 2 + rad2 ** 2)
1101 self.rs.add_restraint(cr)
1102 self.pairs.append((ps1, hrc1, c1, r1, ps2, hrc2, c2, r2, cr))
1115 d1.set_radius(uncertainty1)
1116 d2.set_radius(uncertainty2)
1120 strength = 1 / (uncertainty1 ** 2 + uncertainty2 ** 2)
1127 for i
in range(npoints):
1131 scores.append(cr.unprotected_evaluate(
None))
1132 IMP.pmi.output.plot_xy_data(dists, scores)
1134 def set_label(self, label):
1136 self.rs.set_name(label)
1137 for r
in self.rs.get_restraints():
1140 def add_to_model(self):
1143 def get_restraint_sets(self):
1146 def get_restraint(self):
1149 def set_output_level(self, level="low"):
1151 self.outputlevel = level
1153 def set_weight(self, weight):
1154 self.weight = weight
1155 self.rs.set_weight(weight)
1157 def get_output(self):
1163 score = self.weight * self.rs.unprotected_evaluate(
None)
1164 output[
"_TotalScore"] = str(score)
1165 output[
"ConnectivityCrossLinkMS_Score_" + self.label] = str(score)
1166 for n, p
in enumerate(self.pairs):
1177 for n1, p1
in enumerate(ps1):
1180 for n2, p2
in enumerate(ps2):
1184 label = str(r1) +
":" + name1 +
"_" + str(r2) +
":" + name2
1185 output[
"ConnectivityCrossLinkMS_Distance_" +
1188 label = str(r1) +
":" + c1 +
"_" + str(r2) +
":" + c2
1189 output[
"ConnectivityCrossLinkMS_Score_" +
1190 label] = str(self.weight * cr.unprotected_evaluate(
None))
1194 @
IMP.deprecated_object(
"2.5",
"Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1195 class SimplifiedCrossLinkMS(object):
1205 spheredistancepairscore=
True):
1210 if columnmapping
is None:
1212 columnmapping[
"Protein1"] = 0
1213 columnmapping[
"Protein2"] = 1
1214 columnmapping[
"Residue1"] = 2
1215 columnmapping[
"Residue2"] = 3
1217 self.m = representation.prot.get_model()
1222 self.already_added_pairs = {}
1224 self.outputlevel =
"low"
1225 self.expdistance = expdistance
1226 self.strength = strength
1227 self.truncated = truncated
1228 self.spheredistancepairscore = spheredistancepairscore
1233 protein1 = columnmapping[
"Protein1"]
1234 protein2 = columnmapping[
"Protein2"]
1235 residue1 = columnmapping[
"Residue1"]
1236 residue2 = columnmapping[
"Residue2"]
1238 fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1242 tokens = line.split()
1244 if (tokens[0] ==
"#"):
1246 r1 = int(tokens[residue1])
1247 c1 = tokens[protein1]
1248 r2 = int(tokens[residue2])
1249 c2 = tokens[protein2]
1253 resolution=resolution,
1255 name_is_ambiguous=
False,
1259 resolution=resolution,
1261 name_is_ambiguous=
False,
1266 "residue %d of chain %s selects multiple particles; "
1268 % (r1, c1,
"".join(p.get_name()
for p
in ps1)))
1270 print(
"SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1275 "residue %d of chain %s selects multiple particles; "
1277 % (r2, c2,
"".join(p.get_name()
for p
in ps2)))
1279 print(
"SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1285 if (p1, p2)
in self.already_added_pairs:
1286 dr = self.already_added_pairs[(p1, p2)]
1287 weight = dr.get_weight()
1288 dr.set_weight(weight + 1.0)
1289 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))
1295 limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
1306 if self.spheredistancepairscore:
1311 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2))
1313 self.rs.add_restraint(dr)
1314 self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
1315 self.already_added_pairs[(p1, p2)] = dr
1316 self.already_added_pairs[(p2, p1)] = dr
1318 def set_label(self, label):
1321 def add_to_model(self):
1324 def get_restraint_sets(self):
1327 def get_restraint(self):
1330 def get_restraints(self):
1332 for r
in self.rs.get_restraints():
1333 rlist.append(IMP.core.PairRestraint.get_from(r))
1336 def get_particle_pairs(self):
1338 for i
in range(len(self.pairs)):
1339 p0 = self.pairs[i][0]
1340 p1 = self.pairs[i][1]
1341 ppairs.append((p0, p1))
1344 def set_output_level(self, level="low"):
1346 self.outputlevel = level
1348 def set_weight(self, weight):
1349 self.weight = weight
1350 self.rs.set_weight(weight)
1352 def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
1358 d1.set_radius(radius1)
1359 d2.set_radius(radius2)
1361 limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
1376 for i
in range(npoints):
1380 scores.append(dr.unprotected_evaluate(
None))
1381 IMP.pmi.output.plot_xy_data(dists, scores)
1383 def get_output(self):
1389 score = self.weight * self.rs.unprotected_evaluate(
None)
1390 output[
"_TotalScore"] = str(score)
1391 output[
"SimplifiedCrossLinkMS_Score_" + self.label] = str(score)
1392 for i
in range(len(self.pairs)):
1394 p0 = self.pairs[i][0]
1395 p1 = self.pairs[i][1]
1396 crosslinker =
'standard'
1397 ln = self.pairs[i][2]
1398 resid1 = self.pairs[i][3]
1399 chain1 = self.pairs[i][4]
1400 resid2 = self.pairs[i][5]
1401 chain2 = self.pairs[i][6]
1403 label = str(resid1) +
":" + chain1 +
"_" + \
1404 str(resid2) +
":" + chain2
1405 output[
"SimplifiedCrossLinkMS_Score_" + crosslinker +
"_" +
1406 label] = str(self.weight * ln.unprotected_evaluate(
None))
1408 if self.spheredistancepairscore:
1414 output[
"SimplifiedCrossLinkMS_Distance_" +
1421 @
IMP.deprecated_object(
"2.5",
"Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1422 class SigmoidalCrossLinkMS(object):
1424 self, representation, restraints_file, inflection, slope, amplitude,
1425 linear_slope, resolution=
None, columnmapping=
None, csvfile=
False,
1426 filters=
None, filelabel=
"None"):
1432 if columnmapping
is None:
1434 columnmapping[
"Protein1"] = 0
1435 columnmapping[
"Protein2"] = 1
1436 columnmapping[
"Residue1"] = 2
1437 columnmapping[
"Residue2"] = 3
1443 db = tools.get_db_from_csv(restraints_file)
1445 db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1447 self.m = representation.prot.get_model()
1453 self.already_added_pairs = {}
1454 self.inflection = inflection
1456 self.amplitude = amplitude
1457 self.linear_slope = linear_slope
1459 self.outputlevel =
"low"
1464 protein1 = columnmapping[
"Protein1"]
1465 protein2 = columnmapping[
"Protein2"]
1466 residue1 = columnmapping[
"Residue1"]
1467 residue2 = columnmapping[
"Residue2"]
1469 indb = open(
"included." + filelabel +
".xl.db",
"w")
1470 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
1471 midb = open(
"missing." + filelabel +
".xl.db",
"w")
1475 tokens = entry.split()
1477 if (tokens[0] ==
"#"):
1479 r1 = int(tokens[residue1])
1480 c1 = tokens[protein1]
1481 r2 = int(tokens[residue2])
1482 c2 = tokens[protein2]
1484 if filters
is not None:
1485 if eval(tools.cross_link_db_filter_parser(filters)) ==
False:
1486 exdb.write(str(entry) +
"\n")
1488 indb.write(str(entry) +
"\n")
1489 r1 = int(entry[residue1])
1490 c1 = entry[protein1]
1491 r2 = int(entry[residue2])
1492 c2 = entry[protein2]
1496 resolution=resolution,
1498 name_is_ambiguous=
False,
1502 resolution=resolution,
1504 name_is_ambiguous=
False,
1508 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
1510 print(
"SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1511 midb.write(str(entry) +
"\n")
1515 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
1517 print(
"SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1518 midb.write(str(entry) +
"\n")
1524 if (p1, p2)
in self.already_added_pairs:
1525 dr = self.already_added_pairs[(p1, p2)]
1526 weight = dr.get_weight()
1527 dr.increment_amplitude(amplitude)
1528 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()))
1529 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2)
1530 +
"-ampl:" + str(dr.get_amplitude()))
1543 dr.set_name(c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2)
1544 +
"-ampl:" + str(dr.get_amplitude()))
1546 self.rs.add_restraint(dr)
1548 self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
1549 self.already_added_pairs[(p1, p2)] = dr
1550 self.already_added_pairs[(p2, p1)] = dr
1552 def set_label(self, label):
1555 def add_to_model(self):
1558 def get_restraint_sets(self):
1561 def get_restraint(self):
1564 def get_restraints(self):
1566 for r
in self.rs.get_restraints():
1567 rlist.append(IMP.core.PairRestraint.get_from(r))
1570 def get_particle_pairs(self):
1572 for i
in range(len(self.pairs)):
1573 p0 = self.pairs[i][0]
1574 p1 = self.pairs[i][1]
1575 ppairs.append((p0, p1))
1578 def set_output_level(self, level="low"):
1580 self.outputlevel = level
1582 def set_weight(self, weight):
1583 self.weight = weight
1584 self.rs.set_weight(weight)
1586 def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
1591 d1.set_radius(radius1)
1592 d2.set_radius(radius2)
1603 for i
in range(npoints):
1607 scores.append(dr.unprotected_evaluate(
None))
1608 IMP.pmi.output.plot_xy_data(dists, scores)
1610 def get_output(self):
1616 score = self.weight * self.rs.unprotected_evaluate(
None)
1617 output[
"_TotalScore"] = str(score)
1618 output[
"SigmoidalCrossLinkMS_Score_" + self.label] = str(score)
1619 for i
in range(len(self.pairs)):
1621 p0 = self.pairs[i][0]
1622 p1 = self.pairs[i][1]
1623 crosslinker =
'standard'
1624 ln = self.pairs[i][2]
1625 resid1 = self.pairs[i][3]
1626 chain1 = self.pairs[i][4]
1627 resid2 = self.pairs[i][5]
1628 chain2 = self.pairs[i][6]
1630 label = str(resid1) +
":" + chain1 +
"_" + \
1631 str(resid2) +
":" + chain2
1632 output[
"SigmoidalCrossLinkMS_Score_" + crosslinker +
"_" +
1633 label +
"_" + self.label] = str(self.weight * ln.unprotected_evaluate(
None))
1636 output[
"SigmoidalCrossLinkMS_Distance_" +
1641 @
IMP.deprecated_object(
"2.5",
"Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1642 class ISDCrossLinkMS(IMP.pmi.restraints._NuisancesBase):
1643 def __init__(self, representation,
1663 automatic_sigma_classification=
False,
1664 attributes_for_label=
None):
1674 if type(representation) != list:
1675 representations = [representation]
1677 representations = representation
1679 if columnmapping
is None:
1681 columnmapping[
"Protein1"] = 0
1682 columnmapping[
"Protein2"] = 1
1683 columnmapping[
"Residue1"] = 2
1684 columnmapping[
"Residue2"] = 3
1685 columnmapping[
"IDScore"] = 4
1686 columnmapping[
"XLUniqueID"] = 5
1692 db = IMP.pmi.tools.get_db_from_csv(restraints_file)
1694 db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1696 indb = open(
"included." + filelabel +
".xl.db",
"w")
1697 exdb = open(
"excluded." + filelabel +
".xl.db",
"w")
1698 midb = open(
"missing." + filelabel +
".xl.db",
"w")
1700 self.m = representations[0].prot.get_model()
1709 self.sigma_dictionary = {}
1710 self.psi_dictionary = {}
1711 self.psi_is_sampled =
True
1712 self.sigma_is_sampled =
True
1717 self.ids_map = IMP.pmi.tools.map()
1718 self.ids_map.set_map_element(20.0, 0.05)
1719 self.ids_map.set_map_element(65.0, 0.01)
1721 self.ids_map = ids_map
1723 if radius_map
is None:
1724 self.radius_map = IMP.pmi.tools.map()
1725 if automatic_sigma_classification:
1726 self.radius_map.set_map_element(10, 10)
1727 self.radius_map.set_map_element(1, 1)
1729 self.radius_map = radius_map
1731 self.outputlevel =
"low"
1735 self.linear.set_slope(slope)
1738 protein1 = columnmapping[
"Protein1"]
1739 protein2 = columnmapping[
"Protein2"]
1740 residue1 = columnmapping[
"Residue1"]
1741 residue2 = columnmapping[
"Residue2"]
1742 idscore = columnmapping[
"IDScore"]
1744 xluniqueid = columnmapping[
"XLUniqueID"]
1754 uniqueid_restraints_map = {}
1756 for nxl, entry
in enumerate(db):
1758 if not jackknifing
is None:
1765 raise NotImplementedError(
"jackknifing not yet implemented")
1768 tokens = entry.split()
1773 if (tokens[0] ==
"#"):
1776 r1 = int(tokens[residue1])
1777 c1 = tokens[protein1]
1778 r2 = int(tokens[residue2])
1779 c2 = tokens[protein2]
1781 if offset_dict
is not None:
1782 if c1
in offset_dict: r1+=offset_dict[c1]
1783 if c2
in offset_dict: r2+=offset_dict[c2]
1785 if rename_dict
is not None:
1786 if c1
in rename_dict: c1=rename_dict[c1]
1787 if c2
in rename_dict: c2=rename_dict[c2]
1792 ids = float(tokens[idscore])
1793 if xluniqueid
is None:
1796 xlid = tokens[xluniqueid]
1798 print(
"this line was not accessible " + str(entry))
1799 if residue1
not in entry: print(str(residue1)+
" keyword not in database")
1800 if residue2
not in entry: print(str(residue2)+
" keyword not in database")
1801 if protein1
not in entry: print(str(protein1)+
" keyword not in database")
1802 if protein2
not in entry: print(str(protein2)+
" keyword not in database")
1803 if idscore
not in entry: print(str(idscore)+
" keyword not in database")
1804 if xluniqueid
not in entry: print(str(xluniqueid)+
" keyword not in database")
1808 if filters
is not None:
1810 exdb.write(str(entry) +
"\n")
1814 r1 = int(entry[residue1])
1815 c1 = entry[protein1]
1816 r2 = int(entry[residue2])
1817 c2 = entry[protein2]
1819 if offset_dict
is not None:
1820 if c1
in offset_dict: r1+=offset_dict[c1]
1821 if c2
in offset_dict: r2+=offset_dict[c2]
1823 if rename_dict
is not None:
1824 if c1
in rename_dict: c1=rename_dict[c1]
1825 if c2
in rename_dict: c2=rename_dict[c2]
1831 ids = float(entry[idscore])
1833 ids = entry[idscore]
1834 if xluniqueid
is None:
1837 xlid = entry[xluniqueid]
1840 print(
"this line was not accessible " + str(entry))
1841 if residue1
not in entry: print(str(residue1)+
" keyword not in database")
1842 if residue2
not in entry: print(str(residue2)+
" keyword not in database")
1843 if protein1
not in entry: print(str(protein1)+
" keyword not in database")
1844 if protein2
not in entry: print(str(protein2)+
" keyword not in database")
1845 if idscore
not in entry: print(str(idscore)+
" keyword not in database")
1846 if xluniqueid
not in entry: print(str(xluniqueid)+
" keyword not in database")
1849 for nstate, r
in enumerate(representations):
1854 resolution=resolution,
1856 name_is_ambiguous=
False,
1860 resolution=resolution,
1862 name_is_ambiguous=
False,
1866 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
1868 print(
"ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1869 midb.write(str(entry) +
"\n")
1873 raise ValueError(
"residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
1875 print(
"ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1876 midb.write(str(entry) +
"\n")
1882 if (p1 == p2)
and (r1 == r2) :
1883 print(
"ISDCrossLinkMS Restraint: WARNING> on the identical bead particles and the identical residues, thus skipping this cross-link.")
1886 if xlid
in uniqueid_restraints_map:
1887 print(
"Appending a crosslink restraint into the uniqueID %s" % str(xlid))
1888 dr = uniqueid_restraints_map[xlid]
1890 print(
"Generating a NEW crosslink restraint with a uniqueID %s" % str(xlid))
1895 restraints.append(dr)
1896 uniqueid_restraints_map[xlid] = dr
1898 mappedr1 = self.radius_map.get_map_element(
1900 sigma1 = self.get_sigma(mappedr1)[0]
1901 mappedr2 = self.radius_map.get_map_element(
1903 sigma2 = self.get_sigma(mappedr2)[0]
1905 psival = self.ids_map.get_map_element(ids)
1906 psi = self.get_psi(psival)[0]
1909 p1i = p1.get_particle_index()
1910 p2i = p2.get_particle_index()
1918 dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
1919 print(
"--------------")
1920 print(
"ISDCrossLinkMS: generating cross-link restraint between")
1921 print(
"ISDCrossLinkMS: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
1922 print(
"ISDCrossLinkMS: with sigma1 %f sigma2 %f psi %s" % (mappedr1, mappedr2, psival))
1923 print(
"ISDCrossLinkMS: between particles %s and %s" % (p1.get_name(), p2.get_name()))
1924 print(
"==========================================\n")
1925 indb.write(str(entry) +
"\n")
1932 xlattribute =
"intrarb"
1934 xlattribute =
"interrb"
1937 if not attributes_for_label
is None:
1938 for a
in attributes_for_label:
1939 xlattribute = xlattribute +
"_" + str(entry[a])
1941 xlattribute = xlattribute +
"-State:" + str(nstate)
1944 xlattribute +
"-" + c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2) +
"_" + self.label)
1949 xlattribute +
"-" + c1 +
":" + str(r1) +
"-" + c2 +
":" + str(r2) +
"_" + self.label)
1950 self.rslin.add_restraint(pr)
1970 self.rs.add_restraint(lw)
1973 def set_slope_linear_term(self, slope):
1974 self.linear.set_slope(slope)
1976 def set_label(self, label):
1979 def add_to_model(self):
1986 def get_hierarchies(self):
1989 def get_restraint_sets(self):
1992 def get_restraint(self):
1995 def get_restraint_for_rmf(self):
1998 def get_restraints(self):
2000 for r
in self.rs.get_restraints():
2001 rlist.append(IMP.core.PairRestraint.get_from(r))
2004 def get_particle_pairs(self):
2006 for i
in range(len(self.pairs)):
2007 p0 = self.pairs[i][0]
2008 p1 = self.pairs[i][1]
2009 ppairs.append((p0, p1))
2012 def set_output_level(self, level="low"):
2014 self.outputlevel = level
2016 def set_psi_is_sampled(self, is_sampled=True):
2017 self.psi_is_sampled = is_sampled
2019 def set_sigma_is_sampled(self, is_sampled=True):
2020 self.sigma_is_sampled = is_sampled
2022 def get_label(self,pairs_index):
2023 resid1 = self.pairs[pairs_index][3]
2024 chain1 = self.pairs[pairs_index][4]
2025 resid2 = self.pairs[pairs_index][5]
2026 chain2 = self.pairs[pairs_index][6]
2027 attribute = self.pairs[pairs_index][7]
2028 rad1 = self.pairs[pairs_index][8]
2029 rad2 = self.pairs[pairs_index][9]
2030 psi = self.pairs[pairs_index][10]
2031 xlid= self.pairs[pairs_index][11]
2032 label = attribute +
"-" + \
2033 str(resid1) +
":" + chain1 +
"_" + str(resid2) +
":" + \
2034 chain2 +
"-" + str(rad1) +
"-" + str(rad2) +
"-" + str(psi)
2037 def write_db(self,filename):
2038 cldb=IMP.pmi.output.CrossLinkIdentifierDatabase()
2040 for pairs_index
in range(len(self.pairs)):
2042 resid1 = self.pairs[pairs_index][3]
2043 chain1 = self.pairs[pairs_index][4]
2044 resid2 = self.pairs[pairs_index][5]
2045 chain2 = self.pairs[pairs_index][6]
2046 attribute = self.pairs[pairs_index][7]
2047 rad1 = self.pairs[pairs_index][8]
2048 rad2 = self.pairs[pairs_index][9]
2049 psi = self.pairs[pairs_index][10]
2050 xlid= self.pairs[pairs_index][11]
2051 nstate=self.pairs[pairs_index][12]
2052 ids=self.pairs[pairs_index][13]
2054 label=self.get_label(pairs_index)
2055 cldb.set_unique_id(label,xlid)
2056 cldb.set_protein1(label,chain1)
2057 cldb.set_protein2(label,chain2)
2058 cldb.set_residue1(label,resid1)
2059 cldb.set_residue2(label,resid2)
2060 cldb.set_idscore(label,ids)
2061 cldb.set_state(label,nstate)
2062 cldb.set_sigma1(label,rad1)
2063 cldb.set_sigma2(label,rad2)
2064 cldb.set_psi(label,psi)
2065 cldb.write(filename)
2068 def get_output(self):
2074 score = self.rs.unprotected_evaluate(
None)
2075 output[
"_TotalScore"] = str(score)
2076 output[
"ISDCrossLinkMS_Data_Score_" + self.label] = str(score)
2077 output[
"ISDCrossLinkMS_PriorSig_Score_" +
2078 self.label] = self.rssig.unprotected_evaluate(
None)
2079 output[
"ISDCrossLinkMS_PriorPsi_Score_" +
2080 self.label] = self.rspsi.unprotected_evaluate(
None)
2081 output[
"ISDCrossLinkMS_Linear_Score_" +
2082 self.label] = self.rslin.unprotected_evaluate(
None)
2083 for i
in range(len(self.pairs)):
2085 label=self.get_label(i)
2086 ln = self.pairs[i][2]
2087 p0 = self.pairs[i][0]
2088 p1 = self.pairs[i][1]
2089 output[
"ISDCrossLinkMS_Score_" +
2090 label +
"_" + self.label] = str(-log(ln.unprotected_evaluate(
None)))
2094 output[
"ISDCrossLinkMS_Distance_" +
2098 for psiindex
in self.psi_dictionary:
2099 output[
"ISDCrossLinkMS_Psi_" +
2100 str(psiindex) +
"_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2102 for resolution
in self.sigma_dictionary:
2103 output[
"ISDCrossLinkMS_Sigma_" +
2104 str(resolution) +
"_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
2109 def get_particles_to_sample(self):
2111 for resolution
in self.sigma_dictionary:
2112 if self.sigma_dictionary[resolution][2]
and self.sigma_is_sampled:
2113 ps[
"Nuisances_ISDCrossLinkMS_Sigma_" + str(resolution) +
"_" + self.label] =\
2114 ([self.sigma_dictionary[resolution][0]],
2115 self.sigma_dictionary[resolution][1])
2116 if self.psi_is_sampled:
2117 for psiindex
in self.psi_dictionary:
2118 if self.psi_dictionary[psiindex][2]:
2119 ps[
"Nuisances_ISDCrossLinkMS_Psi_" +
2120 str(psiindex) +
"_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2124 class CysteineCrossLinkRestraint(object):
2125 def __init__(self, representations, filename, cbeta=False,
2126 betatuple=(0.03, 0.1),
2127 disttuple=(0.0, 25.0, 1000),
2128 omegatuple=(1.0, 1000.0, 50),
2129 sigmatuple=(0.3, 0.3, 1),
2131 sigmaissampled=
False,
2132 weightissampled=
True,
2133 epsilonissampled=
True
2140 self.representations = representations
2141 self.m = self.representations[0].prot.get_model()
2144 self.epsilonmaxtrans = 0.01
2145 self.sigmamaxtrans = 0.1
2146 self.betamaxtrans = 0.01
2147 self.weightmaxtrans = 0.1
2149 self.outputlevel =
"low"
2150 self.betaissampled = betaissampled
2151 self.sigmaissampled = sigmaissampled
2152 self.weightissampled = weightissampled
2153 self.epsilonissampled = epsilonissampled
2155 betalower = betatuple[0]
2156 betaupper = betatuple[1]
2163 self.beta = IMP.pmi.tools.SetupNuisance(
2168 betaissampled).get_particle(
2171 self.sigma = IMP.pmi.tools.SetupNuisance(
2176 sigmaissampled).get_particle()
2178 self.weight = IMP.pmi.tools.SetupWeight(
2180 weightissampled).get_particle(
2184 fl = IMP.pmi.tools.open_file_or_inline_text(filename)
2193 if t[5]
in self.epsilons:
2194 if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
2195 self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
2197 self.epsilons[t[5]] = IMP.pmi.tools.SetupNuisance(self.m,
2198 0.01, 0.01, 1.0 - float(t[4]), epsilonissampled).get_particle()
2199 up = self.epsilons[t[5]].get_upper()
2200 low = self.epsilons[t[5]].get_lower()
2202 self.epsilons[t[5]].set_upper(low)
2204 data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
2208 crossdata = IMP.pmi.tools.get_cross_link_data(
2209 "cysteine",
"cysteine_CA_FES.txt.standard",
2210 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2212 crossdata = IMP.pmi.tools.get_cross_link_data(
2213 "cysteine",
"cysteine_CB_FES.txt.standard",
2214 disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2217 fmod_grid = IMP.pmi.tools.get_grid(0.0, 1.0, 300,
True)
2218 omega2_grid = IMP.pmi.tools.get_log_grid(0.001, 10000.0, 100)
2219 beta_grid = IMP.pmi.tools.get_log_grid(betalower, betaupper, betangrid)
2222 print(
"--------------")
2223 print(
"CysteineCrossLink: attempting to create a restraint " + str(d))
2243 self.epsilons[epslabel],
2249 for i, representation
in enumerate(self.representations):
2256 resolution=1, name=chain1,
2257 name_is_ambiguous=
False, residue=resid1)[0]
2263 resolution=1, name=chain2,
2264 name_is_ambiguous=
False, residue=resid2)[0]
2273 for t
in range(-1, 2):
2275 resolution=1, name=chain1,
2276 name_is_ambiguous=
False, residue=resid1 + t)
2282 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
2285 resolution=1, name=chain2,
2286 name_is_ambiguous=
False, residue=resid2 + t)
2292 print(
"\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
2295 if (p1
is not None and p2
is not None):
2296 ccl.add_contribution(p1, p2)
2300 print(
"Distance_" + str(resid1) +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2,
IMP.core.get_distance(d1, d2))
2303 if (len(p1) == 3
and len(p2) == 3):
2304 p11n = p1[0].get_name()
2305 p12n = p1[1].get_name()
2306 p13n = p1[2].get_name()
2307 p21n = p2[0].get_name()
2308 p22n = p2[1].get_name()
2309 p23n = p2[2].get_name()
2311 print(
"CysteineCrossLink: generating CB cysteine cross-link restraint between")
2312 print(
"CysteineCrossLink: residue %d of chain %s and residue %d of chain %s" % (resid1, chain1, resid2, chain2))
2313 print(
"CysteineCrossLink: between particles %s %s %s and %s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
2315 ccl.add_contribution(p1, p2)
2318 self.rs.add_restraint(ccl)
2319 ccl.set_name(
"CysteineCrossLink_" + str(resid1)
2320 +
"_" + chain1 +
":" + str(resid2) +
"_" + chain2)
2322 def set_label(self, label):
2325 def add_to_model(self):
2328 def get_particles_to_sample(self):
2330 if self.epsilonissampled:
2331 for eps
in self.epsilons.keys():
2332 ps[
"Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
2333 self.label] = ([self.epsilons[eps]], self.epsilonmaxtrans)
2334 if self.betaissampled:
2335 ps[
"Nuisances_CysteineCrossLinkRestraint_beta_" +
2336 self.label] = ([self.beta], self.betamaxtrans)
2337 if self.weightissampled:
2338 ps[
"Weights_CysteineCrossLinkRestraint_" +
2339 self.label] = ([self.weight], self.weightmaxtrans)
2340 if self.sigmaissampled:
2341 ps[
"Nuisances_CysteineCrossLinkRestraint_" +
2342 self.label] = ([self.sigma], self.sigmamaxtrans)
2345 def set_output_level(self, level="low"):
2347 self.outputlevel = level
2349 def get_restraint(self):
2352 def get_sigma(self):
2355 def get_output(self):
2358 score = self.rs.unprotected_evaluate(
None)
2359 output[
"_TotalScore"] = str(score)
2360 output[
"CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
2361 output[
"CysteineCrossLinkRestraint_sigma_" +
2362 self.label] = str(self.sigma.get_scale())
2363 for eps
in self.epsilons.keys():
2364 output[
"CysteineCrossLinkRestraint_epsilon_" + eps +
"_" +
2365 self.label] = str(self.epsilons[eps].get_scale())
2366 output[
"CysteineCrossLinkRestraint_beta_" +
2367 self.label] = str(self.beta.get_scale())
2368 for n
in range(self.weight.get_number_of_states()):
2369 output[
"CysteineCrossLinkRestraint_weights_" +
2370 str(n) +
"_" + self.label] = str(self.weight.get_weight(n))
2372 if self.outputlevel ==
"high":
2373 for rst
in self.rs.get_restraints():
2374 output[
"CysteineCrossLinkRestraint_Total_Frequency_" +
2375 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2376 "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_model_frequency()
2377 output[
"CysteineCrossLinkRestraint_Standard_Error_" +
2378 IMP.isd.CysteineCrossLinkRestraint.get_from(
2381 + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_standard_error()
2382 if len(self.representations) > 1:
2383 for i
in range(len(self.prots)):
2384 output[
"CysteineCrossLinkRestraint_Frequency_Contribution_" +
2385 IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2386 "_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.
def get_restraint
get the restraint set (redundant with get_restraint_sets)
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.
Object used to hold a set of restraints.
Class for storing model, its restraints, constraints, and particles.
def set_output_level
Set the output level of the output.
def set_label
Set the restraint output label.
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_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.
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 model particles.
def add_to_model
Add the restraint to the model so that it is evaluated.
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)
def load_nuisances_from_stat_file
Read a stat file and load all the sigmas.
Select hierarchy particles identified by the biological name.
def set_psi_is_sampled
Switch on/off the sampling of psi particles.
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.