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