extract sequence from DNA bases like 'DG' (JAL-866)
[jalview.git] / src / MCview / PDBChain.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.6)
3  * Copyright (C) 2010 J Procter, AM Waterhouse, G Barton, M Clamp, S Searle
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
10  * 
11  * Jalview is distributed in the hope that it will be useful, but 
12  * WITHOUT ANY WARRANTY; without even the implied warranty 
13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
14  * PURPOSE.  See the GNU General Public License for more details.
15  * 
16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
17  */
18 package MCview;
19
20 import java.util.*;
21
22 import java.awt.*;
23
24 import jalview.analysis.*;
25 import jalview.datamodel.*;
26 import jalview.schemes.*;
27 import jalview.structure.StructureMapping;
28
29 public class PDBChain
30 {
31   /**
32    * SequenceFeature group for PDB File features added to sequences
33    */
34   private static final String PDBFILEFEATURE = "PDBFile";
35
36   private static final String IEASTATUS = "IEA:jalview";
37
38   public String id;
39
40   public Vector bonds = new Vector();
41
42   public Vector atoms = new Vector();
43
44   public Vector residues = new Vector();
45
46   public int offset;
47
48   public Sequence sequence;
49
50   public boolean isNa = false;
51
52   public boolean isVisible = true;
53
54   public int pdbstart = 0;
55
56   public int pdbend = 0;
57
58   public int seqstart = 0;
59
60   public int seqend = 0;
61
62   public String pdbid = "";
63
64   public PDBChain(String pdbid, String id)
65   {
66     this.pdbid = pdbid.toLowerCase();
67     this.id = id;
68   }
69
70   /**
71    * character used to write newlines
72    */
73   protected String newline = System.getProperty("line.separator");
74
75   public void setNewlineString(String nl)
76   {
77     newline = nl;
78   }
79
80   public String getNewlineString()
81   {
82     return newline;
83   }
84
85   public String print()
86   {
87     String tmp = "";
88
89     for (int i = 0; i < bonds.size(); i++)
90     {
91       tmp = tmp + ((Bond) bonds.elementAt(i)).at1.resName + " "
92               + ((Bond) bonds.elementAt(i)).at1.resNumber + " " + offset
93               + newline;
94     }
95
96     return tmp;
97   }
98
99   /**
100    * Annotate the residues with their corresponding positions in s1 using the
101    * alignment in as
102    * 
103    * @param as
104    * @param s1
105    */
106   public void makeExactMapping(AlignSeq as, SequenceI s1)
107   {
108     int pdbpos = as.getSeq2Start() - 2;
109     int alignpos = s1.getStart() + as.getSeq1Start() - 3;
110
111     for (int i = 0; i < as.astr1.length(); i++)
112     {
113       if (as.astr1.charAt(i) != '-')
114       {
115         alignpos++;
116       }
117
118       if (as.astr2.charAt(i) != '-')
119       {
120         pdbpos++;
121       }
122
123       if (as.astr1.charAt(i) == as.astr2.charAt(i))
124       {
125         Residue res = (Residue) residues.elementAt(pdbpos);
126         Enumeration en = res.atoms.elements();
127         while (en.hasMoreElements())
128         {
129           Atom atom = (Atom) en.nextElement();
130           atom.alignmentMapping = alignpos;
131         }
132       }
133     }
134   }
135
136   /**
137    * copy over the RESNUM seqfeatures from the internal chain sequence to the
138    * mapped sequence
139    * 
140    * @param seq
141    * @param status
142    *          The Status of the transferred annotation
143    * @return the features added to sq (or its dataset)
144    */
145   public SequenceFeature[] transferRESNUMFeatures(SequenceI seq,
146           String status)
147   {
148     SequenceI sq = seq;
149     while (sq != null && sq.getDatasetSequence() != null)
150     {
151       sq = sq.getDatasetSequence();
152       if (sq == sequence)
153       {
154         return null;
155       }
156     }
157     /**
158      * Remove any existing features for this chain if they exist ?
159      * SequenceFeature[] seqsfeatures=seq.getSequenceFeatures(); int
160      * totfeat=seqsfeatures.length; // Remove any features for this exact chain
161      * ? for (int i=0; i<seqsfeatures.length; i++) { }
162      */
163     if (status == null)
164     {
165       status = PDBChain.IEASTATUS;
166     }
167     SequenceFeature[] features = sequence.getSequenceFeatures();
168     for (int i = 0; i < features.length; i++)
169     {
170       if (features[i].getFeatureGroup().equals(pdbid))
171       {
172         SequenceFeature tx = new SequenceFeature(features[i]);
173         tx.setBegin(1 + ((Atom) ((Residue) residues.elementAt(tx.getBegin()
174                 - offset)).atoms.elementAt(0)).alignmentMapping);
175         tx.setEnd(1 + ((Atom) ((Residue) residues.elementAt(tx.getEnd()
176                 - offset)).atoms.elementAt(0)).alignmentMapping);
177         tx.setStatus(status
178                 + ((tx.getStatus() == null || tx.getStatus().length() == 0) ? ""
179                         : ":" + tx.getStatus()));
180         if (tx.begin != 0 && tx.end != 0)
181           sq.addSequenceFeature(tx);
182       }
183     }
184     return features;
185   }
186
187   public void makeCaBondList()
188   {
189     boolean na = false;
190     int numNa = 0;
191     for (int i = 0; i < (residues.size() - 1); i++)
192     {
193       Residue tmpres = (Residue) residues.elementAt(i);
194       Residue tmpres2 = (Residue) residues.elementAt(i + 1);
195       Atom at1 = tmpres.findAtom("CA");
196       Atom at2 = tmpres2.findAtom("CA");
197       na = false;
198       if ((at1 == null) && (at2 == null))
199       {
200         na = true;
201         at1 = tmpres.findAtom("P");
202         at2 = tmpres2.findAtom("P");
203       }
204       if ((at1 != null) && (at2 != null))
205       {
206         if (at1.chain.equals(at2.chain))
207         {
208           if (na)
209           {
210             numNa++;
211           }
212           makeBond(at1, at2);
213         }
214       }
215       else
216       {
217         System.out.println("not found " + i);
218       }
219     }
220     if (numNa > 0 && ((numNa / residues.size()) > 0.99))
221     {
222       isNa = true;
223     }
224   }
225
226   public void makeBond(Atom at1, Atom at2)
227   {
228     float[] start = new float[3];
229     float[] end = new float[3];
230
231     start[0] = at1.x;
232     start[1] = at1.y;
233     start[2] = at1.z;
234
235     end[0] = at2.x;
236     end[1] = at2.y;
237     end[2] = at2.z;
238
239     bonds.addElement(new Bond(start, end, at1, at2));
240   }
241
242   public void makeResidueList()
243   {
244     int count = 0;
245     Object symbol;
246     boolean deoxyn=false;
247     boolean nucleotide = false;
248     StringBuffer seq = new StringBuffer();
249     Vector resFeatures = new Vector();
250     Vector resAnnotation = new Vector();
251     int i, iSize = atoms.size() - 1;
252     int resNumber = -1;
253     for (i = 0; i <= iSize; i++)
254     {
255       Atom tmp = (Atom) atoms.elementAt(i);
256       resNumber = tmp.resNumber;
257       int res = resNumber;
258
259       if (i == 0)
260       {
261         offset = resNumber;
262       }
263
264       Vector resAtoms = new Vector();
265       // Add atoms to a vector while the residue number
266       // remains the same as the first atom's resNumber (res)
267       while ((resNumber == res) && (i < atoms.size()))
268       {
269         resAtoms.addElement((Atom) atoms.elementAt(i));
270         i++;
271
272         if (i < atoms.size())
273         {
274           resNumber = ((Atom) atoms.elementAt(i)).resNumber;
275         }
276         else
277         {
278           resNumber++;
279         }
280       }
281
282       // We need this to keep in step with the outer for i = loop
283       i--;
284
285       // Make a new Residue object with the new atoms vector
286       residues.addElement(new Residue(resAtoms, resNumber - 1, count));
287
288       Residue tmpres = (Residue) residues.lastElement();
289       Atom tmpat = (Atom) tmpres.atoms.elementAt(0);
290       // Make A new SequenceFeature for the current residue numbering
291       SequenceFeature sf = new SequenceFeature("RESNUM", tmpat.resName
292               + ":" + tmpat.resNumIns + " " + pdbid + id, "", offset
293               + count, offset + count, pdbid);
294       // MCview.PDBChain.PDBFILEFEATURE);
295       resFeatures.addElement(sf);
296       resAnnotation.addElement(new Annotation(tmpat.tfactor));
297       // Keep totting up the sequence
298       if ((symbol = ResidueProperties.getAA3Hash().get(tmpat.resName)) == null)
299       {
300         String nucname = tmpat.resName.trim();
301         // use the aaIndex rather than call 'toLower' - which would take a bit more time.
302         deoxyn=nucname.length()==2 && ResidueProperties.aaIndex[nucname.charAt(0)]==ResidueProperties.aaIndex['D'];
303         if (tmpat.name.equalsIgnoreCase("CA")
304                 || ResidueProperties.nucleotideIndex[nucname.charAt((deoxyn ? 1 : 0))] == -1)
305         {
306           seq.append("X");
307           // System.err.println("PDBReader:Null aa3Hash for " +
308           // tmpat.resName);
309         }
310         else
311         {
312           // nucleotide flag
313           nucleotide = true;
314           seq.append(nucname.charAt((deoxyn ? 1 : 0)));
315         }
316       }
317       else
318       {
319         if (nucleotide)
320         {
321           System.err
322                   .println("Warning: mixed nucleotide and amino acid chain.. its gonna do bad things to you!");
323         }
324         seq.append(ResidueProperties.aa[((Integer) symbol).intValue()]);
325       }
326       count++;
327     }
328
329     if (id.length() < 1)
330     {
331       id = " ";
332     }
333     isNa = nucleotide;
334     sequence = new Sequence(id, seq.toString(), offset, resNumber - 1); // Note:
335     // resNumber-offset
336     // ~=
337     // seq.size()
338     // Add normalised feature scores to RESNUM indicating start/end of sequence
339     // sf.setScore(offset+count);
340
341     // System.out.println("PDB Sequence is :\nSequence = " + seq);
342     // System.out.println("No of residues = " + residues.size());
343     for (i = 0, iSize = resFeatures.size(); i < iSize; i++)
344     {
345       sequence.addSequenceFeature((SequenceFeature) resFeatures
346               .elementAt(i));
347       resFeatures.setElementAt(null, i);
348     }
349     Annotation[] annots = new Annotation[resAnnotation.size()];
350     float max = 0;
351     for (i = 0, iSize = annots.length; i < iSize; i++)
352     {
353       annots[i] = (Annotation) resAnnotation.elementAt(i);
354       if (annots[i].value > max)
355         max = annots[i].value;
356       resAnnotation.setElementAt(null, i);
357     }
358     AlignmentAnnotation tfactorann = new AlignmentAnnotation(
359             "PDB.TempFactor", "Temperature Factor for "
360                     + sequence.getName(), annots, 0, max,
361             AlignmentAnnotation.LINE_GRAPH);
362     tfactorann.setSequenceRef(sequence);
363     sequence.addAlignmentAnnotation(tfactorann);
364   }
365
366   public void setChargeColours()
367   {
368     for (int i = 0; i < bonds.size(); i++)
369     {
370       try
371       {
372         Bond b = (Bond) bonds.elementAt(i);
373
374         if (b.at1.resName.equalsIgnoreCase("ASP")
375                 || b.at1.resName.equalsIgnoreCase("GLU"))
376         {
377           b.startCol = Color.red;
378         }
379         else if (b.at1.resName.equalsIgnoreCase("LYS")
380                 || b.at1.resName.equalsIgnoreCase("ARG"))
381         {
382           b.startCol = Color.blue;
383         }
384         else if (b.at1.resName.equalsIgnoreCase("CYS"))
385         {
386           b.startCol = Color.yellow;
387         }
388         else
389         {
390           b.startCol = Color.lightGray;
391         }
392
393         if (b.at2.resName.equalsIgnoreCase("ASP")
394                 || b.at2.resName.equalsIgnoreCase("GLU"))
395         {
396           b.endCol = Color.red;
397         }
398         else if (b.at2.resName.equalsIgnoreCase("LYS")
399                 || b.at2.resName.equalsIgnoreCase("ARG"))
400         {
401           b.endCol = Color.blue;
402         }
403         else if (b.at2.resName.equalsIgnoreCase("CYS"))
404         {
405           b.endCol = Color.yellow;
406         }
407         else
408         {
409           b.endCol = Color.lightGray;
410         }
411       } catch (Exception e)
412       {
413         Bond b = (Bond) bonds.elementAt(i);
414         b.startCol = Color.gray;
415         b.endCol = Color.gray;
416       }
417     }
418   }
419
420   public void setChainColours(jalview.schemes.ColourSchemeI cs)
421   {
422     Bond b;
423     int index;
424     for (int i = 0; i < bonds.size(); i++)
425     {
426       try
427       {
428         b = (Bond) bonds.elementAt(i);
429
430         index = ((Integer) ResidueProperties.aa3Hash.get(b.at1.resName))
431                 .intValue();
432         b.startCol = cs.findColour(ResidueProperties.aa[index].charAt(0));
433
434         index = ((Integer) ResidueProperties.aa3Hash.get(b.at2.resName))
435                 .intValue();
436         b.endCol = cs.findColour(ResidueProperties.aa[index].charAt(0));
437
438       } catch (Exception e)
439       {
440         b = (Bond) bonds.elementAt(i);
441         b.startCol = Color.gray;
442         b.endCol = Color.gray;
443       }
444     }
445   }
446
447   public void setChainColours(Color col)
448   {
449     for (int i = 0; i < bonds.size(); i++)
450     {
451       Bond tmp = (Bond) bonds.elementAt(i);
452       tmp.startCol = col;
453       tmp.endCol = col;
454     }
455   }
456
457   public AlignmentAnnotation[] transferResidueAnnotation(SequenceI seq,
458           String status)
459   {
460     AlignmentAnnotation[] transferred = null;
461
462     return transferred;
463
464   }
465
466   /**
467    * copy any sequence annotation onto the sequence mapped using the provided
468    * StructureMapping
469    * 
470    * @param mapping
471    */
472   public void transferResidueAnnotation(StructureMapping mapping)
473   {
474     SequenceI sq = mapping.getSequence();
475     if (sq != null)
476     {
477       if (sequence != null && sequence.getAnnotation() != null)
478       {
479
480       }
481       float min = -1, max = 0;
482       Annotation[] an = new Annotation[sq.getEnd() - sq.getStart() + 1];
483       for (int i = sq.getStart(), j = sq.getEnd(), k = 0; i <= j; i++, k++)
484       {
485         int prn = mapping.getPDBResNum(k + 1);
486
487         an[k] = new Annotation((float) prn);
488         if (min == -1)
489         {
490           min = k;
491           max = k;
492         }
493         else
494         {
495           if (min > k)
496           {
497             min = k;
498           }
499           else if (max < k)
500           {
501             max = k;
502           }
503         }
504       }
505       sq.addAlignmentAnnotation(new AlignmentAnnotation("PDB.RESNUM",
506               "PDB Residue Numbering for " + this.pdbid + ":" + this.id,
507               an, (float) min, (float) max, AlignmentAnnotation.LINE_GRAPH));
508     }
509   }
510 }