import IMP.atom import IMP.core import IMP.domino import IMP.container import math model=IMP.Model() model.set_log_level(IMP.SILENT) sqrt3=math.sqrt(3.0) # various parameters ds=40.0 # initial bin size niter=3 # number of iterations do_statistics=True # print statistics do_random=True # stochastic filter do_save_ass=False # save assignments to rmf file skip=100 # skip solutions cell_type="hexagon" # cell type side=80.0 # cell size # cell dependent parameters if(cell_type=="rhombus"): num_cells=21 num_copies=2 error_bound=1.45*ds**2 if(cell_type=="hexagon"): num_cells=7 num_copies=6 error_bound=1.45*ds**2 if(cell_type=="square"): num_cells=9 num_copies=6 side=math.sqrt(side*side*1.5*sqrt3) error_bound=ds**2 #print "error bound is", error_bound def do_compress(points, xf, yf): ret=[] for xx in points: ret.append([xx[0]*xf,xx[1]*yf]) return ret def do_shear(points, bh, h): b=bh/h ret=[] for xx in points: ret.append([xx[0]-b*xx[1],xx[1]]) return ret def grid_cell(side, ds, z): bb= IMP.algebra.BoundingBox2D(IMP.algebra.Vector2D(0,0), IMP.algebra.Vector2D(side,side)) cur=IMP.algebra.get_grid_interior_cover_by_spacing(bb, ds) if(cell_type!="square"): cur_comp=do_compress(cur,1.0,sqrt3/2.0) cur=do_shear(cur_comp,side/2.0,side*sqrt3/2.0) positions=[] for xx in cur: positions.append(IMP.algebra.Vector3D(xx[0],xx[1],z)) if(cell_type=="hexagon"): positions2=[] tra=IMP.algebra.Vector3D(0.0,0.0,0.0) for i in range(1,3): rot=IMP.algebra.get_rotation_about_axis(IMP.algebra.Vector3D(0,0,1), float(i)*math.radians(120.0)) tr=IMP.algebra.Transformation3D(rot,tra) for v in positions: positions2.append(tr.get_transformed(v)) positions=positions+positions2 return positions def create_hierarchies(ncells,name): h=[] for i in range(0,ncells): p1= IMP.Particle(model) h1= IMP.atom.Hierarchy.setup_particle(p1) h1.set_name(name+" hierarchy, cell " + str(i)) h.append(h1) return h def create_protein(name, mass, nbeads, copy, start_residue=1, length=-1): if length==-1: length=int(mass*1000.0/110) p=IMP.Particle(model) protein=IMP.atom.Molecule.setup_particle(p) protein.set_name(name) vol=IMP.atom.get_volume_from_mass(1000.0*mass)/float(nbeads) ms=1000.0*mass/float(nbeads) rg=IMP.algebra.get_ball_radius_from_volume_3d(vol) for i in range(0,nbeads): pp=IMP.Particle(model) first=start_residue+i*int(length/nbeads) last= start_residue+(i+1)*int(length/nbeads) dom=IMP.atom.Domain.setup_particle(pp, (first, last)) dom.set_name(name+str(i)+"-"+str(copy)) d=IMP.core.XYZR.setup_particle(pp) d.set_radius(rg) mm=IMP.atom.Mass.setup_particle(pp,ms) protein.add_child(dom) if ( nbeads > 1 and copy==0 ): con=IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c) \ for c in protein.get_children()],1.0) con.set_name("Connectivity Restraint for "+name) model.add_restraint(con) model.set_maximum_score(con, error_bound) return protein def create_merged_protein(name, protein_a, protein_b, copy, dist=-1): p= IMP.Particle(model) h= IMP.atom.Molecule.setup_particle(p) h.set_name(name) if (copy==0 and dist >=0.0): add_internal_restraint(name, protein_a, protein_b, dist) for ps in protein_a.get_leaves(): protein_a.remove_child(IMP.atom.Domain(ps)) h.add_child(IMP.atom.Domain(ps)) for ps in protein_b.get_leaves(): protein_b.remove_child(IMP.atom.Domain(ps)) h.add_child(IMP.atom.Domain(ps)) return h # connect the two proteins with a chain of length dist def add_internal_restraint(name, protein_a, protein_b, dist): ps= get_pair_score((-500.0, dist)) sa=IMP.atom.Selection(hierarchy=protein_a, terminus=IMP.atom.Selection.C) sb=IMP.atom.Selection(hierarchy=protein_b, terminus=IMP.atom.Selection.N) pa=sa.get_selected_particles()[0] pb=sb.get_selected_particles()[0] r= IMP.core.PairRestraint(ps, IMP.ParticlePair(pa, pb)) r.set_name("IR " + name) model.add_restraint(r) model.set_maximum_score(r, error_bound) def add_distance_restraint(name, p1, p2, tbr, dist): ps= get_pair_score(dist) kps= IMP.core.KClosePairsPairScore(ps, tbr, 1.0) mpr= IMP.core.PairRestraint(kps, IMP.ParticlePair(p1,p2)) mpr.set_name(name) model.add_restraint(mpr) model.set_maximum_score(mpr, error_bound) def add_fret_restraint(ha, protein_a, residues_a, hb, protein_b, residues_b, r_value, reference=[]): if (residues_a=="C"): s1=IMP.atom.Selection(hierarchies=ha, molecule=protein_a, terminus=IMP.atom.Selection.C) elif (residues_a=="N"): s1=IMP.atom.Selection(hierarchies=ha, molecule=protein_a, terminus=IMP.atom.Selection.N) else: s1=IMP.atom.Selection(hierarchies=ha, molecule=protein_a, residue_index=residues_a) if (residues_b=="C"): s2=IMP.atom.Selection(hierarchies=hb, molecule=protein_b, terminus=IMP.atom.Selection.C) elif (residues_b=="N"): s2=IMP.atom.Selection(hierarchies=hb, molecule=protein_b, terminus=IMP.atom.Selection.N) else: s2=IMP.atom.Selection(hierarchies=hb, molecule=protein_b, residue_index=residues_b) p1=s1.get_selected_particles() p2=s2.get_selected_particles() #ra=IMP.core.XYZR(p1[0]).get_radius() #rb=IMP.core.XYZR(p2[0]).get_radius() range=get_range_from_fret_value(r_value) #range=(range0[0]-(ra+rb),range0[1]) if ( protein_a != protein_b): lp1=IMP.container.ListSingletonContainer(p1) lp2=IMP.container.ListSingletonContainer(p2) pc=IMP.container.AllBipartitePairContainer(lp1,lp2) elif (protein_a == protein_b and residues_a == residues_b): lp=IMP.container.ListSingletonContainer(p1) pc=IMP.container.AllPairContainer(lp) #f=IMP.atom.SameResiduePairFilter() #pc.add_pair_filter(f) elif (protein_a == protein_b and residues_a != residues_b): lp1=IMP.container.ListSingletonContainer(p1) lp2=IMP.container.ListSingletonContainer(p2) pc=IMP.container.AllBipartitePairContainer(lp1,lp2) #f=IMP.atom.SameResiduePairFilter() #pc.add_pair_filter(f) ps=get_pair_score(range) mpr=IMP.container.MinimumPairRestraint(ps,pc,1) model.add_restraint(mpr) model.set_maximum_score(mpr, error_bound) def get_range_from_fret_value(r_value): if (r_value >= 2.0) : return get_range_from_fret_class("High") if (r_value >= 1.5 and r_value < 2.0) : return get_range_from_fret_class("Mod") if (r_value >= 1.2 and r_value < 1.5) : return get_range_from_fret_class("Low") if (r_value >= 1.05 and r_value < 1.2) : return get_range_from_fret_class("Lowest") if (r_value < 1.05) : return get_range_from_fret_class("None") # take into account the size of the sphere def get_range_from_fret_class(r_class): if (r_class=="High"): return (-500.0, 41.0) if (r_class=="Mod"): return (41.0, 55.5) if (r_class=="Low"): return (55.5, 66.0) if (r_class=="Lowest"): return (66.0, 70.0) if (r_class=="None"): return (70.0, 100000.0) def get_pair_score(dist): hw=IMP.core.HarmonicWell(dist, 1.0) ps= IMP.core.DistancePairScore(hw) return ps def get_sphere_pair_score(dist): hw=IMP.core.HarmonicWell(dist, 1.0) ps= IMP.core.SphereDistancePairScore(hw) return ps def add_y2h_restraint(ha, protein_a, residues_a, hb, protein_b, residues_b): if (residues_a=="C"): s1=IMP.atom.Selection(hierarchies=ha, molecule=protein_a, terminus=IMP.atom.Selection.C) elif (residues_a=="N"): s1=IMP.atom.Selection(hierarchies=ha, molecule=protein_a, terminus=IMP.atom.Selection.N) elif (residues_a=="ALL"): s1=IMP.atom.Selection(hierarchies=ha, molecule=protein_a) else: s1=IMP.atom.Selection(hierarchies=ha, molecule=protein_a, residue_indexes=[residues_a]) if (residues_b=="C"): s2=IMP.atom.Selection(hierarchies=hb, molecule=protein_b, terminus=IMP.atom.Selection.C) elif (residues_b=="N"): s2=IMP.atom.Selection(hierarchies=hb, molecule=protein_b, terminus=IMP.atom.Selection.N) elif (residues_b=="ALL"): s2=IMP.atom.Selection(hierarchies=hb, molecule=protein_b) else: s2=IMP.atom.Selection(hierarchies=hb, molecule=protein_b, residue_indexes=[residues_b]) p1=s1.get_selected_particles() p2=s2.get_selected_particles() #print "y2h" #print [p.get_name() for p in p1] #print [p.get_name() for p in p2] range=(-500,0) ll= IMP.DEFAULT if ( protein_a != protein_b): lp1=IMP.container.ListSingletonContainer(p1, "C"+ protein_a+str(residues_a)) lp2=IMP.container.ListSingletonContainer(p2,"C"+ protein_a+str(residues_a)) pc=IMP.container.AllBipartitePairContainer(lp1,lp2, "ABPC"+protein_a+str(residues_a)+"-"\ +protein_b+str(residues_b)) elif (protein_a == protein_b and residues_a == residues_b): lp=IMP.container.ListSingletonContainer(p1, "C"+protein_a+str(residues_a)) pc=IMP.container.AllPairContainer(lp, "APC"+protein_a+str(residues_a)) #f=IMP.atom.SameResiduePairFilter() #ll=IMP.VERBOSE #pc.add_pair_filter(f) elif (protein_a == protein_b and residues_a != residues_b): lp1=IMP.container.ListSingletonContainer(p1, "C"+protein_a+str(residues_a)) lp2=IMP.container.ListSingletonContainer(p2, "C"+protein_b+str(residues_b)) pc=IMP.container.AllBipartitePairContainer(lp1,lp2, "ABPC"+protein_a+str(residues_a)+"-"\ +protein_b+str(residues_b)) #f=IMP.atom.SameResiduePairFilter() #pc.add_pair_filter(f) ps=get_sphere_pair_score(range) mpr=IMP.container.MinimumPairRestraint(ps,pc,1) mpr.set_log_level(ll) mpr.set_name("yh2:"+protein_a+str(residues_a)+": "\ +protein_b+str(residues_b)) model.add_restraint(mpr) model.set_maximum_score(mpr, error_bound) def add_symmetry_restraint(h): translations=[] transformations=[] if(cell_type!="square"): translations.append(IMP.algebra.Vector3D( 0.0, 0.0, 0.0 )) translations.append(IMP.algebra.Vector3D( 0.0, side*sqrt3, 0.0 )) translations.append(IMP.algebra.Vector3D( 1.5*side, side*sqrt3/2.0, 0.0 )) translations.append(IMP.algebra.Vector3D( 1.5*side, -side*sqrt3/2.0, 0.0 )) translations.append(IMP.algebra.Vector3D( 0.0, -side*sqrt3, 0.0 )) translations.append(IMP.algebra.Vector3D( -1.5*side, -side*sqrt3/2.0, 0.0 )) translations.append(IMP.algebra.Vector3D( -1.5*side, side*sqrt3/2.0, 0.0 )) num_rotations=3*num_cells/21 angle=360.0/num_rotations else: translations.append(IMP.algebra.Vector3D( 0.0, 0.0, 0.0 )) translations.append(IMP.algebra.Vector3D( 0.0, side, 0.0 )) translations.append(IMP.algebra.Vector3D( 0.0, -side, 0.0 )) translations.append(IMP.algebra.Vector3D( side, 0.0, 0.0 )) translations.append(IMP.algebra.Vector3D( side, side, 0.0 )) translations.append(IMP.algebra.Vector3D( side, -side, 0.0 )) translations.append(IMP.algebra.Vector3D( -side, 0.0, 0.0 )) translations.append(IMP.algebra.Vector3D( -side, side, 0.0 )) translations.append(IMP.algebra.Vector3D( -side, -side, 0.0 )) num_rotations=1 angle=0.0 # assembling translations and rotations (for rhombus cell) for tr in translations: for i in range(0,num_rotations): rot=IMP.algebra.get_rotation_about_axis(IMP.algebra.Vector3D(0,0,1), float(i)*math.radians(angle)) transformations.append(IMP.algebra.Transformation3D(rot,tr)) # applying restraints ps0=h[0].get_leaves() for i in range(1,num_cells): sm= IMP.core.TransformationSymmetry(transformations[i]) for j,ps in enumerate(h[i].get_leaves()): IMP.core.Reference.setup_particle(ps, ps0[j]) lc= IMP.container.ListSingletonContainer(h[i].get_leaves(), "copy "+str(i)) c= IMP.container.SingletonsConstraint(sm, None, lc, "Symmetry"+str(i)) model.add_score_state(c) def add_equivalence_filters(all, h): f=IMP.domino.EquivalenceSubsetFilterTable() eq_ps=[] for name in all: ps=[] for mol in h.get_children(): if (mol.get_name()==name): pp=mol.get_leaves() ps.append(pp) ncopies=len(ps) # number of copies of protein called "name" if( ncopies == 1): for i in range(0,len(ps[0])): eq_ps.append(ps[0][i]) if( ncopies > 1 ): for i in range(0,len(ps[0])): # cycle on number of particles per copy pps=[] for j in range(0,ncopies): # cycle of number of copies pps.append(ps[j][i]) f.add_set(pps) eq_ps.append(pps[0]) return (f,eq_ps) def setup_sampler(h0,cover0,h1,cover1,scale,filters): pst= IMP.domino.ParticleStatesTable() states0= IMP.domino.XYZStates(cover0) states1= IMP.domino.XYZStates(cover1) for ps in h0.get_leaves(): pst.set_particle_states(ps,states0) for ps in h1.get_leaves(): pst.set_particle_states(ps,states1) for r in model.get_restraints(): if(cell_type=="square"): r.set_maximum_score(scale**2) else: r.set_maximum_score(1.45*scale**2) opt=IMP.domino.OptimizeRestraints(model.get_root_restraint_set(),pst) ig=IMP.domino.get_interaction_graph(model.get_root_restraint_set(),pst) jt=IMP.domino.get_junction_tree(ig) mt=IMP.domino.get_merge_tree(jt) #IMP.show_graphviz(mt) lsft= IMP.domino.ListSubsetFilterTable(pst) rssft=IMP.domino.RestraintScoreSubsetFilterTable(model, pst) #esft= IMP.domino.ExclusionSubsetFilterTable(pst) filters.append(lsft) filters.append(rssft) filters[-1].set_use_caching(False) sampler= IMP.domino.DominoSampler(model, pst) sampler.set_subset_filter_tables(filters) sampler.set_log_level(IMP.SILENT) return (sampler, lsft, pst) def get_mapping(cover0, cover1): nn= IMP.algebra.NearestNeighbor3D(cover0) ret=[[] for c in cover0] for i,p in enumerate(cover1): nns=nn.get_nearest_neighbor(p) ret[nns].append(i) return ret def print_statistics(ps0,ass,subs,pst): nps=len(ps0) mean_dist=[] mean_dist2=[] for i in range(0,nps*(nps-1)/2): mean_dist.append(0.0) mean_dist2.append(0.0) for a in ass: IMP.domino.load_particle_states(subs, a, pst); for i in range(0,nps-1): d0=IMP.core.XYZ(ps0[i]) for j in range(i+1,nps): d1=IMP.core.XYZ(ps0[j]) idpair=i*nps-i*(i+1)/2+j-(i+1) dist=IMP.core.get_distance(d0,d1) mean_dist[idpair]+=dist mean_dist2[idpair]+=dist*dist # print results print "STATISTICS ENSEMBLE SOLUTIONS" for i in range(0,nps-1): name0=(IMP.atom.Hierarchy.decorate_particle(ps0[i])).get_name() for j in range(i+1,nps): name1=(IMP.atom.Hierarchy.decorate_particle(ps0[j])).get_name() idpair=i*nps-i*(i+1)/2+j-(i+1) ave=mean_dist[idpair]/len(ass) sig=math.sqrt(mean_dist2[idpair]/len(ass)-ave*ave) print name0,name1,ave,sig def display(grid1, h1,w): #for r in model.get_restraints(): # g=IMP.display.RestraintGeometry(r) # w.add_geometry(g) for point in grid1: g= IMP.display.SphereGeometry(IMP.algebra.Sphere3D(point,1.0)) g.set_name("grid_CP") w.add_geometry(g) for hs in h1: for ps in hs.get_leaves(): name=((IMP.atom.Hierarchy.decorate_particle(ps)).get_parent()).get_name() if(name=="Spc42p_n"): c=IMP.display.Color(219.0/255.0, 241.0/255.0, 250.0/255.0) if(name=="Spc29p"): c=IMP.display.Color(252.0/255.0, 210.0/255.0, 88.0/255.0) if(name=="Cmd1p"): c=IMP.display.Color(255.0/255.0, 255.0/255.0, 0.0/255.0) if(name=="Spc110p_c"): c=IMP.display.Color(255.0/255.0, 0.0/255.0, 0.0/255.0) g= IMP.display.XYZRGeometry(ps) g.set_name(name) g.set_color(c) w.add_geometry(g) for r in h1[0].get_model().get_restraints(): g= IMP.display.RestraintGeometry(r) w.add_geometry(g) # for hs in h2: # for ps in hs.get_leaves(): # name=((IMP.atom.Hierarchy.decorate_particle(ps)).get_parent()).get_name() # if(name=="Spc42p_c"): c=IMP.display.Color(0.0/255.0, 0.0/255.0, 255.0/255.0) # if(name=="Cnm67p_c"): c=IMP.display.Color(0.0/255.0, 210.0/255.0, 155.0/255.0) # g= IMP.display.XYZRGeometry(ps) # g.set_name(name) # g.set_color(c) # w.add_geometry(g)