pKa/pdb2pka_interface.py
def pdb2pka_sugelm(self):
        """Explore all possible mutations and calculate a phimap for each using pdb2pka (APBS)"""
        import Protool
        P=Protool.structureIO()
        P.readpdb(self.pdbfile)
        P.RemoveALT()
        #import Protool.mutate
        #MUT=Protool.mutate.Mutate(P)
        #
        # Construct arrays
        #
        import pKD_dict
        self.data=pKD_dict.pKD_dict()
        self.atom_data=pKD_dict.pKD_dict()
        #
        # Create dir for mutant PDB files
        #
        import os
        mutdir=os.path.join(self.topdir,self.pdbfile+'.pdbs')
        if not os.path.isdir(mutdir):
            os.mkdir(mutdir)
        #
        # Loop over all residues
        #
        residues=P.residues.keys()
        residues.sort()
        for residue in residues:
            orgres=P.resname(residue)
            print 'Calculating for %s %s' %(residue,P.resname(residue))
            #
            # If neutral mutate to Asp, Glu, Lys, Arg, His
            #
            targets=[]
            for res in ['ARG','LYS','HIS','ASP','GLU']:
                if P.resname(residue)!=res:
                    targets.append(res)
            #if orgres=='GLU':
            #    targets.append('GLN')
            #elif orgres=='ASP':
            #    targets.append('ASN')
            #elif orgres=='HIS':
            #    targets.append('PHE')
            #elif orgres=='ARG' or P.resname(residue)=='LYS':
            #    targets.append('MET')
            #
            # Target identified. Now model each
            #
            for target in targets:
                import pKD_tools
                resid=pKD_tools.get_resid_from_res(residue)
                orgres=P.resname(residue)
                filename=os.path.join(mutdir,'%s:%s:%s.pdb' %(residue,orgres,target))
                mutation='%s:%s:%s' %(residue,orgres,target)
                if not os.path.isfile(filename):
                    import Design_pKa_help
                    Design_pKa_help.make_mutation(self.pdbfile,mutation)
                NP=Protool.structureIO()
                NP.readpdb(filename)
                NP.writepdb(filename,TER=None)
                #
                # Calculate the interaction energies
                #
                protein,routines,forcefield,apbs_setup,lig_titgrps = pdb2pka.pre_init(pdbfilename=filename,
                                                                                      ff='parse',
                                                                                      ligand=None,
                                                                                      verbose=1)
                mypkaRoutines = pdb2pka.pKaRoutines(protein, routines, forcefield,apbs_setup)
                #
                # Find our group
                #
                sp=residue.split(':')
                chainid=sp[0]
                resnum=int(sp[1])
                mypkaRoutines.findTitratableGroups()
                this_pKa=None
                for pKa in mypkaRoutines.pKas:
                    print pKa.residue.resSeq,resnum
                    print pKa.residue.chainID,chainid
                    print pKa.residue.name,target
                    print pKa.pKaGroup.name,target
                    print '--------------'
                    print 'ChainID',pKa.residue.chainID
                    if pKa.residue.resSeq==resnum and pKa.residue.chainID==chainid and pKa.residue.name==target and pKa.pKaGroup.name==target:
                        #print 'Found group',pKa.residue.resSeq,pKa.pKaGroup.name
                        this_pKa=pKa
                        break
                if not this_pKa:
                    raise Exception,'Could not find inserted titratable group'
                mypkaRoutines.get_interaction_energies_setup(this_pKa,mode='pKD')
                matrix=mypkaRoutines.matrix
                #
                # Dig the interaction energies out of the pdb2pka array
                #
                for titration1 in matrix[this_pKa].keys():
                    for state1 in matrix[this_pKa][titration1].keys():
                        grp_sub=matrix[this_pKa][titration1][state1]
                        if mypkaRoutines.is_charged(this_pKa,titration1,state1):
                            for pKa2 in grp_sub.keys(): 
                                import string
                                chainID2=pKa.residue.chainID
                                resid2='%s:%s' %(chainID2,string.zfill(pKa2.residue.resSeq,4))
                                for titration2 in grp_sub[pKa2].keys():
                                    for state2 in grp_sub[pKa2][titration2].keys():
                                        if mypkaRoutines.is_charged(pKa2,titration2,state2):
                                            #
                                            # Both states are charged, so now we can pull the
                                            # interaction energies out
                                            #
                                            if not self.data.has_key(mutation):
                                                self.data[mutation]={}
                                            self.data[mutation][resid2]=grp_sub[pKa2][titration2][state2]
                                            #
                                            # Get the potentials at all atoms too
                                            #
                                            all_pots=mypkaRoutines.all_potentials[this_pKa][titration1][state1]
                                            sub_all_pots=all_pots[pKa2][titration2][state2]
                                            for atom in sub_all_pots.keys():
                                                resid=mutation
                                                import pKD_tools
                                                resid2=pKD_tools.get_resid_from_res(atom)
                                                atomname=atom.split(':')[-1] #atom.name
                                                if atomname[0]=='H' or atomname in ['N','C','O']:
                                                    continue # Skip all H atoms and all non-CA backbone atoms to save memory
                                                if not self.atom_data.has_key(resid):
                                                    self.atom_data[resid]={}
                                                if not self.atom_data[resid].has_key(resid2):
                                                    self.atom_data[resid][resid2]={}
                                                self.atom_data[resid][resid2][atomname]=abs(sub_all_pots[atom])
        return self.data,self.atom_data
Explore all possible mutations and calculate a phimap for each using pdb2pka (APBS)
Explore all possible mutations and calculate a phimap for each using pdb2pka (APBS)
Explore all possible mutations and calculate a phimap for each using pdb2pka ( APBS )
def pdb2pka_sugelm(self):
        import Protool
        P=Protool.structureIO()
        P.readpdb(self.pdbfile)
        P.RemoveALT()
         Construct arrays
        import pKD_dict
        self.data=pKD_dict.pKD_dict()
        self.atom_data=pKD_dict.pKD_dict()
         Create dir for mutant PDB files
        import os
        mutdir=os.path.join(self.topdir,self.pdbfile+'.pdbs')
        if not os.path.isdir(mutdir):
            os.mkdir(mutdir)
         Loop over all residues
        residues=P.residues.keys()
        residues.sort()
        for residue in residues:
            orgres=P.resname(residue)
            print 'Calculating for %s %s' %(residue,P.resname(residue))
             If neutral mutate to Asp, Glu, Lys, Arg, His
            targets=[]
            for res in ['ARG','LYS','HIS','ASP','GLU']:
                if P.resname(residue)!=res:
                    targets.append(res)
            if orgres=='GLU':
                targets.append('GLN')
            elif orgres=='ASP':
                targets.append('ASN')
            elif orgres=='HIS':
                targets.append('PHE')
            elif orgres=='ARG' or P.resname(residue)=='LYS':
                targets.append('MET')
             Target identified. Now model each
            for target in targets:
                import pKD_tools
                resid=pKD_tools.get_resid_from_res(residue)
                orgres=P.resname(residue)
                filename=os.path.join(mutdir,'%s:%s:%s.pdb' %(residue,orgres,target))
                mutation='%s:%s:%s' %(residue,orgres,target)
                if not os.path.isfile(filename):
                    import Design_pKa_help
                    Design_pKa_help.make_mutation(self.pdbfile,mutation)
                NP=Protool.structureIO()
                NP.readpdb(filename)
                NP.writepdb(filename,TER=None)
                 Calculate the interaction energies
                protein,routines,forcefield,apbs_setup,lig_titgrps = pdb2pka.pre_init(pdbfilename=filename,
                                                                                      ff='parse',
                                                                                      ligand=None,
                                                                                      verbose=1)
                mypkaRoutines = pdb2pka.pKaRoutines(protein, routines, forcefield,apbs_setup)
                 Find our group
                sp=residue.split(':')
                chainid=sp[0]
                resnum=int(sp[1])
                mypkaRoutines.findTitratableGroups()
                this_pKa=None
                for pKa in mypkaRoutines.pKas:
                    print pKa.residue.resSeq,resnum
                    print pKa.residue.chainID,chainid
                    print pKa.residue.name,target
                    print pKa.pKaGroup.name,target
                    print '--------------'
                    print 'ChainID',pKa.residue.chainID
                    if pKa.residue.resSeq==resnum and pKa.residue.chainID==chainid and pKa.residue.name==target and pKa.pKaGroup.name==target:
                        print 'Found group',pKa.residue.resSeq,pKa.pKaGroup.name
                        this_pKa=pKa
                        break
                if not this_pKa:
                    raise Exception,'Could not find inserted titratable group'
                mypkaRoutines.get_interaction_energies_setup(this_pKa,mode='pKD')
                matrix=mypkaRoutines.matrix
                 Dig the interaction energies out of the pdb2pka array
                for titration1 in matrix[this_pKa].keys():
                    for state1 in matrix[this_pKa][titration1].keys():
                        grp_sub=matrix[this_pKa][titration1][state1]
                        if mypkaRoutines.is_charged(this_pKa,titration1,state1):
                            for pKa2 in grp_sub.keys(): 
                                import string
                                chainID2=pKa.residue.chainID
                                resid2='%s:%s' %(chainID2,string.zfill(pKa2.residue.resSeq,4))
                                for titration2 in grp_sub[pKa2].keys():
                                    for state2 in grp_sub[pKa2][titration2].keys():
                                        if mypkaRoutines.is_charged(pKa2,titration2,state2):
                                             Both states are charged, so now we can pull the
                                             interaction energies out
                                            if not self.data.has_key(mutation):
                                                self.data[mutation]={}
                                            self.data[mutation][resid2]=grp_sub[pKa2][titration2][state2]
                                             Get the potentials at all atoms too
                                            all_pots=mypkaRoutines.all_potentials[this_pKa][titration1][state1]
                                            sub_all_pots=all_pots[pKa2][titration2][state2]
                                            for atom in sub_all_pots.keys():
                                                resid=mutation
                                                import pKD_tools
                                                resid2=pKD_tools.get_resid_from_res(atom)
                                                atomname=atom.split(':')[-1] atom.name
                                                if atomname[0]=='H' or atomname in ['N','C','O']:
                                                    continue  Skip all H atoms and all non-CA backbone atoms to save memory
                                                if not self.atom_data.has_key(resid):
                                                    self.atom_data[resid]={}
                                                if not self.atom_data[resid].has_key(resid2):
                                                    self.atom_data[resid][resid2]={}
                                                self.atom_data[resid][resid2][atomname]=abs(sub_all_pots[atom])
        return self.data,self.atom_data
def pdb2pka_sugelm ( self ) : import Protool P = Protool . structureIO ( ) P . readpdb ( self . pdbfile ) P . RemoveALT ( ) import pKD_dict self . data = pKD_dict . pKD_dict ( ) self . atom_data = pKD_dict . pKD_dict ( ) import os mutdir = os . path . join ( self . topdir , self . pdbfile + '.pdbs' ) if not os . path . isdir ( mutdir ) : os . mkdir ( mutdir ) residues = P . residues . keys ( ) residues . sort ( ) for residue in residues : orgres = P . resname ( residue ) print 'Calculating for %s %s' % ( residue , P . resname ( residue ) ) targets = [ ] for res in [ 'ARG' , 'LYS' , 'HIS' , 'ASP' , 'GLU' ] : if P . resname ( residue ) != res : targets . append ( res ) for target in targets : import pKD_tools resid = pKD_tools . get_resid_from_res ( residue ) orgres = P . resname ( residue ) filename = os . path . join ( mutdir , '%s:%s:%s.pdb' % ( residue , orgres , target ) ) mutation = '%s:%s:%s' % ( residue , orgres , target ) if not os . path . isfile ( filename ) : import Design_pKa_help Design_pKa_help . make_mutation ( self . pdbfile , mutation ) NP = Protool . structureIO ( ) NP . readpdb ( filename ) NP . writepdb ( filename , TER = None ) protein , routines , forcefield , apbs_setup , lig_titgrps = pdb2pka . pre_init ( pdbfilename = filename , ff = 'parse' , ligand = None , verbose = 1 ) mypkaRoutines = pdb2pka . pKaRoutines ( protein , routines , forcefield , apbs_setup ) sp = residue . split ( ':' ) chainid = sp [ 0 ] resnum = int ( sp [ 1 ] ) mypkaRoutines . findTitratableGroups ( ) this_pKa = None for pKa in mypkaRoutines . pKas : print pKa . residue . resSeq , resnum print pKa . residue . chainID , chainid print pKa . residue . name , target print pKa . pKaGroup . name , target print '--------------' print 'ChainID' , pKa . residue . chainID if pKa . residue . resSeq == resnum and pKa . residue . chainID == chainid and pKa . residue . name == target and pKa . pKaGroup . name == target : this_pKa = pKa break if not this_pKa : raise Exception , 'Could not find inserted titratable group' mypkaRoutines . get_interaction_energies_setup ( this_pKa , mode = 'pKD' ) matrix = mypkaRoutines . matrix for titration1 in matrix [ this_pKa ] . keys ( ) : for state1 in matrix [ this_pKa ] [ titration1 ] . keys ( ) : grp_sub = matrix [ this_pKa ] [ titration1 ] [ state1 ] if mypkaRoutines . is_charged ( this_pKa , titration1 , state1 ) : for pKa2 in grp_sub . keys ( ) : import string chainID2 = pKa . residue . chainID resid2 = '%s:%s' % ( chainID2 , string . zfill ( pKa2 . residue . resSeq , 4 ) ) for titration2 in grp_sub [ pKa2 ] . keys ( ) : for state2 in grp_sub [ pKa2 ] [ titration2 ] . keys ( ) : if mypkaRoutines . is_charged ( pKa2 , titration2 , state2 ) : if not self . data . has_key ( mutation ) : self . data [ mutation ] = { } self . data [ mutation ] [ resid2 ] = grp_sub [ pKa2 ] [ titration2 ] [ state2 ] all_pots = mypkaRoutines . all_potentials [ this_pKa ] [ titration1 ] [ state1 ] sub_all_pots = all_pots [ pKa2 ] [ titration2 ] [ state2 ] for atom in sub_all_pots . keys ( ) : resid = mutation import pKD_tools resid2 = pKD_tools . get_resid_from_res ( atom ) atomname = atom . split ( ':' ) [ - 1 ] if atomname [ 0 ] == 'H' or atomname in [ 'N' , 'C' , 'O' ] : continue if not self . atom_data . has_key ( resid ) : self . atom_data [ resid ] = { } if not self . atom_data [ resid ] . has_key ( resid2 ) : self . atom_data [ resid ] [ resid2 ] = { } self . atom_data [ resid ] [ resid2 ] [ atomname ] = abs ( sub_all_pots [ atom ] ) return self . data , self . atom_data
Explore all possible mutations and calculate a phimap for each using pdb2pka (APBS)
Explore all possible mutations and calculate a phimap for each using pdb2pka ( APBS )
"""Explore all possible mutations and calculate a phimap for each using pdb2pka (APBS)""" #import Protool.mutate #MUT=Protool.mutate.Mutate(P) # # Construct arrays # # # Create dir for mutant PDB files # # # Loop over all residues # # # If neutral mutate to Asp, Glu, Lys, Arg, His # #if orgres=='GLU': #    targets.append('GLN') #elif orgres=='ASP': #    targets.append('ASN') #elif orgres=='HIS': #    targets.append('PHE') #elif orgres=='ARG' or P.resname(residue)=='LYS': #    targets.append('MET') # # Target identified. Now model each # # # Calculate the interaction energies # # # Find our group # #print 'Found group',pKa.residue.resSeq,pKa.pKaGroup.name # # Dig the interaction energies out of the pdb2pka array # # # Both states are charged, so now we can pull the # interaction energies out # # # Get the potentials at all atoms too # #atom.name # Skip all H atoms and all non-CA backbone atoms to save memory
self