Python restraints.EDrestraint类代码示例

def calcCC(crdlist,pts,mapfile,esmin,esmax):
    crdlist1 = VecVecFloat(crdlist)
    pis = VecInt( range(crdlist1.size()) )
    from restraints import EDrestraint
    er  = EDrestraint.makeEDrestraintFromMap(pis, "Xranker", mapfile, esmin, esmax, esmax)
    xcorr = EDrestraint.findBadfit1(mtzfn, f1label, pts, ptinds)
    print "xcor",xcorr,resids[ri]
    return xcorr

def makeMaps(mtzfile,pdbfile,sg):
    a,b,c,alpha,beta,gamma = getCRYST(pdbfile)
    cif2mtz(mtzfile, "base.mtz", a, b, c, alpha, beta, gamma, sg)
    sfallA(pdbfile,"base.mtz" , "%s.phased.mtz" %pdbfile)
    phasedmtz = "%s.phased.mtz"%pdbfile
    
    
    ## 0 : 2F1-F2 , 1: F1 , 2: F2
    from restraints import EDrestraint
    esmin, esmax, esmean, rcmult = 0.0001, 5.0, 0.05, 10
    EDrestraint.makeMap(phasedmtz,"FP" ,"FC", "PHIC", 1 , esmin, esmax, esmax)

 def score(s, pdbfile, mtzfn, f1label="FP", f2label="FC", philabel="PHIC", maptype="2F1-F2", aasc="mcsc") :
     assert aasc in ["pept","mc","sc","mcsc"]
     prot = protein(pdbfile, read_hydrogens=0, read_waters=1, read_hets=1)
     res, resids, resnums, resns, chids, inscodes, pts = prot2res.readProtRes(prot)
     pts = VecVecFloat(pts)
     retkeys = []
     badlines = []
     allkeys = list(res.keys()) ; allkeys.sort()
     for ai in range(len(allkeys)) :
         ri = allkeys[ai]
         key = resids[ri]
         if not isAAres(resns[ri]) : continue
         if aasc=="sc" and resns[ri] in ["GLY","ALA"] : continue
         if aasc=="sc" and isAAres(resns[ri]) :
             ptinds = []
             for k,v in res[ri].items() :
                 if not k in [" N  "," CA "," C  "," O  "," CB "] : ptinds.append(v)
             ptinds = VecInt(ptinds)
         elif aasc=="pept" and isAAres(resns[ri]) and ai+1 < len(allkeys) and isAAres(resns[allkeys[ai+1]]) :
             ptinds = []
             for k,v in res[ri].items() :
                 if k in [" C  "," O  "] : ptinds.append(v)
             for k,v in res[ allkeys[ai+1] ].items() :
                 if k in [" N  ",] : ptinds.append(v)
             ptinds = VecInt(ptinds)
         elif aasc=="mc" and isAAres(resns[ri]) :
             ptinds = []
             for k,v in res[ri].items() :
                 if k in [" N  "," CA "," C  "," O  "," CB "] : ptinds.append(v)
             ptinds = VecInt(ptinds)
         elif aasc=="ca" and isAAres(resns[ri]) :
             ptinds = []
             for k,v in res[ri].items() :
                 if k in [" CA "] : ptinds.append(v)
             ptinds = VecInt(ptinds)                
         else :
             ptinds = VecInt( res[ri].values() )
         if re.compile("\.map$").search(mtzfn) and re.compile("\.map$").search(f1label) :
             xcorr = EDrestraint.findBadfit1(mtzfn, f1label, pts, ptinds)
             print "xcor",xcorr,resids[ri]
             
         else :
             xcorr = EDrestraint.findBadfit(mtzfn, f1label, f2label, philabel, pts, ptinds)
         #print "XcorrZ [%s]"%key, xcorr
         if not s.scores.has_key(key) : s.scores[key] = [] ; s.keyorder.append(key)
         s.scores[key].append(xcorr)
         if xcorr < s.cutoff :
             retkeys.append(key) ;
             bl = "[%s]" % key
             badlines.append(bl)
     print len(badlines), "BADRESIDS", aasc, len(s.scores.values()[0]), s.cutoff, "------------------------------------" ;
     for bl in badlines : print bl
     s.gplot() 
     return retkeys

 def score(s, crdlist, pis=None) :
     crdlist1 = VecVecFloat(crdlist)
     #print "-----------------------------------"
     #for ci in range(crdlist1.size()) :
     #    print crdlist1[ci][0] , crdlist1[ci][0] , crdlist1[ci][0]
     if not pis : pis = VecInt( range(crdlist1.size()) )
     else : pis = VecInt( range(pis) )
     if re.compile("\.map$").search(s.phasedmtz) :
         er = EDrestraint.makeEDrestraintFromMap(pis, "Xranker", s.phasedmtz, s.esmin, s.esmax, s.esmax)
     else :
         er = EDrestraint.makeEDrestraintFromMTZ(pis, "Xranker", s.phasedmtz, s.fplabel, s.fclabel, s.philabel, s.maptype, s.esmin, s.esmax, s.esmax)
     return -1*er.scoreAround(crdlist1)

def main(pdbfile, mtzfn, f1label, f2label, philabel, maptype, aasc=None) :
    mt = {"2F1-F2":0, "F1":1}
    prot = protein(pdbfile, read_hydrogens=0, read_waters=1, read_hets=1)
    res, resids, resnums, resns, chids, inscodes, pts = prot2res.readProtRes(prot)
    pts = VecVecFloat(pts)
    Xcorrs = {}
#    print len(res.keys()), len(set(res.keys()))
    allkeys = list(res.keys()) ; allkeys.sort()
    for ri in allkeys :
        key = resids[ri]
        if aasc and (not isAAres(resns[ri]) or resns[ri] in ["GLY","ALA"]) : continue
        if aasc :
            ptinds = []
            for k,v in res[ri].items() :
                if not k in [" N  "," CA "," C  "," O  "," CB "] : ptinds.append(v)
            ptinds = VecInt(ptinds)
        else :
            ptinds = VecInt( res[ri].values() )

        xcorr = EDrestraint.findBadfit(mtzfn, f1label, f2label, philabel, pts, ptinds)
        Xcorrs[ key ] = xcorr
#        print ri, "RESID", key        print "XcorrZ [%s]"%key, xcorr
    #mean, stdev = findMeanStdev(Xcorrs.values())
    #for key,val in Xcorrs.items() : print "XcorrZ", key, val, (val-mean) / stdev
    return Xcorrs

    def score_section(s, pdbfile, mtzfn, f1label="FP", f2label="FC", philabel="PHIC", maptype="2F1-F2", aasc="mcsc",ptinds=[]) :
        assert aasc in ["pept","mc","sc","mcsc"]
        prot = protein(pdbfile, read_hydrogens=0, read_waters=1, read_hets=1)
        res, resids, resnums, resns, chids, inscodes, pts = prot2res.readProtRes(prot)
        pts = VecVecFloat(pts)
        xcorr = EDrestraint.findBadfit1(mtzfn, f1label, pts, ptinds)
        print "XcorrZ = ", xcorr

        return []

def getTime(sf,pdb):
    xrayRestGen = []
    prot = protein(pdb, read_hydrogens=0, read_waters=0, read_hets=0)
    res, resids, resnums, resns, chids, inscodes, pts = prot2res.readProtRes(prot)
    esmin, esmax, esmean, rcmult, xscoreCutoff = .000, 5., .0, 5, 0.9
    pis = VecInt(resids.keys())
    er = EDrestraint.makeEDrestraintFromMap(pis, "Xranker", phasedmtz, esmin, esmax, smax)
    score = er.scoreAll(VecVecFloat(pts))
    print score

def createXEnvelopeRestraints(blist, mapfile, aiSC) :
    from restraints import EDrestraint
    erlist = []
    for b in blist :
        bop = b.getOP()
        edrIP = []
        for i in range(bop.size()) :
            if aiSC[bop[i]] != -1 : edrIP.append(bop[i])
        rstr = "XEnvelopeRestraint on op of %s" % b.name()
        esmin, esmax, esmean, rcmult = .000, 5., .0, 5
        erlist.append( EDrestraint.makeEDrestraintFromMap(VecInt(edrIP), rstr, mapfile,esmin, esmax, esmean) )
    return erlist

 def generateMTZ(s, blist, aiSC) :
     xrlist, optional = [], []
     for b in blist :
         bop = b.getOP()
         edrIP = []
         for i in range(bop.size()) :
             if aiSC[bop[i]] != -1 : edrIP.append(bop[i])
         if len(edrIP) > 0 :
             xrlist.append( EDrestraint.makeEDrestraintFromMTZ(VecInt(edrIP), "EDrestraint on o/p of %s" % b.name(), \
                 s.mtzfn, s.folabel, s.fclabel, s.philabel, s.maptype, s.min, s.max, s.mean) )
             if not "LigandBuilder" in b.name() : optional.append(len(xrlist)-1)
     return xrlist, optional

 def generateMAP(s, blist, aiSC) :
     print s.mapfn, s.min, s.max, s.mean
     erlist = []
     for b in blist :
         bop = b.getOP()
         edrIP = []
         for i in range(bop.size()) :
             if aiSC[bop[i]] != -1 : edrIP.append(bop[i])
         rstr = "XEnvelopeRestraint on op of %s" % b.name()
         print s.mapfn, s.min, s.max, s.mean
         erlist.append( EDrestraint.makeEDrestraintFromMap(VecInt(edrIP), rstr, s.mapfn, s.min, s.max, s.mean) )
     return erlist

def createXEnvelopeRestraintsCA(blist, mapfile, aiSC) :
    from restraints import EDrestraint
    erlist = []
    for b in blist :
        bop = b.getOP()
        edrIP = []
        if "Pept" in b.name() or "anchor" in b.name(): 
            for i in range(bop.size()) :
                print bop[i]
                #if aiSC[bop[i]] != -1 :
                edrIP.append(bop[i])
                rstr = "XEnvelopeRestraint on op of %s" % b.name()
            esmin, esmax, esmean, rcmult = .000, 5., .0, 5
            erlist.append( EDrestraint.makeEDrestraintFromMap(VecInt(edrIP), rstr, mapfile,esmin, esmax, esmean) )
        else :
            print "DDD", b.name() 
            continue

    return erlist

def createXEnvelopeRestraints2(blist, mapfile, aiSC,ncDummies) :
    from restraints import EDrestraint
    erlist = []
    
    print "aiSC",aiSC
    print "NCD",ncDummies
    for b in blist :
        bop = b.getOP()
        edrIP = []
        rstr = "XEnvelopeRestraint on op of %s" % b.name()
        print rstr
        for i in range(bop.size()) :
            print i , bop[i]

            if aiSC[bop[i]] != -1 :
                edrIP.append(bop[i])

        esmin, esmax, esmean, rcmult = .000, 5., .0, 5
        erlist.append( EDrestraint.makeEDrestraintFromMap(VecInt(edrIP), rstr, mapfile,esmin, esmax, esmean) )
    return erlist

def getTime(sf,pdb,sg):
    a,b,c,alpha,beta,gamma = getCRYST(pdb)
    cif2mtz(sf, "base.mtz", a, b, c, alpha, beta, gamma, sg)
    uniqueify("base.mtz", "rfree.mtz")
    sfall(pdb, "rfree.mtz", "phased.mtz")
    xrayRestGen = []
    prot = protein(pdb, read_hydrogens=0, read_waters=0, read_hets=0)
    res, resids, resnums, resns, chids, inscodes, pts = prot2res.readProtRes(prot)
    esmin, esmax, esmean, rcmult, xscoreCutoff = .000, 5., .0, 5, 0.9
    pis = VecInt(resids.keys())
    folabel="FP"
    fclabel="FC"
    philabel="PHIC"
    maptype="2F1-F2"

    if maptype == '2F1-F2' : maptype = 0
    elif maptype == 'F1' : maptype = 1
    else : print "unknown maptype ", maptype ; sys.exit(1)
    er = EDrestraint.makeEDrestraintFromMTZ(pis, "EDrestraint", "phased.mtz", folabel, fclabel, philabel, maptype, esmin, esmax, esmean) 
    score = er.scoreAll(VecVecFloat(pts))
    print score
    sys.exit()

#.........这里部分代码省略.........
    l =  "--gamma = %s "%str(options.gamma) ;      print >> fp, l 
    l =  "--sg = %s "%str(options.sg) ;      print >> fp, l 
    l =  "--resolution = %s "%str(options.resolution) ;      print >> fp, l 
    l =  "--ca-restraint-radius = %s "%str(options.caRad) ;      print >> fp, l 
    l =  "--sc-centroid-restraint-radius = %s "%str(options.scRad) ;      print >> fp, l 
    l =  "--sidechain-vdw-reduction = %s "%str(options.scReduction) ;      print >> fp, l 
    l =  "--population-size = %s "%str(options.popsize) ;      print >> fp, l 
    l =  "--verbose = %s "%str(options.verbose) ;      print >> fp, l 
    l =  "--backtrack = %s "%str(options.backtrack) ;      print >> fp, l 
    l =  "--randomize = %s "%str(options.randomize) ;      print >> fp, l 
    l =  "--loopres = %s "%str(options.loopres) ;      print >> fp, l 
    l =  "--framework = %s "%str(options.framework) ;      print >> fp, l 
    l =  "--framework-ca-threshold = %s "%str(options.fcarad) ;      print >> fp, l 
    l =  "--framework-sc-threshold = %s "%str(options.fscrad) ;      print >> fp, l 
    l =  "--rebuild-poor-regions-only = %s "%str(options.poorOnly) ;      print >> fp, l 
    l =  "--poor-fit-threshold = %s "%str(options.poorThreshold) ;      print >> fp, l 
    l =  "--mconly = %s "%str(options.mconly) ;      print >> fp, l 
    l =  "--sconly = %s "%str(options.sconly) ;      print >> fp, l 
    l =  "--default-mainchain-b-factor = %s "%str(options.mcBfac) ;      print >> fp, l 
    l =  "--default-sidechain-b-factor = %s "%str(options.scBfac) ;      print >> fp, l 
    l =  "--models-get-native-bfactors = %s "%str(options.nativeBfac) ;      print >> fp, l 
    l =  "--dontusesc = %s "%str(options.dontusesc) ;      print >> fp, l 
    l =  "--num-models = %s "%str(options.nmodels) ;      print >> fp, l 
    l =  "--f1label = %s "%str(options.f1label) ;      print >> fp, l 
    l =  "--f2label = %s "%str(options.f2label) ;      print >> fp, l 
    l =  "--philabel = %s "%str(options.philabel) ;      print >> fp, l 
    l =  "--maptype = %s "%str(options.maptype) ;      print >> fp, l 
    l =  "--cacaCutoff = %s "%str(options.cacaCutoff) ;      print >> fp, l 
    l =  "--make-ed-optional = %s "%str(options.edOpt) ;      print >> fp, l
#    l =  "--chids = %s "%str(options.chids) ;      print >> fp, l 
    fp.close()
    
    guidedsampling = None
    from restraints import EDrestraint ; EDrestraint.setPenalty(2.) ; EDrestraint.setScatRad(1.) ;
    esmin, esmax, esmean, rcmult, xscoreCutoff = .000, 5., .0, 5, options.poorThreshold
    multiPrepC = prepareChain.PrepareChain("SCL1.0")
    xscorer = XrayScorer(None, xscoreCutoff)

    
    modelIn =  options.pdbfile
    rtkmodel = options.pdbout # rappertk model to be generated in this cycle


    if options.mapfn1 == None and options.mtzfn == None:
        print "Warning no mtz or map input, no electron density restraints will be used"
    
    if options.poorOnly == 1 :
        if  options.mtzfn == None and options.mapfn1 == None and options.mapfn2 == None:
            print "No mapfn1,mapfn2 or mtzfn  set"
            sys.exit()
        if options.mapfn2 == None and options.mapfn1 !=None and options.mtzfn == None:
            print "Both mapfn1 and mapfn2 need to be  set"
            sys.exit()            

        if options.mapfn1 == None and options.mapfn2 !=None and options.mtzfn == None :
            print "Both mapfn1 and mapfn2 need to be  set"
            sys.exit()            


        if options.mtzfn != None :
            badmcsc = xscorer.score(modelIn, options.mtzfn, options.f1label, options.f2label, options.philabel, options.maptype, "mcsc")
            badmc = xscorer.score(modelIn, options.mtzfn, options.f1label, options.f2label, options.philabel, options.maptype, "mc")
            badpept = xscorer.score(modelIn, options.mtzfn, options.f1label, options.f2label, options.philabel, options.maptype,"pept")
            badscs = xscorer.score(modelIn, options.mtzfn, options.f1label, options.f2label, options.philabel, options.maptype,"sc")
            
        elif options.mapfn1 !=None and options.mapfn2 != None  :

def main() :

    import optparse ; parser = optparse.OptionParser()
    parser.add_option("--scratchdir", action='store', type='string', dest='scratchdir', help='to create all the files during refinement. it shdnt be already present.')
    parser.add_option("--pdb", action='store', type='string', dest='pdbfile', help='starting pdb containing a model of pdb-ligand complex')

    parser.add_option("--sf", action='store', type='string', dest='sf', help='structure factors file')
    parser.add_option("--a", action='store', type='float', dest='a', help='cell dimension a')
    parser.add_option("--b", action='store', type='float', dest='b', help='cell dimension b')
    parser.add_option("--c", action='store', type='float', dest='c', help='cell dimension c')
    parser.add_option("--alpha", action='store', type='float', dest='alpha', help='cell angle alpha')
    parser.add_option("--beta", action='store', type='float', dest='beta', help='cell angle beta')
    parser.add_option("--gamma", action='store', type='float', dest='gamma', help='cell angle gamma')
    parser.add_option("--sg", action='store', type='string', dest='sg', help='cell spacegroup, in CCP4 notation')
    parser.add_option("--resolution", action='store', type='float', dest='resolution', help='resolution of the data')

    parser.add_option("--ca-restraint-radius", action='store', type='float', dest='caRad', help='radius of spherical restraint on CA position', default=1)
    parser.add_option("--sc-centroid-restraint-radius", action='store', type='float', dest='scRad', help='radius of spherical restraint on sidechain centroid', default=2)
    parser.add_option("--sidechain-vdw-reduction", action='store', type='float', dest='scReduction', help='factor to reduce effective vdw dist in case of sidechains', default=1)
    parser.add_option("--population-size", action='store', type='int', dest='popsize', help='population size for PopulationStrategy', default=100)
    parser.add_option("--verbose", action='store', type='int', dest='verbose', help='0 means least verbosity etc.', default=0)
    parser.add_option("--backtrack", action='store', type='string', dest='backtrack', help='use backtracking version of PopulationStrategy. eg 4X5 will set backtrack numsteps and stepsize to 4,5 respectively. not used by default.', default=None)
    parser.add_option("--noRTK", action='store', type='int', dest='noRTK', help='dont rebuild bad-fits with rtk', default=0)
    parser.add_option("--randomize", action='store', type='int', dest='randomize', help='seed for randomizing', default=None)
    parser.add_option("--loopres", action='store', type='string', dest='loopres', help='filename containing resids for starting perturbation', default=None)
    parser.add_option("--framework", action='store', type='int', dest='framework', help='to be used in conjunction with loopres. it puts a 1/3 ca/sc pos restr on non-loopres and perturbs them too', default=None)


    (options, args) = parser.parse_args()



    import misc; misc.setVerbosity(options.verbose)
    randomize(options.randomize)
    if not os.path.isdir(options.scratchdir) : os.mkdir(options.scratchdir)
    shutil.copyfile(options.pdbfile, "%s/model0.pdb" % options.scratchdir)
    shutil.copyfile(options.sf, "%s/rfree.mtz" % options.scratchdir)
    os.chdir(options.scratchdir)

    #cif2mtz("strfactors.mtz", "base.mtz", options.a, options.b, options.c, options.alpha, options.beta, options.gamma, options.sg)
    #uniqueify("base.mtz", "rfree.mtz")


    esmin, esmax, esmean, rcmult, xscoreCutoff = .000, 5., .0, 5, 0.9
    guidedsampling = None
    multiPrepC = prepareChain.PrepareChain("SCL1.0")
    numRefCycles = 10 ; startCycle = 0

    from stump import getCRYST , getRESO
    if (options.a == None or options.b == None or options.c == None or options.alpha == None or options.beta == None or options.gamma == None) :
        
        print "Getting cell paramters from coordinate file....."
        options.a,options.b,options.c,options.alpha , options.beta , options.gamma,d1  = getCRYST(options.pdbfile)
        
        if (options.a == None or options.b == None or options.c == None or options.alpha== None or options.beta==None or options.gamma == None ):
            print "CRYST card cannot be read from coordinate file. Please input cell paramater a, b , c , alpha, beta , gamma = ",options.a , options.b , options.c , options.alpha , options.beta  , options.gamma 
            import sys ; sys.exit()

        print "Found a b c alpha beta gamma  ", options.a , options.b , options.c , options.alpha , options.beta  , options.gamma 
            
    if options.sg == None : 
        print "Getting space group from coordinate file....."
        d1,d2,d3,d4 , d5 , d6, options.sg  = getCRYST(options.pdbfile)
        if options.sg == None : 
            print "Please input space group " , options.sg ; import sys ; sys.exit()
        ss = ""
        for sg1 in options.sg:
            if sg1 in ["\n","\t","\s"]:
                continue
            else :
                ss = ss+sg1
        options.sg = ss
        print "Space Group ",options.sg

    if options.sg  in long2shortHM.keys():
        shortsg = long2shortHM[options.sg]
        options.sg = shortsg
    if options.sg not in sgtable.keys():
        print "Check --sg , Not recognised [%s][%d]"%( options.sg, len(options.sg))
        import sys ; sys.exit()

    if options.resolution == None : 
        print "Getting resolution limit from coordinate file........"
        options.resolution = getRESO(options.pdbfile)
        if (options.resolution == None):
            print "Please input resolution " , options.resolution
            import sys ; sys.exit()
                
        print "Resolution = [ " , options.resolution, " ] "


    from restraints import EDrestraint ; EDrestraint.setPenalty(2.) ; EDrestraint.setScatRad(1.) ;
    for cycle in range(startCycle, numRefCycles) :

        if cycle < 10 :
            xscoreCutoff = 0.8
        else :
            xscoreCutoff = 0.9

        xscorer = XrayScorer(None, xscoreCutoff) ; xrayRestGen = []
#.........这里部分代码省略.........