Package Bio :: Package PDB :: Module DSSP'
[hide private]
[frames] | no frames]

Source Code for Module Bio.PDB.DSSP'

  1  # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk) 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  5   
  6  from __future__ import generators 
  7   
  8  import os 
  9  import tempfile 
 10  from Bio.PDB import * 
 11  from PDBExceptions import PDBException 
 12  from AbstractPropertyMap import AbstractResiduePropertyMap 
 13  import re 
 14   
 15   
 16  __doc__=""" 
 17  Use the DSSP program to calculate secondary structure and accessibility. 
 18  You need to have a working version of DSSP (and a license, free for  
 19  academic use) in order to use this. For DSSP, see U{http://www.cmbi.kun.nl/gv/dssp/}. 
 20   
 21  The DSSP codes for secondary structure used here are: 
 22   
 23      - H        Alpha helix (4-12) 
 24      - B        Isolated beta-bridge residue 
 25      - E        Strand 
 26      - G        3-10 helix 
 27      - I        pi helix 
 28      - T        Turn 
 29      - S        Bend 
 30      - -        None 
 31  """ 
 32   
 33  # Match C in DSSP 
 34  _dssp_cys=re.compile('[a-z]') 
 35   
 36  # Maximal ASA of amino acids 
 37  # Values from Sander & Rost, (1994), Proteins, 20:216-226 
 38  # Used for relative accessibility 
 39  MAX_ACC={} 
 40  MAX_ACC["ALA"]=106.0 
 41  MAX_ACC["CYS"]=135.0 
 42  MAX_ACC["ASP"]=163.0 
 43  MAX_ACC["GLU"]=194.0 
 44  MAX_ACC["PHE"]=197.0 
 45  MAX_ACC["GLY"]=84.0 
 46  MAX_ACC["HIS"]=184.0 
 47  MAX_ACC["ILE"]=169.0 
 48  MAX_ACC["LYS"]=205.0 
 49  MAX_ACC["LEU"]=164.0 
 50  MAX_ACC["MET"]=188.0 
 51  MAX_ACC["ASN"]=157.0 
 52  MAX_ACC["PRO"]=136.0 
 53  MAX_ACC["GLN"]=198.0 
 54  MAX_ACC["ARG"]=248.0 
 55  MAX_ACC["SER"]=130.0 
 56  MAX_ACC["THR"]=142.0 
 57  MAX_ACC["VAL"]=142.0 
 58  MAX_ACC["TRP"]=227.0 
 59  MAX_ACC["TYR"]=222.0 
 60   
 61   
62 -def ss_to_index(ss):
63 """ 64 Secondary structure symbol to index. 65 H=0 66 E=1 67 C=2 68 """ 69 if ss=='H': 70 return 0 71 if ss=='E': 72 return 1 73 if ss=='C': 74 return 2 75 assert(0)
76
77 -def dssp_dict_from_pdb_file(in_file, DSSP="dssp"):
78 """ 79 Create a DSSP dictionary from a PDB file. 80 81 Example: 82 >>> dssp_dict=dssp_dict_from_pdb_file("1fat.pdb") 83 >>> aa, ss, acc=dssp_dict[('A', 1)] 84 85 @param in_file: pdb file 86 @type in_file: string 87 88 @param DSSP: DSSP executable (argument to os.system) 89 @type DSSP: string 90 91 @return: a dictionary that maps (chainid, resid) to 92 amino acid type, secondary structure code and 93 accessibility. 94 @rtype: {} 95 """ 96 out_file=tempfile.mktemp() 97 os.system(DSSP+" %s > %s" % (in_file, out_file)) 98 dict, keys=make_dssp_dict(out_file) 99 # This can be dangerous... 100 #os.system("rm "+out_file) 101 return dict, keys
102
103 -def make_dssp_dict(filename):
104 """ 105 Return a DSSP dictionary that maps (chainid, resid) to 106 aa, ss and accessibility, from a DSSP file. 107 108 @param filename: the DSSP output file 109 @type filename: string 110 """ 111 dssp={} 112 fp=open(filename, "r") 113 start=0 114 keys=[] 115 for l in fp.readlines(): 116 sl=l.split() 117 if sl[1]=="RESIDUE": 118 # start 119 start=1 120 continue 121 if not start: 122 continue 123 if l[9]==" ": 124 # skip -- missing residue 125 continue 126 resseq=int(l[5:10]) 127 icode=l[10] 128 chainid=l[11] 129 aa=l[13] 130 ss=l[16] 131 if ss==" ": 132 ss="-" 133 acc=int(l[34:38]) 134 res_id=(" ", resseq, icode) 135 dssp[(chainid, res_id)]=(aa, ss, acc) 136 keys.append((chainid, res_id)) 137 fp.close() 138 return dssp, keys
139 140
141 -class DSSP(AbstractResiduePropertyMap):
142 """ 143 Run DSSP on a pdb file, and provide a handle to the 144 DSSP secondary structure and accessibility. 145 146 Note that DSSP can only handle one model. 147 148 Example: 149 >>> p=PDBParser() 150 >>> structure=parser.get_structure("1fat.pdb") 151 >>> model=structure[0] 152 >>> dssp=DSSP(model, "1fat.pdb") 153 >>> # print dssp data for a residue 154 >>> secondary_structure, accessibility=dssp[(chain_id, res_id)] 155 """
156 - def __init__(self, model, pdb_file, dssp="dssp"):
157 """ 158 @param model: the first model of the structure 159 @type model: L{Model} 160 161 @param pdb_file: a PDB file 162 @type pdb_file: string 163 164 @param dssp: the dssp executable (ie. the argument to os.system) 165 @type dssp: string 166 """ 167 # create DSSP dictionary 168 dssp_dict, dssp_keys=dssp_dict_from_pdb_file(pdb_file, dssp) 169 dssp_map={} 170 dssp_list=[] 171 # Now create a dictionary that maps Residue objects to 172 # secondary structure and accessibility, and a list of 173 # (residue, (secondary structure, accessibility)) tuples 174 for key in dssp_keys: 175 chain_id, res_id=key 176 chain=model[chain_id] 177 res=chain[res_id] 178 aa, ss, acc=dssp_dict[key] 179 res.xtra["SS_DSSP"]=ss 180 res.xtra["EXP_DSSP_ASA"]=acc 181 # relative accessibility 182 resname=res.get_resname() 183 rel_acc=acc/MAX_ACC[resname] 184 if rel_acc>1.0: 185 rel_acc=1.0 186 res.xtra["EXP_DSSP_RASA"]=rel_acc 187 # Verify if AA in DSSP == AA in Structure 188 # Something went wrong if this is not true! 189 resname=to_one_letter_code[resname] 190 if resname=="C": 191 # DSSP renames C in C-bridges to a,b,c,d,... 192 # - we rename it back to 'C' 193 if _dssp_cys.match(aa): 194 aa='C' 195 if not (resname==aa): 196 raise PDBException, "Structure/DSSP mismatch at "+str(res) 197 dssp_map[key]=((res, ss, acc, rel_acc)) 198 dssp_list.append((res, ss, acc, rel_acc)) 199 AbstractResiduePropertyMap.__init__(self, dssp_map, dssp_keys, dssp_list)
200 201 202 if __name__=="__main__": 203 204 import sys 205 206 p=PDBParser() 207 s=p.get_structure('X', sys.argv[1]) 208 209 model=s[0] 210 211 d=DSSP(model, sys.argv[1]) 212 213 for r in d: 214 print r 215 216 print d.keys() 217 218 print len(d) 219 220 print d.has_key(('A', 1)) 221 222 print d[('A', 1)] 223 224 print s[0]['A'][1].xtra 225