Back to index

python-biopython  1.60
NACCESS.py
Go to the documentation of this file.
00001 # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
00002 # This code is part of the Biopython distribution and governed by its
00003 # license.  Please see the LICENSE file that should have been included
00004 # as part of this package.
00005 
00006 # NACCESS interface adapted from Bio/PDB/DSSP.py
00007 
00008 import os, sys, tempfile
00009 from Bio.PDB.PDBIO import PDBIO
00010 from Bio.PDB.AbstractPropertyMap import AbstractResiduePropertyMap, AbstractAtomPropertyMap
00011 
00012 """Interface for the program NACCESS.
00013 
00014 See: http://wolf.bms.umist.ac.uk/naccess/
00015 
00016 errors likely to occur with the binary:
00017 default values are often due to low default settings in accall.pars
00018 - e.g. max cubes error: change in accall.pars and recompile binary
00019 
00020 use naccess -y, naccess -h or naccess -w to include HETATM records
00021 """
00022 
00023 def run_naccess(model, pdb_file, probe_size = None, z_slice = None, \
00024                 naccess = 'naccess', temp_path = '/tmp/'):
00025     
00026     # make temp directory; chdir to temp directory, 
00027     # as NACCESS writes to current working directory
00028     tmp_path = tempfile.mktemp(dir = temp_path)
00029     os.mkdir(tmp_path)
00030     old_dir = os.getcwd()
00031     os.chdir(tmp_path)
00032     
00033     # file name must end with '.pdb' to work with NACCESS
00034     # -> create temp file of existing pdb
00035     #    or write model to temp file
00036     tmp_pdb_file = tempfile.mktemp('.pdb', dir = tmp_path)
00037     if pdb_file:
00038         os.system('cp %s %s' % (pdb_file, tmp_pdb_file))
00039     else:
00040         writer = PDBIO()
00041         writer.set_structure(model.get_parent())
00042         writer.save(tmp_pdb_file)
00043 
00044     # create the command line and run
00045     # catch standard out & err
00046     command = '%s %s ' % (naccess, tmp_pdb_file)
00047     if probe_size:
00048         command += '-p %s ' % probe_size
00049     if z_slice:
00050         command += '-z %s ' % z_slice
00051     in_, out, err = os.popen3(command)
00052     in_.close()
00053     stdout = out.readlines()
00054     out.close()
00055     stderr = err.readlines()
00056     err.close()
00057 
00058     # get the output, then delete the temp directory
00059     rsa_file = tmp_pdb_file[:-4] + '.rsa'
00060     rf = open(rsa_file)
00061     rsa_data = rf.readlines()
00062     rf.close()
00063     asa_file = tmp_pdb_file[:-4] + '.asa'
00064     af = open(asa_file)
00065     asa_data = af.readlines()
00066     af.close()
00067     os.chdir(old_dir)
00068     os.system('rm -rf %s >& /dev/null' % tmp_path)
00069     return rsa_data, asa_data
00070 
00071 def process_rsa_data(rsa_data):
00072     # process the .rsa output file: residue level SASA data
00073     naccess_rel_dict = {}
00074     for line in rsa_data:
00075         if line.startswith('RES'):
00076             res_name = line[4:7]
00077             chain_id = line[8]
00078             resseq = int(line[9:13])
00079             icode = line[13]
00080             res_id = (' ', resseq, icode)
00081             naccess_rel_dict[(chain_id, res_id)] = { \
00082                 'res_name': res_name,
00083                 'all_atoms_abs': float(line[16:22]),
00084                 'all_atoms_rel': float(line[23:28]),
00085                 'side_chain_abs': float(line[29:35]),
00086                 'side_chain_rel': float(line[36:41]),
00087                 'main_chain_abs': float(line[42:48]),
00088                 'main_chain_rel': float(line[49:54]),
00089                 'non_polar_abs': float(line[55:61]),
00090                 'non_polar_rel': float(line[62:67]),
00091                 'all_polar_abs': float(line[68:74]),
00092                 'all_polar_rel': float(line[75:80]) } 
00093     return naccess_rel_dict
00094 
00095 def process_asa_data(rsa_data):
00096     # process the .asa output file: atomic level SASA data
00097     naccess_atom_dict = {}
00098     for line in rsa_data:
00099         atom_serial = line[6:11]
00100         full_atom_id = line[12:16]
00101         atom_id = full_atom_id.strip()
00102         altloc = line[16]
00103         resname = line[17:20]
00104         chainid = line[21]
00105         resseq = int(line[22:26])
00106         icode = line[26]
00107         res_id = (' ', resseq, icode)
00108         id = (chainid, res_id, atom_id)
00109         asa = line[54:62]               # solvent accessibility in Angstrom^2
00110         vdw = line[62:68]               # van der waal radius
00111         naccess_atom_dict[id] = asa
00112     return naccess_atom_dict
00113 
00114 
00115 class NACCESS(AbstractResiduePropertyMap):
00116     
00117     def __init__(self, model, pdb_file = None,
00118                  naccess_binary = 'naccess', tmp_directory = '/tmp'):
00119         res_data, atm_data = run_naccess(model, pdb_file, naccess = naccess_binary,
00120                                          temp_path = tmp_directory)
00121         naccess_dict = process_rsa_data(res_data)
00122         res_list = []
00123         property_dict={}
00124         property_keys=[]
00125         property_list=[]
00126         # Now create a dictionary that maps Residue objects to accessibility
00127         for chain in model:
00128             chain_id=chain.get_id()
00129             for res in chain:
00130                 res_id=res.get_id()
00131                 if (chain_id, res_id) in naccess_dict:
00132                     item = naccess_dict[(chain_id, res_id)]
00133                     res_name = item['res_name']
00134                     assert (res_name == res.get_resname())
00135                     property_dict[(chain_id, res_id)] = item
00136                     property_keys.append((chain_id, res_id))
00137                     property_list.append((res, item))
00138                     res.xtra["EXP_NACCESS"]=item
00139                 else:
00140                     pass
00141         AbstractResiduePropertyMap.__init__(self, property_dict, property_keys, 
00142                 property_list)
00143 
00144 class NACCESS_atomic(AbstractAtomPropertyMap):
00145 
00146     def __init__(self, model, pdb_file = None,
00147                  naccess_binary = 'naccess', tmp_directory = '/tmp'):
00148         res_data, atm_data = run_naccess(model, pdb_file, naccess = naccess_binary,
00149                                          temp_path = tmp_directory)
00150         self.naccess_atom_dict = process_asa_data(atm_data)
00151         atom_list = []
00152         property_dict={}
00153         property_keys=[]
00154         property_list=[]
00155         # Now create a dictionary that maps Atom objects to accessibility
00156         for chain in model:
00157             chain_id = chain.get_id()
00158             for residue in chain:
00159                 res_id = residue.get_id()
00160                 for atom in residue:
00161                     atom_id = atom.get_id()
00162                     full_id=(chain_id, res_id, atom_id)
00163                     if full_id in self.naccess_atom_dict:
00164                         asa = self.naccess_atom_dict[full_id]
00165                         property_dict[full_id]=asa
00166                         property_keys.append((full_id))
00167                         property_list.append((atom, asa))
00168                         atom.xtra['EXP_NACCESS']=asa
00169         AbstractAtomPropertyMap.__init__(self, property_dict, property_keys, 
00170                 property_list)
00171 
00172 
00173 if __name__=="__main__":
00174     
00175     import sys
00176     from Bio.PDB import PDBParser
00177     
00178     p=PDBParser()
00179     s=p.get_structure('X', sys.argv[1])
00180     model=s[0]
00181 
00182     n = NACCESS(model, sys.argv[1])
00183     for e in n.get_iterator():
00184         print e
00185