general pairwise substitution score models for pairwise alignments and trees.
[jalview.git] / src / jalview / analysis / AlignSeq.java
1 /*
2 * Jalview - A Sequence Alignment Editor and Viewer
3 * Copyright (C) 2006 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 jalview.datamodel.SequenceI;
22
23 import jalview.schemes.*;
24
25 import jalview.util.*;
26
27 import java.awt.*;
28
29 import java.util.*;
30
31
32 /**
33  *
34  *
35  * @author $author$
36  * @version $Revision$
37  */
38 public class AlignSeq
39 {
40     public static final String PEP = "pep";
41     public static final String DNA = "dna";
42     /** DOCUMENT ME!! */
43     public static java.util.Hashtable dnaHash = new java.util.Hashtable();
44
45     static
46     {
47         dnaHash.put("C", new Integer(0));
48         dnaHash.put("T", new Integer(1));
49         dnaHash.put("A", new Integer(2));
50         dnaHash.put("G", new Integer(3));
51         dnaHash.put("-", new Integer(4));
52     }
53
54     static String[] dna = { "C", "T", "A", "G", "-" };
55     static String[] pep =
56     {
57         "A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F",
58         "P", "S", "T", "W", "Y", "V", "B", "Z", "X", "-"
59     };
60     int[][] score;
61     int[][] E;
62     int[][] F;
63     int[][] traceback;
64     int[] seq1;
65     int[] seq2;
66     SequenceI s1;
67     SequenceI s2;
68     public String s1str;
69     public String s2str;
70     int maxi;
71     int maxj;
72     int[] aseq1;
73     int[] aseq2;
74     public String astr1="";
75     public String astr2="";
76
77     /** DOCUMENT ME!! */
78     public int seq1start;
79
80     /** DOCUMENT ME!! */
81     public int seq1end;
82
83     /** DOCUMENT ME!! */
84     public int seq2start;
85
86     /** DOCUMENT ME!! */
87     public int seq2end;
88     int count;
89
90     /** DOCUMENT ME!! */
91     public int maxscore;
92     float pid;
93     int prev = 0;
94     int gapOpen = 120;
95     int gapExtend = 20;
96     int[][] lookup = ResidueProperties.getBLOSUM62();
97     String[] intToStr = pep;
98     int defInt = 23;
99     StringBuffer output = new StringBuffer();
100     String type;
101
102     /**
103      * Creates a new AlignSeq object.
104      *
105      * @param s1 DOCUMENT ME!
106      * @param s2 DOCUMENT ME!
107      * @param type DOCUMENT ME!
108      */
109     public AlignSeq(SequenceI s1, SequenceI s2, String type)
110     {
111         SeqInit(s1, s1.getSequenceAsString(), s2,  s2.getSequenceAsString(), type);
112     }
113
114     /**
115      * Creates a new AlignSeq object.
116      *
117      * @param s1 DOCUMENT ME!
118      * @param s2 DOCUMENT ME!
119      * @param type DOCUMENT ME!
120      */
121     public AlignSeq(SequenceI s1,
122                     String string1,
123                     SequenceI s2,
124                     String string2,
125                     String type)
126     {
127         SeqInit(s1, string1, s2,  string2,  type);
128     }
129
130     /**
131      * DOCUMENT ME!
132      *
133      * @return DOCUMENT ME!
134      */
135     public int getMaxScore()
136     {
137         return maxscore;
138     }
139
140     /**
141      * DOCUMENT ME!
142      *
143      * @return DOCUMENT ME!
144      */
145     public int getSeq2Start()
146     {
147         return seq2start;
148     }
149
150     /**
151      * DOCUMENT ME!
152      *
153      * @return DOCUMENT ME!
154      */
155     public int getSeq2End()
156     {
157         return seq2end;
158     }
159
160     /**
161      * DOCUMENT ME!
162      *
163      * @return DOCUMENT ME!
164      */
165     public int getSeq1Start()
166     {
167         return seq1start;
168     }
169
170     /**
171      * DOCUMENT ME!
172      *
173      * @return DOCUMENT ME!
174      */
175     public int getSeq1End()
176     {
177         return seq1end;
178     }
179
180     /**
181      * DOCUMENT ME!
182      *
183      * @return DOCUMENT ME!
184      */
185     public String getOutput()
186     {
187         return output.toString();
188     }
189
190     /**
191      * DOCUMENT ME!
192      *
193      * @return DOCUMENT ME!
194      */
195     public String getAStr1()
196     {
197         return astr1;
198     }
199
200     /**
201      * DOCUMENT ME!
202      *
203      * @return DOCUMENT ME!
204      */
205     public String getAStr2()
206     {
207         return astr2;
208     }
209
210     /**
211      * DOCUMENT ME!
212      *
213      * @return DOCUMENT ME!
214      */
215     public int[] getASeq1()
216     {
217         return aseq1;
218     }
219
220     /**
221      * DOCUMENT ME!
222      *
223      * @return DOCUMENT ME!
224      */
225     public int[] getASeq2()
226     {
227         return aseq2;
228     }
229
230     /**
231      * DOCUMENT ME!
232      *
233      * @return DOCUMENT ME!
234      */
235     public SequenceI getS1()
236     {
237         return s1;
238     }
239
240     /**
241      * DOCUMENT ME!
242      *
243      * @return DOCUMENT ME!
244      */
245     public SequenceI getS2()
246     {
247         return s2;
248     }
249
250     /**
251      * DOCUMENT ME!
252      *
253      * @param s1 DOCUMENT ME!
254      * @param string1 -  string to align for sequence1
255      * @param s2 sequence 2
256      * @param string2 -  string to align for sequence2
257      * @param type DNA or PEPTIDE
258      */
259     public void SeqInit(SequenceI s1,
260                         String string1,
261                         SequenceI s2,
262                         String string2,
263                         String type)
264     {
265       this.s1 = s1;
266       this.s2 = s2;
267       setDefaultParams(type);
268       SeqInit(string1,string2);
269     }
270     public void SeqInit(SequenceI s1,
271         String string1,
272         SequenceI s2,
273         String string2,
274         ScoreMatrix scoreMatrix) { 
275       this.s1=s1;
276       this.s2=s2;
277       setType(scoreMatrix.isDNA() ? AlignSeq.DNA : AlignSeq.PEP);
278       lookup=scoreMatrix.getMatrix();
279     }
280     /**
281      * construct score matrix for string1 and string2 (after removing any existing gaps
282      * @param string1
283      * @param string2
284      */
285     private void SeqInit(String string1, String string2) {
286       s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
287       s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
288
289       if(s1str.length()==0 || s2str.length()==0)
290       {
291         output.append("ALL GAPS: " +
292                            (s1str.length()==0?s1.getName():" ")
293                            +(s2str.length()==0?s2.getName():""));
294         return;
295       }
296
297       //System.out.println("lookuip " + rt.freeMemory() + " "+  rt.totalMemory());
298       seq1 = new int[s1str.length()];
299
300       //System.out.println("seq1 " + rt.freeMemory() +" "  + rt.totalMemory());
301       seq2 = new int[s2str.length()];
302
303       //System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());
304       score = new int[s1str.length()][s2str.length()];
305
306       //System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());
307       E = new int[s1str.length()][s2str.length()];
308
309       //System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());
310       F = new int[s1str.length()][s2str.length()];
311       traceback = new int[s1str.length()][s2str.length()];
312
313       //System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());
314       seq1 = stringToInt(s1str, type);
315
316       //System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());
317       seq2 = stringToInt(s2str, type);
318
319       //System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());
320       //   long tstart = System.currentTimeMillis();
321       //    calcScoreMatrix();
322       //long tend = System.currentTimeMillis();
323       //System.out.println("Time take to calculate score matrix = " + (tend-tstart) + " ms");
324       //   printScoreMatrix(score);
325       //System.out.println();
326       //printScoreMatrix(traceback);
327       //System.out.println();
328       //  printScoreMatrix(E);
329       //System.out.println();
330       ///printScoreMatrix(F);
331       //System.out.println();
332       // tstart = System.currentTimeMillis();
333       //traceAlignment();
334       //tend = System.currentTimeMillis();
335       //System.out.println("Time take to traceback alignment = " + (tend-tstart) + " ms");
336     }
337
338     private void setDefaultParams(String type) {
339         setType(type);
340         
341         if (type.equals(AlignSeq.PEP))
342         {
343             lookup = ResidueProperties.getDefaultPeptideMatrix();
344         }
345         else if (type.equals(AlignSeq.DNA))
346         {
347             lookup = ResidueProperties.getDefaultDnaMatrix();
348         }
349     }
350
351     private void setType(String type2) {
352       this.type = type2;
353       if (type.equals(AlignSeq.PEP))
354       {
355           intToStr = pep;
356           defInt = 23;
357       }
358       else if (type.equals(AlignSeq.DNA))
359       {
360           intToStr = dna;
361           defInt = 4;
362       }
363       else
364       {
365           output.append("Wrong type = dna or pep only");
366           throw new Error("Unknown Type "+type2+" - dna or pep are the only allowed values.");
367       }
368     }
369
370     /**
371      * DOCUMENT ME!
372      */
373     public void traceAlignment()
374     {
375         // Find the maximum score along the rhs or bottom row
376         int max = -9999;
377
378         for (int i = 0; i < seq1.length; i++)
379         {
380             if (score[i][seq2.length - 1] > max)
381             {
382                 max = score[i][seq2.length - 1];
383                 maxi = i;
384                 maxj = seq2.length - 1;
385             }
386         }
387
388         for (int j = 0; j < seq2.length; j++)
389         {
390             if (score[seq1.length - 1][j] > max)
391             {
392                 max = score[seq1.length - 1][j];
393                 maxi = seq1.length - 1;
394                 maxj = j;
395             }
396         }
397
398         //  System.out.println(maxi + " " + maxj + " " + score[maxi][maxj]);
399         int i = maxi;
400         int j = maxj;
401         int trace;
402         maxscore = score[i][j] / 10;
403
404         seq1end = maxi + 1;
405         seq2end = maxj + 1;
406
407         aseq1 = new int[seq1.length + seq2.length];
408         aseq2 = new int[seq1.length + seq2.length];
409
410         count = (seq1.length + seq2.length) - 1;
411
412         while ((i > 0) && (j > 0))
413         {
414             if ((aseq1[count] != defInt) && (i >= 0))
415             {
416                 aseq1[count] = seq1[i];
417                 astr1 = intToStr[seq1[i]] + astr1;
418             }
419
420             if ((aseq2[count] != defInt) && (j > 0))
421             {
422                 aseq2[count] = seq2[j];
423                 astr2 = intToStr[seq2[j]] + astr2;
424             }
425
426             trace = findTrace(i, j);
427
428             if (trace == 0)
429             {
430                 i--;
431                 j--;
432             }
433             else if (trace == 1)
434             {
435                 j--;
436                 aseq1[count] = defInt;
437                 astr1 = "-" + astr1.substring(1);
438             }
439             else if (trace == -1)
440             {
441                 i--;
442                 aseq2[count] = defInt;
443                 astr2 = "-" + astr2.substring(1);
444             }
445
446             count--;
447         }
448
449         seq1start = i + 1;
450         seq2start = j + 1;
451
452         if (aseq1[count] != defInt)
453         {
454             aseq1[count] = seq1[i];
455             astr1 = intToStr[seq1[i]] + astr1;
456         }
457
458         if (aseq2[count] != defInt)
459         {
460             aseq2[count] = seq2[j];
461             astr2 = intToStr[seq2[j]] + astr2;
462         }
463     }
464
465     /**
466      * DOCUMENT ME!
467      */
468     public void printAlignment(java.io.PrintStream os)
469     {
470         // Find the biggest id length for formatting purposes
471         int maxid = s1.getName().length();
472
473         if (s2.getName().length() > maxid)
474         {
475             maxid = s2.getName().length();
476         }
477
478         int len = 72 - maxid - 1;
479         int nochunks = ((aseq1.length - count) / len) + 1;
480         pid = 0;
481
482         output.append("Score = " + score[maxi][maxj] + "\n");
483         output.append("Length of alignment = " + (aseq1.length - count) + "\n");
484         output.append("Sequence ");
485         output.append(new Format("%" + maxid + "s").form(s1.getName()));
486         output.append(" :  " + s1.getStart() + " - " + s1.getEnd() + " (Sequence length = " +
487             s1str.length() + ")\n");
488         output .append("Sequence ");
489         output.append(new Format("%" + maxid + "s").form(s2.getName()));
490         output.append(" :  " + s2.getStart() + " - " + s2.getEnd() + " (Sequence length = " +
491             s2str.length() + ")\n\n");
492
493         for (int j = 0; j < nochunks; j++)
494         {
495             // Print the first aligned sequence
496             output.append(new Format("%" + (maxid) + "s").form(s1.getName()) + " ");
497
498             for (int i = 0; i < len; i++)
499             {
500                 if ((count + i + (j * len)) < aseq1.length)
501                 {
502                     output.append(new Format("%s").form(intToStr[aseq1[count + i +
503                             (j * len)]]));
504                 }
505             }
506
507             output.append("\n");
508             output.append(new Format("%" + (maxid) + "s").form(" ") + " ");
509
510             // Print out the matching chars
511             for (int i = 0; i < len; i++)
512             {
513                 if ((count + i + (j * len)) < aseq1.length)
514                 {
515                     if (intToStr[aseq1[count + i + (j * len)]].equals(
516                                 intToStr[aseq2[count + i + (j * len)]]) &&
517                             !intToStr[aseq1[count + i + (j * len)]].equals("-"))
518                     {
519                         pid++;
520                         output.append("|");
521                     }
522                     else if (type.equals("pep"))
523                     {
524                         if (ResidueProperties.getPAM250(
525                                     intToStr[aseq1[count + i + (j * len)]],
526                                     intToStr[aseq2[count + i + (j * len)]]) > 0)
527                         {
528                             output.append(".");
529                         }
530                         else
531                         {
532                             output.append(" ");
533                         }
534                     }
535                     else
536                     {
537                         output.append(" ");
538                     }
539                 }
540             }
541
542             // Now print the second aligned sequence
543             output = output.append("\n");
544             output = output.append(new Format("%" + (maxid) + "s").form(s2.getName()) + " ");
545
546             for (int i = 0; i < len; i++)
547             {
548                 if ((count + i + (j * len)) < aseq1.length)
549                 {
550                     output .append(new Format("%s").form(intToStr[aseq2[count + i +
551                             (j * len)]]));
552                 }
553             }
554
555             output = output .append("\n\n");
556         }
557
558         pid = pid / (float) (aseq1.length - count) * 100;
559         output = output.append(new Format("Percentage ID = %2.2f\n\n").form(pid));
560
561         try{
562           os.print(output.toString());
563         }catch(Exception ex){}
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
706     /**
707      * DOCUMENT ME!
708      *
709      * @param gapChar DOCUMENT ME!
710      * @param seq DOCUMENT ME!
711      *
712      * @return DOCUMENT ME!
713      */
714     public static String extractGaps(String gapChar, String seq)
715     {
716         StringTokenizer str = new StringTokenizer(seq, gapChar);
717         StringBuffer newString = new StringBuffer();
718
719         while (str.hasMoreTokens())
720         {
721             newString.append( str.nextToken() );
722         }
723
724         return newString.toString();
725     }
726
727     /**
728      * DOCUMENT ME!
729      *
730      * @param i1 DOCUMENT ME!
731      * @param i2 DOCUMENT ME!
732      * @param i3 DOCUMENT ME!
733      *
734      * @return DOCUMENT ME!
735      */
736     public int max(int i1, int i2, int i3)
737     {
738         int max = i1;
739
740         if (i2 > i1)
741         {
742             max = i2;
743         }
744
745         if (i3 > max)
746         {
747             max = i3;
748         }
749
750         return max;
751     }
752
753     /**
754      * DOCUMENT ME!
755      *
756      * @param i1 DOCUMENT ME!
757      * @param i2 DOCUMENT ME!
758      *
759      * @return DOCUMENT ME!
760      */
761     public int max(int i1, int i2)
762     {
763         int max = i1;
764
765         if (i2 > i1)
766         {
767             max = i2;
768         }
769
770         return max;
771     }
772
773     /**
774      * DOCUMENT ME!
775      *
776      * @param s DOCUMENT ME!
777      * @param type DOCUMENT ME!
778      *
779      * @return DOCUMENT ME!
780      */
781     public int[] stringToInt(String s, String type)
782     {
783         int[] seq1 = new int[s.length()];
784
785         for (int i = 0; i < s.length(); i++)
786         {
787            // String ss = s.substring(i, i + 1).toUpperCase();
788             char c = s.charAt(i);
789             if ('a' <= c && c <= 'z')
790             {
791               // TO UPPERCASE !!!
792               c -= ('a' - 'A');
793             }
794
795
796             try
797             {
798                 if (type.equals("pep"))
799                 {
800                     seq1[i] = ResidueProperties.aaIndex[c];
801                 }
802                 else if (type.equals("dna"))
803                 {
804                     seq1[i] = ResidueProperties.nucleotideIndex[c];
805                 }
806
807                 if (seq1[i] > 23)
808                 {
809                     seq1[i] = 23;
810                 }
811             }
812             catch (Exception e)
813             {
814                 if (type.equals("dna"))
815                 {
816                     seq1[i] = 4;
817                 }
818                 else
819                 {
820                     seq1[i] = 23;
821                 }
822             }
823         }
824
825         return seq1;
826     }
827
828     /**
829      * DOCUMENT ME!
830      *
831      * @param g DOCUMENT ME!
832      * @param mat DOCUMENT ME!
833      * @param n DOCUMENT ME!
834      * @param m DOCUMENT ME!
835      * @param psize DOCUMENT ME!
836      */
837     public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
838         int psize)
839     {
840         int max = -1000;
841         int min = 1000;
842
843         for (int i = 0; i < n; i++)
844         {
845             for (int j = 0; j < m; j++)
846             {
847                 if (mat[i][j] >= max)
848                 {
849                     max = mat[i][j];
850                 }
851
852                 if (mat[i][j] <= min)
853                 {
854                     min = mat[i][j];
855                 }
856             }
857         }
858
859         System.out.println(max + " " + min);
860
861         for (int i = 0; i < n; i++)
862         {
863             for (int j = 0; j < m; j++)
864             {
865                 int x = psize * i;
866                 int y = psize * j;
867
868                 //      System.out.println(mat[i][j]);
869                 float score = (float) (mat[i][j] - min) / (float) (max - min);
870                 g.setColor(new Color(score, 0, 0));
871                 g.fillRect(x, y, psize, psize);
872
873                 //      System.out.println(x + " " + y + " " + score);
874             }
875         }
876     }
877 }