Python Protool类代码示例

    def update_elements(self):
        """Insert a new dropdown list for the element"""
        #
        # Get the group type
        #
        elements = None
        group_type = self.group_type_box.getcurselection()[0]
        import Protool

        if group_type == "Residues":
            P = Protool.structureIO()
            P.parsepdb(self.pdblines)
            residues = P.residues.keys()
            residues.sort()
            elements = []
            for res in residues:
                elements.append("%s %s" % (res, P.resname(res)))
        elif group_type == "Atoms":
            P = Protool.structureIO()
            P.parsepdb(self.pdblines)
            atoms = P.atoms.keys()
            for res in P.residues.keys():
                resname = P.resname(res)
                if self.AAdefs.has_key(resname):
                    defatoms = self.AAdefs[resname]["atoms"]
                    # print defatoms
                    for defatom, coord, dummy in defatoms:
                        atom_name = "%s:%s" % (res, defatom)
                        if not P.atoms.has_key(atom_name):
                            atoms.append(atom_name)
                            # print 'Adding',atom_name
            atoms.sort()
            elements = []
            for at in atoms:
                elements.append(at)
        elif group_type == "Titratable groups":
            P = Protool.structureIO()
            P.parsepdb(self.pdblines)
            P.get_titratable_groups()
            titgrps = P.titratable_groups.keys()
            titgrps.sort()
            elements = []
            for res in titgrps:
                for titgrp in P.titratable_groups[res]:
                    name = "%s %s" % (res, titgrp["name"])
                    elements.append(name)
        else:
            print "Unkown group type", group_type
        #
        # Make the new dropdown list
        #
        if elements:
            self.group_elements_box.setlist(elements)
        return

def get_net_charge(pdbfile,HIS):
    """Get the net charge within 20 A of the HIS"""
    import Protool
    X=Protool.structureIO()
    X.readpdb(pdbfile)
    close=[]
    HIS_ND1='%s:ND1' %HIS
    HIS_NE2='%s:NE2' %HIS
    for residue in X.residues.keys():
        for atom in X.residues[residue]:
            #print atom
            mdist=min(X.dist(HIS_ND1,atom),X.dist(HIS_NE2,atom))
            if mdist<50.0:
                close.append(residue)
                break
            elif mdist>355.0:
                break
    # Got all close residues, now count charge
    charge=0.0
    nc={'ASP':-1,'GLU':-1,'LYS':+1,'ARG':+1,'HIS':+1}
    close.sort()
    print close
    for res in close:
        restype=X.resname(res)
        if nc.has_key(restype):
            charge=charge+nc[restype]
            print res,restype,nc[restype],charge
    print 'Net charge',charge
    return charge

    def setSequencesfromMutationCodes(self, DB=None, callback=None, selected=None):
        """Set the aa sequence using wt ref aa and mutation code
           Assumes mutation code is consistent with ref aa seq"""

        if DB == None:
            return
        proteins = DB.getRecs()
        refprot = DB.meta.refprotein
        refseq = DB[refprot].aaseq
        refaa = self.AAList2String(refseq)
        refpdb = DB[refprot].Structure
        #Create protool oinstance for ref pdb
        import Protool
        Xref = Protool.structureIO()
        Xref.parsepdb(refpdb)        
        for protein in selected:
            rec = DB.get(protein)
            if rec.hasStructure() == 'available':
                continue                
            print 'Protein:', protein            
            #if no sequence try create one from mutation code
            if rec.aaseq == None and rec.Mutations != None:               
                print 'no sequence, using mutation code and ref protein seq'                 
                import PEATSA.Core as Core
                print 'Record has mutation code %s' %rec.Mutations
                mutationSet = Core.Data.MutationSet(rec.Mutations)
                Xref.Remove_All_NonAminoAcids()
                refaa = Core.Data.GetChainSequences(Xref)['A']
                #print refaa
                mutseq = mutationSet.applyToSequence(refaa, id='A', offset=None, pdb=Xref)                 
                rec.aaseq = self.string2AAseq(mutseq)                
        
        return

 def __init__(self):
     import Protool
     self.PI=Protool.structureIO()
     self.aas=self.PI.trueaminoacids.keys()
     self.PI.readpdb('test.pdb')
     #
     import FFF.FFFcontrol as FFFC
     import os, sys
     scriptdir=os.path.split(os.path.abspath(__file__))[0]
     FFFdir=os.path.split(scriptdir)[0]
     Rotamerlib=FFFC.Rotamer_class(os.path.join(FFFdir,'parameters/small_lib'))
     self.FFF=FFFC.FFF()
     self.FFF.read_pdb('test.pdb')
     #self.Model=FFFC.pKa_class(self.FFF,Rotamerlib,os.path.join(FFFdir,'parameters'))
     self.Model=FFFC.model_class(self.FFF,Rotamerlib,os.path.join(FFFdir,'parameters'))
     #
     # Test mutations
     #
     self.mutate_test()
     #self.Model.repair_all()
     #
     # Build all hydrogens - standard protonation state
     #
     #self.Model.build_hydrogens()
     #self.FFF.write_pqr('2lzt.pqr.pdb')
     return

def read_sequences(pdbfiles,newfiles):
    #
    # Load all the PDB files and Make sure that the sequences are ok
    #
    align={}
    allseqs={}
    print 'Reading pdb files and extracting sequences'
    import Protool
    for pdb in pdbfiles:
        pdbfile=newfiles[pdb]
        X=Protool.structureIO_fast()
        print 'Reading: %s' %pdbfile
        X.readpdb(pdbfile)
        X.RemoveALT()
        s_keys=X.residues.keys()
        s_keys.sort()
        #
        # Construct a special list for the sequence
        #
        newlist=[]
        for s in s_keys:
            if X.three_to_one.has_key(X.resname(s)):
                newlist.append([s,X.three_to_one[X.resname(s)]])
        align[pdb]=newlist[:]
        allseqs[pdb]=X.Seq2Pir(None,pdb)
    return align,allseqs

def convert_classic_to_PEAT(operations):
    """Convert a set of classic mutations to a set of PEAT operations
    The classic operations are in the format: A12G+R45V etc."""
    #
    # Do a quick sanity check
    #
    for junk in ['?','unknown','empty']:
        if operations.lower().find(junk)!=-1:
            return False
    #
    # Deal with the operations
    #
    sp=operations.split('+')
    import Protool, string
    P=Protool.structureIO()
    POP=[]
    for op in sp:
        if op=='wt':
            continue
        old=op[0]
        new=op[-1]
        number=int(op[1:-1])
        try:
            POP.append('%s:%s:%s:%s' %('',string.zfill(number,P.length_of_residue_numbers),P.one_to_three[old],P.one_to_three[new]))
        except KeyError:
            return False
    return string.join(POP,'+')

 def checkMutation(self, DB, name, ref=None, X=None):
     """Check mutations based on ref sequence and current mutant
        sequence, should be triggered whenever ref protein is altered so 
        that the mutation codes are updated."""
     prot = DB.get(name)
     if prot.aaseq == None:
         return
     if ref == None:
         ref = self.DB.meta.refprotein            
     refseq = self.AAList2String(DB.get(ref).aaseq)        
     if prot.aaseq == None:
         return
     #get mutations from sequence
     seq = self.AAList2String(prot.aaseq)       
     if seq == refseq:
         return
         
     #get alignment for pdb seq and AA from DNA seq    
     import PEATSA.Core as Core
     if X == None:
         #we need to also provide the ref structure
         import Protool
         X=Protool.structureIO()
         X.parsepdb(DB.get(ref).Structure)
     print X 
     mset = Core.Data.mutationSetFromSequencesAndStructure(refseq, seq, X)
     #prot.Mutations = '+'.join(mset.mutationCodes())
     prot.Mutations = mset.codeString(X)
     return

def Model_Mutations_old(pdbfile,mol2files,mutations,max_overlap=0.5,return_score=False):
    """Model a number of mutations in a pdbfile when one or more ligands are present"""
    #
    # Initialise mutate routines
    #
    MUT=Mutate(max_bump=max_overlap)
    #
    # Read PDB file
    #
    import Protool
    P=Protool.structureIO()
    P.readpdb(pdbfile)
    P.remove_all_hydrogens()
    #
    # Read mol2 file
    #
    L=Protool.ligand(P)
    for mol2file in mol2files:
        print "Added %s with tag 'LIGAND'" %mol2file
        L.readmol2(mol2file,tag='LIGAND')
    #
    # Pass combined pdb file to mutate routines and mutate
    #
    MUT.new_PDB(P)
    import pKa.pKD_tools as pKD_tools
    total_bump=0.0
    #
    # Model
    #
    for mutation in mutations:
        #
        # Get info
        #
        resid=pKD_tools.get_resid_from_mut(mutation)
        newres=pKD_tools.get_newrestyp_from_mut(mutation)
        oldres=pKD_tools.get_oldrestyp_from_mut(mutation)
        bump_score=MUT.mutate(resid,newres,orgtype=oldres)
        if bump_score is None or bump_score is False or bump_score>max_overlap:
            print 'Cannot model this set of mutations - too many bumps'
            return False,20.0
        print 'Bump score for %s: %5.3f' %(mutation,bump_score)
        total_bump=total_bump+bump_score
    print 'Total bump score for all mutations: %5.3f' %(bump_score)
    if return_score:
        return MUT,bump_score
    return MUT

    def remALT(self,pdbfile, environment): 
        import Protool

        x = Protool.structureIO()
        x.readpdb('%s.pdb' % (pdbfile))
        x.RemoveALT()
        x.writepdb('%s.pdb' % (pdbfile), dont_write_HETATMS=1)
        environment.output('Removed alternate residues')

    def remALT(self,pdbfile): 
        import Protool

        x = Protool.structureIO()
        x.readpdb('%s.pdb' % (pdbfile))
        x.RemoveALT()
        x.writepdb('%s.pdb' % (pdbfile), dont_write_HETATMS=1)
        print 'Removed alternate residues'

 def load_groups(self):
     #
     # Get all the titratable groups in the pdb file
     #
     import Protool
     P=Protool.structureIO()
     P.readpdb(self.params['pdb'])
     self.groups=P.get_titratable_groups()
     return

def findSequenceDifferences(child_sequence, parent_sequence,full_parent_sequence,ignoreCterm=False):
    """
    # Find all amino acid differences between child_sequence and parent_sequence
    Child sequence and parent sequence must be aligned and in 1-letter format:
    child_sequence, parent_sequence: AAADEFFG
    full parent sequence is a Protool.sequence object
    """
    #
    # Loop over the sequences - changes are record from parent -> child
    #
    import string
    operations=[]
    import Protool
    PI=Protool.structureIO()
    #
    Cterm_add=0
    insert_num=0
    full_parent_count=0
    #for count in range(len(record_sequence)):
    #    parent_aa=parent_sequence[count]
    #    child_aa=record_sequence[count]
    for parent_aa,child_aa in zip(parent_sequence,child_sequence):
        #
        #print parent_aa,child_aa
        if parent_aa!=child_aa:
            
            # Find the PDB file residue number
            if full_parent_count>=len(full_parent_sequence):
                # If we have an insertion at the Cterm
                aa_identifier=full_parent_sequence[-1][0]
                if ignoreCterm:
                    continue
            else:
                aa_identifier=full_parent_sequence[full_parent_count][0]
            #if aa_identifier[-1]==':':
            #    aa_identifier=aa_identifier[:-1]
            #
            # Convert to 3letter format
            #
            if parent_aa!='-':
                full_parent_count=full_parent_count+1
                parent_aa=PI.one_to_three[parent_aa]
            if child_aa!='-':
                child_aa=PI.one_to_three[child_aa]
            if parent_aa=='-':
                operations.append('insert%d:%s:%s' %(insert_num,aa_identifier,child_aa))
                insert_num=insert_num+1
            elif child_aa=='-':
                insert_num=0
                operations.append('delete:%s:%s' %(aa_identifier,parent_aa))
            else:
                insert_num=0
                operations.append('%s:%s:%s' %(aa_identifier,parent_aa,child_aa))
        else:
            full_parent_count=full_parent_count+1
    return operations

def matrix():
    import os
    dirs=os.listdir('data')
    xs=[]
    ys=[]
    for dir in dirs:
        print 'Processing %s' %dir
        #
        # find the PDB file
        #
        realdir=os.path.join(os.getcwd(),'data',dir)
        files=os.listdir(realdir)
        for file in files:
            realfile=os.path.join(realdir,file)
            if file[-4:]=='.pdb':
                import pKa.pKaTool.pKaIO
                X=pKa.pKaTool.pKaIO.pKaIO(realfile)
                X.assess_status()
                if X.calculation_completed==1:
                    #
                    # Hurra, the calc is complete. Load the matrix
                    #
                    PBEmatrix=X.read_matrix()
                    #
                    # Now calculate the same matrix with Protool
                    #
                    P=Protool.structureIO()
                    P.readpdb(realfile)
                    P.get_titratable_groups()
                    dist_matrix=P.Calculate_matrix(8)
                    #
                    # Plot it
                    #
                    x=[]
                    y=[]
                    for group1 in PBEmatrix.keys():
                        for group2 in PBEmatrix.keys():
                            PBE_ene=PBEmatrix[group1][group2][0]
                            try:
                                new_ene=dist_matrix[group1][group2]
                            except:
                                continue
                            #
                            # Load the values, distances in x, PBE_ene in y
                            #
                            if new_ene and PBE_ene:
                                x.append(abs(new_ene))
                                y.append(abs(PBE_ene))
                    #
                    # Append these result to the big arrays
                    #
                    ys.append(y)
                    xs.append(x)
    plotit(xs,ys,'Matrix',dirs)
    return

def CleanPDB2PQR(inputFile, outputFile, forceField="amber", removeWater=True, removeLigand=True, removeAltLoc=True, addHydrogens=True, correct=True):

	'''Cleans a PDB by using PDB2PQR
	
	See CleanPDB for argument details

	Note: With pdb2pqr you cannot remove-water or ligands.
		
	Errors: 
		Raises an exception if the inputFile is not a valid PDB file.'''

	import Protool

	try:
		command = 'pdb2pqr --chain --ff=%s' % (forceField)

		if removeWater == True:
			print >>sys.stderr, 'Warn: Currently PDB2PQR does not remove waters from pdb files'

		if removeLigand == True:
			print >>sys.stderr, 'Warn: Currently PDB2PQR can not be used to remove heterogens from PDB files'

		if addHydrogens == False:
			print >>sys.stderr, 'Warn: Turning of Hydrogen addition with PDB2PQR automatically turns of rotamer correction'
			command = command + " --clean "
		else:
			if correct == False:
				command = command + " --nodebump "

		#Protool ignores altlocs so we can use it to remove them 
		#Do this first as Protool as when protool reads then writes a pdb2pqr cleaned file 
		#it raises an error on reading it again 
		if removeAltLoc is True:
			pdb = Protool.structureIO()
			pdb.readpdb(inputFile)
			pdb.writepdb(outputFile)

			inputFile = outputFile

		command = command + ' %s %s' % (inputFile, outputFile)
		print 'Using: ', command

		process = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
		stdout, stderr = process.communicate()

		if process.returncode != 0:
			raise ProteinDesignToolException, 'Error using pdb2pqr to clean pdb file %s' % inputFile
		
	except BaseException, data:
		print 'Encountered an exception cleaning pdb %s' % inputFile
		if stdout is not None:
			print 'PDB2PQR output follows:'
			print stdout
		
		raise 			

 def addPDBFile(self, DB=None, name=None, pdbfile=None,
                    pdbdata=None, pdbname=None, gui=True):
     """Add a PDB file to the record given as argument"""
     import os, tkMessageBox
     if pdbdata == None and pdbfile == None:
         savedir=os.getcwd()
         global PDB_code
         pdbfile=tkFileDialog.askopenfilename(defaultextension='.pdb',
                                      initialdir=savedir,
                                      filetypes=[("PDB file","*.pdb"),
                                                 ("PDB file","*.brk"),
                                                 ("All files","*.*")])
         if not pdbfile:
             return
     if pdbfile:    
         pdbname = os.path.basename(pdbfile)
     import Protool
     self.X=Protool.structureIO()
     # Extracting PDB_code from pdbfile
     if pdbdata != None: 
         self.X.readpdb(data=pdbdata)
     elif os.path.isfile(pdbfile):
         PDB_code=pdbfile.split('/').pop().split('.')[0]
         # Try to read it using Protool    
         try:
             self.X.readpdb(filename=pdbfile)
         except:
             tkMessageBox.showwarning('Error reading PDB file',
                                      'I could not read the PDB file. This probably means that the PDB file is corrupt in some way.')
             return
         
     AlignmentMap = None    
     if gui==True:
         if tkMessageBox.askyesno('Reset AA Seq?',
                                  'Do you want to reset the amino acid Sequence?'):      
             AlignmentMap = self.checkPDBSequence(name)               
         
     #store it    
     DB.storePDB(name, self.X, AlignmentMap)
     if hasattr(DB.meta,'refprotein'):                
         ref = DB.meta.refprotein
         #if this is the reference protein remodel mutations and rewrite mut codes   
         if name == ref:
             print name, ref
             print 'rechecking mutation codes, ref prot structure has changed'                
             #get new mutation codes
             import PEATSA.Core as Core
             for p in DB.getRecs():
                 self.checkMutation(DB, p, ref, self.X)
             #self.checkModels(DB)
             
     #add the original pdb name
     DB.data[name]['pdbname'] = pdbname
     return