fixed base lookup bug in dna pairwise alignment and ensured abiguous and non-standard...
[jalview.git] / src / jalview / analysis / AlignSeq.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer
3  * Copyright (C) 2007 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
18  */
19 package jalview.analysis;
20
21 import java.util.*;
22
23 import java.awt.*;
24
25 import jalview.datamodel.*;
26 import jalview.schemes.*;
27 import jalview.util.*;
28
29 /**
30  *
31  *
32  * @author $author$
33  * @version $Revision$
34  */
35 public class AlignSeq
36 {
37   public static final String PEP = "pep";
38   public static final String DNA = "dna";
39   
40   static String[] dna =
41       {
42     "A", "C", "G", "T", "-"};
43     //"C", "T", "A", "G", "-"};
44   static String[] pep =
45       {
46       "A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F",
47       "P", "S", "T", "W", "Y", "V", "B", "Z", "X", "-"
48   };
49   int[][] score;
50   int[][] E;
51   int[][] F;
52   int[][] traceback;
53   int[] seq1;
54   int[] seq2;
55   SequenceI s1;
56   SequenceI s2;
57   public String s1str;
58   public String s2str;
59   int maxi;
60   int maxj;
61   int[] aseq1;
62   int[] aseq2;
63   public String astr1 = "";
64   public String astr2 = "";
65
66   /** DOCUMENT ME!! */
67   public int seq1start;
68
69   /** DOCUMENT ME!! */
70   public int seq1end;
71
72   /** DOCUMENT ME!! */
73   public int seq2start;
74
75   /** DOCUMENT ME!! */
76   public int seq2end;
77   int count;
78
79   /** DOCUMENT ME!! */
80   public int maxscore;
81   float pid;
82   int prev = 0;
83   int gapOpen = 120;
84   int gapExtend = 20;
85   int[][] lookup = ResidueProperties.getBLOSUM62();
86   String[] intToStr = pep;
87   int defInt = 23;
88   StringBuffer output = new StringBuffer();
89   String type;
90
91   /**
92    * Creates a new AlignSeq object.
93    *
94    * @param s1 DOCUMENT ME!
95    * @param s2 DOCUMENT ME!
96    * @param type DOCUMENT ME!
97    */
98   public AlignSeq(SequenceI s1, SequenceI s2, String type)
99   {
100     SeqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(), type);
101   }
102
103   /**
104    * Creates a new AlignSeq object.
105    *
106    * @param s1 DOCUMENT ME!
107    * @param s2 DOCUMENT ME!
108    * @param type DOCUMENT ME!
109    */
110   public AlignSeq(SequenceI s1,
111                   String string1,
112                   SequenceI s2,
113                   String string2,
114                   String type)
115   {
116     SeqInit(s1, string1, s2, string2, type);
117   }
118
119   /**
120    * DOCUMENT ME!
121    *
122    * @return DOCUMENT ME!
123    */
124   public int getMaxScore()
125   {
126     return maxscore;
127   }
128
129   /**
130    * DOCUMENT ME!
131    *
132    * @return DOCUMENT ME!
133    */
134   public int getSeq2Start()
135   {
136     return seq2start;
137   }
138
139   /**
140    * DOCUMENT ME!
141    *
142    * @return DOCUMENT ME!
143    */
144   public int getSeq2End()
145   {
146     return seq2end;
147   }
148
149   /**
150    * DOCUMENT ME!
151    *
152    * @return DOCUMENT ME!
153    */
154   public int getSeq1Start()
155   {
156     return seq1start;
157   }
158
159   /**
160    * DOCUMENT ME!
161    *
162    * @return DOCUMENT ME!
163    */
164   public int getSeq1End()
165   {
166     return seq1end;
167   }
168
169   /**
170    * DOCUMENT ME!
171    *
172    * @return DOCUMENT ME!
173    */
174   public String getOutput()
175   {
176     return output.toString();
177   }
178
179   /**
180    * DOCUMENT ME!
181    *
182    * @return DOCUMENT ME!
183    */
184   public String getAStr1()
185   {
186     return astr1;
187   }
188
189   /**
190    * DOCUMENT ME!
191    *
192    * @return DOCUMENT ME!
193    */
194   public String getAStr2()
195   {
196     return astr2;
197   }
198
199   /**
200    * DOCUMENT ME!
201    *
202    * @return DOCUMENT ME!
203    */
204   public int[] getASeq1()
205   {
206     return aseq1;
207   }
208
209   /**
210    * DOCUMENT ME!
211    *
212    * @return DOCUMENT ME!
213    */
214   public int[] getASeq2()
215   {
216     return aseq2;
217   }
218
219   /**
220    * DOCUMENT ME!
221    *
222    * @return DOCUMENT ME!
223    */
224   public SequenceI getS1()
225   {
226     return s1;
227   }
228
229   /**
230    * DOCUMENT ME!
231    *
232    * @return DOCUMENT ME!
233    */
234   public SequenceI getS2()
235   {
236     return s2;
237   }
238
239   /**
240    * DOCUMENT ME!
241    *
242    * @param s1 DOCUMENT ME!
243    * @param string1 -  string to align for sequence1
244    * @param s2 sequence 2
245    * @param string2 -  string to align for sequence2
246    * @param type DNA or PEPTIDE
247    */
248   public void SeqInit(SequenceI s1,
249                       String string1,
250                       SequenceI s2,
251                       String string2,
252                       String type)
253   {
254     this.s1 = s1;
255     this.s2 = s2;
256     setDefaultParams(type);
257     SeqInit(string1, string2);
258   }
259
260   public void SeqInit(SequenceI s1,
261                       String string1,
262                       SequenceI s2,
263                       String string2,
264                       ScoreMatrix scoreMatrix)
265   {
266     this.s1 = s1;
267     this.s2 = s2;
268     setType(scoreMatrix.isDNA() ? AlignSeq.DNA : AlignSeq.PEP);
269     lookup = scoreMatrix.getMatrix();
270   }
271
272   /**
273    * construct score matrix for string1 and string2 (after removing any existing gaps
274    * @param string1
275    * @param string2
276    */
277   private void SeqInit(String string1, String string2)
278   {
279     s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
280     s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
281
282     if (s1str.length() == 0 || s2str.length() == 0)
283     {
284       output.append("ALL GAPS: " +
285                     (s1str.length() == 0 ? s1.getName() : " ")
286                     + (s2str.length() == 0 ? s2.getName() : ""));
287       return;
288     }
289
290     //System.out.println("lookuip " + rt.freeMemory() + " "+  rt.totalMemory());
291     seq1 = new int[s1str.length()];
292
293     //System.out.println("seq1 " + rt.freeMemory() +" "  + rt.totalMemory());
294     seq2 = new int[s2str.length()];
295
296     //System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());
297     score = new int[s1str.length()][s2str.length()];
298
299     //System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());
300     E = new int[s1str.length()][s2str.length()];
301
302     //System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());
303     F = new int[s1str.length()][s2str.length()];
304     traceback = new int[s1str.length()][s2str.length()];
305
306     //System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());
307     seq1 = stringToInt(s1str, type);
308
309     //System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());
310     seq2 = stringToInt(s2str, type);
311
312     //System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());
313     //   long tstart = System.currentTimeMillis();
314     //    calcScoreMatrix();
315     //long tend = System.currentTimeMillis();
316     //System.out.println("Time take to calculate score matrix = " + (tend-tstart) + " ms");
317     //   printScoreMatrix(score);
318     //System.out.println();
319     //printScoreMatrix(traceback);
320     //System.out.println();
321     //  printScoreMatrix(E);
322     //System.out.println();
323     ///printScoreMatrix(F);
324     //System.out.println();
325     // tstart = System.currentTimeMillis();
326     //traceAlignment();
327     //tend = System.currentTimeMillis();
328     //System.out.println("Time take to traceback alignment = " + (tend-tstart) + " ms");
329   }
330
331   private void setDefaultParams(String type)
332   {
333     setType(type);
334
335     if (type.equals(AlignSeq.PEP))
336     {
337       lookup = ResidueProperties.getDefaultPeptideMatrix();
338     }
339     else if (type.equals(AlignSeq.DNA))
340     {
341       lookup = ResidueProperties.getDefaultDnaMatrix();
342     }
343   }
344
345   private void setType(String type2)
346   {
347     this.type = type2;
348     if (type.equals(AlignSeq.PEP))
349     {
350       intToStr = pep;
351       defInt = 23;
352     }
353     else if (type.equals(AlignSeq.DNA))
354     {
355       intToStr = dna;
356       defInt = 4;
357     }
358     else
359     {
360       output.append("Wrong type = dna or pep only");
361       throw new Error("Unknown Type " + type2 +
362                       " - dna or pep are the only allowed values.");
363     }
364   }
365
366   /**
367    * DOCUMENT ME!
368    */
369   public void traceAlignment()
370   {
371     // Find the maximum score along the rhs or bottom row
372     int max = -9999;
373
374     for (int i = 0; i < seq1.length; i++)
375     {
376       if (score[i][seq2.length - 1] > max)
377       {
378         max = score[i][seq2.length - 1];
379         maxi = i;
380         maxj = seq2.length - 1;
381       }
382     }
383
384     for (int j = 0; j < seq2.length; j++)
385     {
386       if (score[seq1.length - 1][j] > max)
387       {
388         max = score[seq1.length - 1][j];
389         maxi = seq1.length - 1;
390         maxj = j;
391       }
392     }
393
394     //  System.out.println(maxi + " " + maxj + " " + score[maxi][maxj]);
395     int i = maxi;
396     int j = maxj;
397     int trace;
398     maxscore = score[i][j] / 10;
399
400     seq1end = maxi + 1;
401     seq2end = maxj + 1;
402
403     aseq1 = new int[seq1.length + seq2.length];
404     aseq2 = new int[seq1.length + seq2.length];
405
406     count = (seq1.length + seq2.length) - 1;
407
408     while ( (i > 0) && (j > 0))
409     {
410       if ( (aseq1[count] != defInt) && (i >= 0))
411       {
412         aseq1[count] = seq1[i];
413         astr1 = s1str.charAt(i) + astr1;
414       }
415
416       if ( (aseq2[count] != defInt) && (j > 0))
417       {
418         aseq2[count] = seq2[j];
419         astr2 = s2str.charAt(j) + astr2;
420       }
421
422       trace = findTrace(i, j);
423
424       if (trace == 0)
425       {
426         i--;
427         j--;
428       }
429       else if (trace == 1)
430       {
431         j--;
432         aseq1[count] = defInt;
433         astr1 = "-" + astr1.substring(1);
434       }
435       else if (trace == -1)
436       {
437         i--;
438         aseq2[count] = defInt;
439         astr2 = "-" + astr2.substring(1);
440       }
441
442       count--;
443     }
444
445     seq1start = i + 1;
446     seq2start = j + 1;
447
448     if (aseq1[count] != defInt)
449     {
450       aseq1[count] = seq1[i];
451       astr1 = s1str.charAt(i) + astr1;
452     }
453
454     if (aseq2[count] != defInt)
455     {
456       aseq2[count] = seq2[j];
457       astr2 = s2str.charAt(j) + astr2;
458     }
459   }
460
461   /**
462    * DOCUMENT ME!
463    */
464   public void printAlignment(java.io.PrintStream os)
465   {
466     // TODO: Use original sequence characters rather than re-translated characters in output
467     // Find the biggest id length for formatting purposes
468     int maxid = s1.getName().length();
469     if (s2.getName().length() > maxid)
470     {
471       maxid = s2.getName().length();
472     }
473
474     int len = 72 - maxid - 1;
475     int nochunks = ( (aseq1.length - count) / len) + 1;
476     pid = 0;
477
478     output.append("Score = " + score[maxi][maxj] + "\n");
479     output.append("Length of alignment = " + (aseq1.length - count) + "\n");
480     output.append("Sequence ");
481     output.append(new Format("%" + maxid + "s").form(s1.getName()));
482     output.append(" :  " + s1.getStart() + " - " + s1.getEnd() +
483                   " (Sequence length = " +
484                   s1str.length() + ")\n");
485     output.append("Sequence ");
486     output.append(new Format("%" + maxid + "s").form(s2.getName()));
487     output.append(" :  " + s2.getStart() + " - " + s2.getEnd() +
488                   " (Sequence length = " +
489                   s2str.length() + ")\n\n");
490     
491     for (int j = 0; j < nochunks; j++)
492     {
493       // Print the first aligned sequence
494       output.append(new Format("%" + (maxid) + "s").form(s1.getName()) + " ");
495
496       for (int i = 0; i < len; i++)
497       {
498         if ( (i + (j * len)) < astr1.length())
499         {
500           output.append(astr1.charAt(i +
501                                               (j * len)));
502         }
503       }
504
505       output.append("\n");
506       output.append(new Format("%" + (maxid) + "s").form(" ") + " ");
507
508       // Print out the matching chars
509       for (int i = 0; i < len; i++)
510       {
511         if ( (i + (j * len)) < astr1.length())
512         {
513           if (astr1.charAt(i + (j * len))==astr2.charAt(i + (j * len)) &&
514               !jalview.util.Comparison.isGap(astr1.charAt(i + (j * len))))
515           {
516             pid++;
517             output.append("|");
518           }
519           else if (type.equals("pep"))
520           {
521             if (ResidueProperties.getPAM250(
522                     astr1.charAt(i + (j * len)),
523                     astr2.charAt(i + (j * len)))>0)
524             {
525               output.append(".");
526             }
527             else
528             {
529               output.append(" ");
530             }
531           }
532           else
533           {
534             output.append(" ");
535           }
536         }
537       }
538
539       // Now print the second aligned sequence
540       output = output.append("\n");
541       output = output.append(new Format("%" + (maxid) + "s").form(s2.getName()) +
542                              " ");
543
544       for (int i = 0; i < len; i++)
545       {
546         if ( (i + (j * len)) < astr2.length())
547         {
548           output.append(astr2.charAt(i + (j * len)));
549         }
550       }
551
552       output = output.append("\n\n");
553     }
554
555     pid = pid / (float) (aseq1.length - count) * 100;
556     output = output.append(new Format("Percentage ID = %2.2f\n\n").form(pid));
557
558     try
559     {
560       os.print(output.toString());
561     }
562     catch (Exception ex)
563     {}
564   }
565
566   /**
567    * DOCUMENT ME!
568    *
569    * @param mat DOCUMENT ME!
570    */
571   public void printScoreMatrix(int[][] mat)
572   {
573     int n = seq1.length;
574     int m = seq2.length;
575
576     for (int i = 0; i < n; i++)
577     {
578       // Print the top sequence
579       if (i == 0)
580       {
581         Format.print(System.out, "%8s", s2str.substring(0, 1));
582
583         for (int jj = 1; jj < m; jj++)
584         {
585           Format.print(System.out, "%5s", s2str.substring(jj, jj + 1));
586         }
587
588         System.out.println();
589       }
590
591       for (int j = 0; j < m; j++)
592       {
593         if (j == 0)
594         {
595           Format.print(System.out, "%3s", s1str.substring(i, i + 1));
596         }
597
598         Format.print(System.out, "%3d ", mat[i][j] / 10);
599       }
600
601       System.out.println();
602     }
603   }
604
605   /**
606    * DOCUMENT ME!
607    *
608    * @param i DOCUMENT ME!
609    * @param j DOCUMENT ME!
610    *
611    * @return DOCUMENT ME!
612    */
613   public int findTrace(int i, int j)
614   {
615     int t = 0;
616     int max = score[i - 1][j - 1] + (lookup[seq1[i]][seq2[j]] * 10);
617
618     if (F[i][j] > max)
619     {
620       max = F[i][j];
621       t = -1;
622     }
623     else if (F[i][j] == max)
624     {
625       if (prev == -1)
626       {
627         max = F[i][j];
628         t = -1;
629       }
630     }
631
632     if (E[i][j] >= max)
633     {
634       max = E[i][j];
635       t = 1;
636     }
637     else if (E[i][j] == max)
638     {
639       if (prev == 1)
640       {
641         max = E[i][j];
642         t = 1;
643       }
644     }
645
646     prev = t;
647
648     return t;
649   }
650
651   /**
652    * DOCUMENT ME!
653    */
654   public void calcScoreMatrix()
655   {
656     int n = seq1.length;
657     int m = seq2.length;
658
659     // top left hand element
660     score[0][0] = lookup[seq1[0]][seq2[0]] * 10;
661     E[0][0] = -gapExtend;
662     F[0][0] = 0;
663
664     // Calculate the top row first
665     for (int j = 1; j < m; j++)
666     {
667       // What should these values be? 0 maybe
668       E[0][j] = max(score[0][j - 1] - gapOpen, E[0][j - 1] - gapExtend);
669       F[0][j] = -gapExtend;
670
671       score[0][j] = max(lookup[seq1[0]][seq2[j]] * 10, -gapOpen,
672                         -gapExtend);
673
674       traceback[0][j] = 1;
675     }
676
677     // Now do the left hand column
678     for (int i = 1; i < n; i++)
679     {
680       E[i][0] = -gapOpen;
681       F[i][0] = max(score[i - 1][0] - gapOpen, F[i - 1][0] - gapExtend);
682
683       score[i][0] = max(lookup[seq1[i]][seq2[0]] * 10, E[i][0], F[i][0]);
684       traceback[i][0] = -1;
685     }
686
687     // Now do all the other rows
688     for (int i = 1; i < n; i++)
689     {
690       for (int j = 1; j < m; j++)
691       {
692         E[i][j] = max(score[i][j - 1] - gapOpen, E[i][j - 1] -
693                       gapExtend);
694         F[i][j] = max(score[i - 1][j] - gapOpen, F[i - 1][j] -
695                       gapExtend);
696
697         score[i][j] = max(score[i - 1][j - 1] +
698                           (lookup[seq1[i]][seq2[j]] * 10), E[i][j], F[i][j]);
699         traceback[i][j] = findTrace(i, j);
700       }
701     }
702   }
703
704   /**
705    * DOCUMENT ME!
706    *
707    * @param gapChar DOCUMENT ME!
708    * @param seq DOCUMENT ME!
709    *
710    * @return DOCUMENT ME!
711    */
712   public static String extractGaps(String gapChar, String seq)
713   {
714     StringTokenizer str = new StringTokenizer(seq, gapChar);
715     StringBuffer newString = new StringBuffer();
716
717     while (str.hasMoreTokens())
718     {
719       newString.append(str.nextToken());
720     }
721
722     return newString.toString();
723   }
724
725   /**
726    * DOCUMENT ME!
727    *
728    * @param i1 DOCUMENT ME!
729    * @param i2 DOCUMENT ME!
730    * @param i3 DOCUMENT ME!
731    *
732    * @return DOCUMENT ME!
733    */
734   public int max(int i1, int i2, int i3)
735   {
736     int max = i1;
737
738     if (i2 > i1)
739     {
740       max = i2;
741     }
742
743     if (i3 > max)
744     {
745       max = i3;
746     }
747
748     return max;
749   }
750
751   /**
752    * DOCUMENT ME!
753    *
754    * @param i1 DOCUMENT ME!
755    * @param i2 DOCUMENT ME!
756    *
757    * @return DOCUMENT ME!
758    */
759   public int max(int i1, int i2)
760   {
761     int max = i1;
762
763     if (i2 > i1)
764     {
765       max = i2;
766     }
767
768     return max;
769   }
770
771   /**
772    * DOCUMENT ME!
773    *
774    * @param s DOCUMENT ME!
775    * @param type DOCUMENT ME!
776    *
777    * @return DOCUMENT ME!
778    */
779   public int[] stringToInt(String s, String type)
780   {
781     int[] seq1 = new int[s.length()];
782
783     for (int i = 0; i < s.length(); i++)
784     {
785       // String ss = s.substring(i, i + 1).toUpperCase();
786       char c = s.charAt(i);
787       if ('a' <= c && c <= 'z')
788       {
789         // TO UPPERCASE !!!
790         c -= ('a' - 'A');
791       }
792
793       try
794       {
795         if (type.equals("pep"))
796         {
797           seq1[i] = ResidueProperties.aaIndex[c];
798           if (seq1[i] > 23)
799           {
800             seq1[i] = 23;
801           }
802         }
803         else if (type.equals("dna"))
804         {
805           seq1[i] = ResidueProperties.nucleotideIndex[c];
806           if (seq1[i] > 4)
807           {
808             seq1[i] = 4;
809           }
810         }
811
812       }
813       catch (Exception e)
814       {
815         if (type.equals("dna"))
816         {
817           seq1[i] = 4;
818         }
819         else
820         {
821           seq1[i] = 23;
822         }
823       }
824     }
825
826     return seq1;
827   }
828
829   /**
830    * DOCUMENT ME!
831    *
832    * @param g DOCUMENT ME!
833    * @param mat DOCUMENT ME!
834    * @param n DOCUMENT ME!
835    * @param m DOCUMENT ME!
836    * @param psize DOCUMENT ME!
837    */
838   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
839                                    int psize)
840   {
841     int max = -1000;
842     int min = 1000;
843
844     for (int i = 0; i < n; i++)
845     {
846       for (int j = 0; j < m; j++)
847       {
848         if (mat[i][j] >= max)
849         {
850           max = mat[i][j];
851         }
852
853         if (mat[i][j] <= min)
854         {
855           min = mat[i][j];
856         }
857       }
858     }
859
860     System.out.println(max + " " + min);
861
862     for (int i = 0; i < n; i++)
863     {
864       for (int j = 0; j < m; j++)
865       {
866         int x = psize * i;
867         int y = psize * j;
868
869         //      System.out.println(mat[i][j]);
870         float score = (float) (mat[i][j] - min) / (float) (max - min);
871         g.setColor(new Color(score, 0, 0));
872         g.fillRect(x, y, psize, psize);
873
874         //      System.out.println(x + " " + y + " " + score);
875       }
876     }
877   }
878 }