ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/owl/trunk/proteinstructure/PdbfilePdb.java
Revision: 237
Committed: Thu Jul 19 14:23:35 2007 UTC (17 years, 2 months ago) by duarte
File size: 10485 byte(s)
Log Message:
Now rmsd calculation works for 2 Pdb objects where there are missing residues (unobserved) or missing atoms
New Exception ConformationNotSameSizeError
Still not printing out individual distances between 2 conformations (altghough it is possible)
Can't specify yet a restricted range of residues to calculate rmsd only with them
Line User Rev File contents
1 duarte 207 package proteinstructure;
2    
3     import java.io.BufferedReader;
4     import java.io.File;
5     import java.io.FileNotFoundException;
6     import java.io.FileReader;
7     import java.io.IOException;
8     import java.util.ArrayList;
9     import java.util.Collections;
10     import java.util.HashMap;
11 duarte 222 import java.util.TreeMap;
12 duarte 207 import java.util.regex.Matcher;
13     import java.util.regex.Pattern;
14    
15 duarte 226 import javax.vecmath.Point3d;
16    
17 duarte 207 public class PdbfilePdb extends Pdb {
18    
19 duarte 208 private static final String UNKNOWN_STRING ="Unknown";
20     private static final String NULL_chainCode = "A";
21    
22 duarte 207 private String pdbfile;
23    
24     /**
25     * Constructs Pdb object reading from pdb file given pdb chain code. Model will be DEFAULT_MODEL
26     * @param pdbfile
27     * @param pdbChainCode
28     * @throws FileNotFoundException
29     * @throws IOException
30 duarte 208 * @throws PdbfileFormatError
31 stehr 215 * @throws PdbChainCodeNotFoundError
32 duarte 207 */
33 stehr 215 public PdbfilePdb (String pdbfile, String pdbChainCode) throws FileNotFoundException, IOException, PdbfileFormatError, PdbChainCodeNotFoundError {
34 duarte 207 this(pdbfile,pdbChainCode,DEFAULT_MODEL);
35     }
36    
37     /**
38     * Constructs Pdb object reading from pdb file given pdb chain code and model serial
39     * @param pdbfile
40     * @param pdbChainCode
41     * @param model_serial
42     * @throws FileNotFoundException
43     * @throws IOException
44 duarte 208 * @throws PdbfileFormatError
45 stehr 215 * @throws PdbChainCodeNotFoundError
46 duarte 207 */
47 stehr 215 public PdbfilePdb (String pdbfile, String pdbChainCode, int model_serial) throws FileNotFoundException, IOException, PdbfileFormatError, PdbChainCodeNotFoundError {
48 duarte 207 this.pdbfile = pdbfile;
49     this.model=model_serial;
50 duarte 208 this.pdbCode=UNKNOWN_STRING; // we initialise to unknown in case we don't find it in pdb file
51 stehr 217 this.pdbChainCode=pdbChainCode.toUpperCase(); // our convention: chain codes are upper case
52 duarte 208 // we set chainCode to pdbChainCode except for case NULL where we use "A"
53     this.chainCode=pdbChainCode;
54     if (pdbChainCode.equals("NULL")) this.chainCode=NULL_chainCode;
55    
56     this.sequence=""; // we initialise it to empty string, then is set inread_pdb_data_from_file
57    
58 duarte 222 // we initialise the resser2secstruct and secstruct2resinterval Maps to empty, if no sec structure info is found then it remains empty
59 duarte 219 this.resser2secstruct = new HashMap<Integer, String>();
60 duarte 222 this.secstruct2resinterval = new TreeMap<String, Interval>();
61 duarte 219
62 duarte 207 read_pdb_data_from_file();
63 duarte 208
64     // when reading from pdb file we have no information of residue numbers or author's (original) pdb residue number, so we fill the mapping with the residue numbers we know
65     //TODO eventually we could assign our own internal residue numbers when reading from pdb and thus this map would be used
66     this.resser2pdbresser = new HashMap<Integer, String>();
67     this.pdbresser2resser = new HashMap<String, Integer>();
68     for (int resser:resser2restype.keySet()){
69     resser2pdbresser.put(resser, String.valueOf(resser));
70     pdbresser2resser.put(String.valueOf(resser), resser);
71     }
72 duarte 219
73 duarte 237 // initialising atomser2atom from resser_atom2atomserial
74     atomser2atom = new HashMap<Integer, String>();
75     for (String resser_atom:resser_atom2atomserial.keySet()){
76     int atomserial = resser_atom2atomserial.get(resser_atom);
77     String atom = resser_atom.split("_")[1];
78     atomser2atom.put(atomserial,atom);
79     }
80 duarte 207 }
81    
82    
83    
84     /**
85 duarte 219 * To read the pdb data (atom coordinates, residue serials, atom serials) from file.
86     * chainCode gets set to same as pdbChainCode, except if input chain code NULL then chainCode will be 'A'
87 duarte 207 * pdbCode gets set to the one parsed in HEADER or to 'Unknown' if not found
88     * sequence gets set to the sequence read from ATOM lines (i.e. observed resdiues only)
89 duarte 219 * No insertion codes are parsed or taken into account at the moment. Thus files with
90     * insertion codes will be incorrectly read
91 duarte 207 * @param pdbfile
92     * @throws FileNotFoundException
93     * @throws IOException
94 duarte 208 * @throws PdbfileFormatError
95 stehr 215 * @throws PdbChainCodeNotFoundError
96 duarte 207 */
97 stehr 215 private void read_pdb_data_from_file() throws FileNotFoundException, IOException, PdbfileFormatError, PdbChainCodeNotFoundError {
98 duarte 207 resser_atom2atomserial = new HashMap<String,Integer>();
99     resser2restype = new HashMap<Integer,String>();
100 duarte 226 atomser2coord = new HashMap<Integer,Point3d>();
101 duarte 207 atomser2resser = new HashMap<Integer,Integer>();
102     boolean empty = true; // controls whether we don't find any atom line for given pdbChainCode and model
103 duarte 208 // we set chainCodeStr (for regex) to pdbChainCode except for case NULL where we use " " (NULL is a blank chain code in pdb files)
104     String chainCodeStr=pdbChainCode;
105     if (pdbChainCode.equals("NULL")) chainCodeStr=" ";
106    
107 duarte 207 int thismodel=DEFAULT_MODEL; // we initialise to DEFAULT_MODEL, in case file doesn't have MODEL lines
108     BufferedReader fpdb = new BufferedReader(new FileReader(new File(pdbfile)));
109 duarte 208 int linecount=0;
110 duarte 207 String line;
111     while ((line = fpdb.readLine() ) != null ) {
112 duarte 208 linecount++;
113 duarte 219 // HEADER
114 duarte 207 Pattern p = Pattern.compile("^HEADER");
115     Matcher m = p.matcher(line);
116     if (m.find()){
117     Pattern ph = Pattern.compile("^HEADER.{56}(\\d\\w{3})");
118     Matcher mh = ph.matcher(line);
119     if (mh.find()) {
120     pdbCode=mh.group(1).toLowerCase();
121     }
122 duarte 208 } else if (linecount==1) { // header not found and we are in line 1
123     // a HEADER is the minimum we ask at the moment for a pdb file to have, if we don't find it in line 1 we throw an exception
124     throw new PdbfileFormatError("The pdb file "+pdbfile+" doesn't seem to have the right format");
125 duarte 207 }
126 duarte 219 // SECONDARY STRUCTURE
127     // helix
128     //HELIX 1 1 LYS A 17 LEU A 26 1
129     // helix ser beg res ser end res ser
130     p = Pattern.compile("^HELIX..(...).{9}"+chainCodeStr+".(....).{6}"+chainCodeStr+".(....)");
131     m = p.matcher(line);
132     if (m.find()){
133     int serial = Integer.valueOf(m.group(1).trim());
134     int beg = Integer.valueOf(m.group(2).trim());
135     int end = Integer.valueOf(m.group(3).trim());
136     String ssId = "H"+serial;
137 duarte 222 secstruct2resinterval.put(ssId, new Interval(beg,end));
138 duarte 219 for (int i=beg;i<=end;i++){
139     if (resser2secstruct.containsKey(i)){// if already assigned we print a warning and then assign it
140     System.err.println("Inconsistency in secondary structure assignment. " +
141     "Residue "+i+" is getting reassigned from "+resser2secstruct.get(i)+" to "+ssId);
142     }
143     resser2secstruct.put(i,ssId);
144     }
145     }
146     // sheet
147     //SHEET 2 A 5 ILE A 96 THR A 99 -1 N LYS A 98 O THR A 107
148     // strand ser sheet id beg res ser end res ser
149     p = Pattern.compile("^SHEET..(...).(...).{7}"+chainCodeStr+"(....).{6}"+chainCodeStr+"(....)");
150     m = p.matcher(line);
151     if (m.find()){
152     int strandSerial = Integer.valueOf(m.group(1).trim());
153     String sheetId = m.group(2).trim();
154     int beg = Integer.valueOf(m.group(3).trim());
155     int end = Integer.valueOf(m.group(4).trim());
156     String ssId = "S"+sheetId+strandSerial;
157 duarte 222 secstruct2resinterval.put(ssId, new Interval(beg,end));
158 duarte 219 for (int i=beg;i<=end;i++){
159     if (resser2secstruct.containsKey(i)){// if already assigned we print a warning and then assign it
160     System.err.println("Inconsistency in secondary structure assignment. " +
161     "Residue "+i+" is getting reassigned from "+resser2secstruct.get(i)+" to "+ssId);
162     }
163     resser2secstruct.put(i,ssId);
164     }
165     }
166     // we've stored the sec structure info in the strands2begEnd and sheets2strands maps.
167     // the assignment to resser2secstruct is done when we reach the ATOM lines, see below
168     // turn
169     //TURN 1 S1A GLY A 16 GLN A 18 SURFACE
170     // turn ser beg res ser end res ser
171     p = Pattern.compile("^TURN...(...).{9}"+chainCodeStr+"(....).{6}"+chainCodeStr+"(....)");
172     m = p.matcher(line);
173     if (m.find()){
174     int serial = Integer.valueOf(m.group(1).trim());
175     int beg = Integer.valueOf(m.group(2).trim());
176     int end = Integer.valueOf(m.group(3).trim());
177     String ssId = "T"+serial;
178 duarte 222 secstruct2resinterval.put(ssId, new Interval(beg,end));
179 duarte 219 for (int i=beg;i<=end;i++){
180     if (resser2secstruct.containsKey(i)){// if already assigned we print a warning and then assign it
181     System.err.println("Inconsistency in secondary structure assignment. " +
182     "Residue "+i+" is getting reassigned from "+resser2secstruct.get(i)+" to "+ssId);
183     }
184     resser2secstruct.put(i,ssId);
185     }
186     }
187     // MODEL
188 duarte 207 p = Pattern.compile("^MODEL\\s+(\\d+)");
189     m = p.matcher(line);
190     if (m.find()){
191     thismodel=Integer.parseInt(m.group(1));
192     }
193     if (thismodel!=model) continue; // we skip reading of atom lines if we are not in the desired model
194 duarte 219 // ATOM
195 duarte 207 p = Pattern.compile("^ATOM");
196     m = p.matcher(line);
197     if (m.find()){
198 duarte 208 // serial atom res_type chain res_ser x y z
199     Pattern pl = Pattern.compile(".{6}(.....).{2}(...).{1}(...).{1}"+chainCodeStr+"(.{4}).{4}(.{8})(.{8})(.{8})",Pattern.CASE_INSENSITIVE);
200 duarte 207 Matcher ml = pl.matcher(line);
201     if (ml.find()) {
202     empty=false;
203     int atomserial=Integer.parseInt(ml.group(1).trim());
204     String atom = ml.group(2).trim();
205     String res_type = ml.group(3).trim();
206     int res_serial = Integer.parseInt(ml.group(4).trim());
207     double x = Double.parseDouble(ml.group(5).trim());
208     double y = Double.parseDouble(ml.group(6).trim());
209     double z = Double.parseDouble(ml.group(7).trim());
210 duarte 226 Point3d coords = new Point3d(x,y,z);
211 duarte 207 ArrayList<String> aalist=AA.aas();
212     if (aalist.contains(res_type)) {
213     atomser2coord.put(atomserial, coords);
214     atomser2resser.put(atomserial, res_serial);
215     resser2restype.put(res_serial, res_type);
216     ArrayList<String> atomlist = aas2atoms.get(res_type);
217     if (atomlist.contains(atom)){
218     resser_atom2atomserial.put(res_serial+"_"+atom, atomserial);
219     }
220     }
221     }
222     }
223     }
224     fpdb.close();
225 stehr 215 if (empty) {
226     //System.err.println("Couldn't find any atom line for given pdbChainCode: "+pdbChainCode+", model: "+model);
227     throw new PdbChainCodeNotFoundError("Couldn't find any ATOM line for given pdbChainCode: "+pdbChainCode+", model: "+model);
228     }
229 duarte 207 // now we read the sequence from the resser2restype HashMap
230     // NOTE: we must make sure elsewhere that there are no unobserved residues, we can't check that here!
231     ArrayList<Integer> ressers = new ArrayList<Integer>();
232     for (int resser:resser2restype.keySet()) {
233     ressers.add(resser);
234     }
235     Collections.sort(ressers);
236     for (int resser:ressers){
237     String oneletter = AA.threeletter2oneletter(resser2restype.get(resser));
238     sequence += oneletter;
239     }
240     }
241 duarte 219
242 duarte 207 }