various loadTree bits for loading and saving trees.
[jalview.git] / src / jalview / analysis / AlignSeq.java
1 /* Jalview - a java multiple alignment editor\r
2  * Copyright (C) 1998  Michele Clamp\r
3  *\r
4  * This program is free software; you can redistribute it and/or\r
5  * modify it under the terms of the GNU General Public License\r
6  * as published by the Free Software Foundation; either version 2\r
7  * of the License, or (at your option) any later version.\r
8  *\r
9  * This program is distributed in the hope that it will be useful,\r
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
12  * GNU General Public License for more details.\r
13  *\r
14  * You should have received a copy of the GNU General Public License\r
15  * along with this program; if not, write to the Free Software\r
16  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.\r
17  */\r
18 package jalview.analysis;\r
19 \r
20 import jalview.schemes.*;\r
21 import jalview.datamodel.SequenceI;\r
22 import jalview.util.*;\r
23 import jalview.io.*;\r
24 \r
25 import java.util.*;\r
26 import java.io.*;\r
27 import java.awt.*;\r
28 \r
29 public class AlignSeq {\r
30 \r
31   int[][] score;\r
32   int[][] E;\r
33   int[][] F;\r
34   int[][] traceback;\r
35 \r
36   int[] seq1;\r
37   int[] seq2;\r
38 \r
39   SequenceI s1;\r
40   SequenceI s2;\r
41 \r
42   String s1str;\r
43   String s2str;\r
44 \r
45   int maxi;\r
46   int maxj;\r
47 \r
48   int[] aseq1;\r
49   int[] aseq2;\r
50 \r
51   String astr1 = "";\r
52   String astr2 = "";\r
53 \r
54   public int seq1start;\r
55   public int seq1end;\r
56   public int seq2start;\r
57   public int seq2end;\r
58 \r
59   int count;\r
60 \r
61   public int maxscore;\r
62   float pid;\r
63   int prev = 0;\r
64 \r
65   public static java.util.Hashtable dnaHash = new java.util.Hashtable();\r
66 \r
67   static  {\r
68     dnaHash.put("C", new Integer(0));\r
69     dnaHash.put("T", new Integer(1));\r
70     dnaHash.put("A", new Integer(2));\r
71     dnaHash.put("G", new Integer(3));\r
72     dnaHash.put("-", new Integer(4));\r
73   }\r
74 \r
75   static String dna[] = {"C","T","A","G","-"};\r
76   static String pep[] = {"A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V","B","Z","X","-"};\r
77 \r
78   int gapOpen = 120;\r
79   int gapExtend = 20;\r
80 \r
81   int lookup[][] = ResidueProperties.getBLOSUM62();\r
82   String intToStr[] = pep;\r
83   int defInt = 23;\r
84 \r
85   String output = "";\r
86 \r
87   String type;\r
88   Runtime rt;\r
89   public AlignSeq() {}\r
90 \r
91   public AlignSeq(SequenceI s1, SequenceI s2,String type) {\r
92     rt  = Runtime.getRuntime();\r
93     SeqInit(s1,s2,type);\r
94   }\r
95 \r
96   public int getMaxScore() {\r
97     return maxscore;\r
98   }\r
99 \r
100   public int getSeq2Start() {\r
101     return seq2start;\r
102   }\r
103 \r
104   public int getSeq2End() {\r
105     return seq2end;\r
106   }\r
107 \r
108   public int getSeq1Start() {\r
109     return seq1start;\r
110   }\r
111 \r
112   public int getSeq1End() {\r
113     return seq1end;\r
114   }\r
115 \r
116   public String getOutput() {\r
117     return output;\r
118   }\r
119 \r
120   public String getAStr1() {\r
121     return astr1;\r
122   }\r
123   public String getAStr2() {\r
124     return astr2;\r
125   }\r
126   public int [] getASeq1() {\r
127     return aseq1;\r
128   }\r
129   public int [] getASeq2() {\r
130     return aseq2;\r
131   }\r
132   public SequenceI getS1() {\r
133     return s1;\r
134   }\r
135   public SequenceI getS2() {\r
136     return s2;\r
137   }\r
138 \r
139   public void SeqInit(SequenceI s1, SequenceI s2,String type) {\r
140     s1str = extractGaps(".",s1.getSequence());\r
141     s2str = extractGaps(".",s2.getSequence());\r
142     s1str = extractGaps("-",s1str);\r
143     s2str = extractGaps("-",s2str);\r
144     s1str = extractGaps(" ",s1str);\r
145     s2str = extractGaps(" ",s2str);\r
146 \r
147     this.s1 = s1;\r
148     this.s2 = s2;\r
149 \r
150     this.type = type;\r
151 \r
152     if (type.equals("pep")) {\r
153       lookup = ResidueProperties.getBLOSUM62();\r
154       intToStr = pep;\r
155       defInt = 23;\r
156     } else if (type.equals("dna")) {\r
157       lookup = ResidueProperties.getDNA();\r
158       intToStr = dna;\r
159       defInt = 4;\r
160     } else {\r
161       output = output + ("Wrong type = dna or pep only");\r
162       System.exit(0);\r
163     }\r
164 \r
165 \r
166     //System.out.println("lookuip " + rt.freeMemory() + " "+  rt.totalMemory());\r
167     seq1 = new int[s1str.length()];\r
168     //System.out.println("seq1 " + rt.freeMemory() +" "  + rt.totalMemory());\r
169     seq2 = new int[s2str.length()];\r
170     //System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());\r
171     score = new int[s1str.length()][s2str.length()];\r
172     //System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());\r
173     E = new int[s1str.length()][s2str.length()];\r
174     //System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());\r
175     F = new int[s1str.length()][s2str.length()];\r
176     traceback = new int[s1str.length()][s2str.length()];\r
177     //System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());\r
178     seq1 = stringToInt(s1str,type);\r
179     //System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());\r
180     seq2 = stringToInt(s2str,type);\r
181     //System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());\r
182 \r
183     //   long tstart = System.currentTimeMillis();\r
184     //    calcScoreMatrix();\r
185     //long tend = System.currentTimeMillis();\r
186 \r
187     //System.out.println("Time take to calculate score matrix = " + (tend-tstart) + " ms");\r
188 \r
189 \r
190     //   printScoreMatrix(score);\r
191     //System.out.println();\r
192 \r
193     //printScoreMatrix(traceback);\r
194     //System.out.println();\r
195 \r
196     //  printScoreMatrix(E);\r
197     //System.out.println();\r
198 \r
199     ///printScoreMatrix(F);\r
200     //System.out.println();\r
201     // tstart = System.currentTimeMillis();\r
202     //traceAlignment();\r
203     //tend = System.currentTimeMillis();\r
204     //System.out.println("Time take to traceback alignment = " + (tend-tstart) + " ms");\r
205   }\r
206 \r
207   public void traceAlignment() {\r
208 \r
209     // Find the maximum score along the rhs or bottom row\r
210     int max = -9999;\r
211     for (int i = 0; i < seq1.length; i++) {\r
212       if (score[i][seq2.length - 1] > max ) {\r
213         max = score[i][seq2.length - 1];\r
214         maxi = i;\r
215         maxj = seq2.length-1;\r
216       }\r
217     }\r
218     for (int j = 0; j < seq2.length; j++) {\r
219       if (score[seq1.length - 1][j] > max ) {\r
220         max = score[seq1.length - 1][j];\r
221         maxi = seq1.length-1;\r
222         maxj = j;\r
223       }\r
224     }\r
225 \r
226     //  System.out.println(maxi + " " + maxj + " " + score[maxi][maxj]);\r
227 \r
228     int i = maxi;\r
229     int j = maxj;\r
230     int trace;\r
231     maxscore = score[i][j] / 10;\r
232 \r
233     seq1end = maxi+1;\r
234     seq2end = maxj+1;\r
235 \r
236     aseq1 = new int[seq1.length + seq2.length];\r
237     aseq2 = new int[seq1.length + seq2.length];\r
238 \r
239     count = seq1.length + seq2.length - 1;\r
240 \r
241     while (i>0 && j >0) {\r
242 \r
243       if (aseq1[count] != defInt && i >=0) {\r
244         aseq1[count] = seq1[i];\r
245         astr1 = intToStr[seq1[i]] + astr1;\r
246       }\r
247 \r
248       if (aseq2[count] != defInt && j > 0) {\r
249         aseq2[count] = seq2[j];\r
250         astr2 = intToStr[seq2[j]] + astr2;\r
251       }\r
252       trace = findTrace(i,j);\r
253       if (trace == 0) {\r
254         i--;\r
255         j--;\r
256 \r
257       } else  if (trace == 1) {\r
258         j--;\r
259         aseq1[count] = defInt;\r
260         astr1 = "-" + astr1.substring(1);\r
261       } else  if (trace == -1) {\r
262         i--;\r
263         aseq2[count] = defInt;\r
264         astr2 = "-" + astr2.substring(1);\r
265       }\r
266       count--;\r
267     }\r
268 \r
269     seq1start = i+1;\r
270     seq2start = j+1;\r
271 \r
272     if (aseq1[count] != defInt) {\r
273       aseq1[count] = seq1[i];\r
274       astr1 = intToStr[seq1[i]] + astr1;\r
275     }\r
276 \r
277     if (aseq2[count] != defInt) {\r
278       aseq2[count] = seq2[j];\r
279       astr2 = intToStr[seq2[j]] + astr2;\r
280     }\r
281   }\r
282 \r
283   public void printAlignment() {\r
284     // Find the biggest id length for formatting purposes\r
285     int maxid = s1.getName().length();\r
286 \r
287     if (s2.getName().length() > maxid) {\r
288       maxid = s2.getName().length();\r
289     }\r
290 \r
291     int len = 72 - maxid - 1;\r
292     int nochunks = ((aseq1.length - count) / len) + 1;\r
293     pid = 0;\r
294     int overlap = 0;\r
295 \r
296     output = output + ("Score = " + score[maxi][maxj] + "\n");\r
297     output = output + ("Length of alignment = " + (aseq1.length-count) + "\n");\r
298     output = output + ("Sequence ");\r
299     output = output + (new Format("%" + maxid + "s").form(s1.getName()));\r
300     output = output + (" :  " + seq1start + " - " + seq1end + " (Sequence length = " + s1str.length() + ")\n");\r
301     output = output + ("Sequence ");\r
302     output = output + (new Format("%" + maxid + "s").form(s2.getName()));\r
303     output = output + (" :  " + seq2start + " - " + seq2end + " (Sequence length = " + s2str.length() + ")\n\n");\r
304 \r
305     for (int j = 0; j < nochunks; j++) {\r
306       // Print the first aligned sequence\r
307       output = output + (new Format("%" + (maxid) + "s").form(s1.getName()) + " ");\r
308       for (int i = 0; i < len ; i++) {\r
309 \r
310         if ((count + i + j*len) < aseq1.length) {\r
311           output = output + (new Format("%s").form(intToStr[aseq1[count + i + j*len]]));\r
312         }\r
313       }\r
314 \r
315       output = output + ("\n");\r
316       output = output + (new Format("%" + (maxid) + "s").form(" ") + " ");\r
317       // Print out the matching chars\r
318       for (int i = 0; i < len ; i++) {\r
319 \r
320         if ((count + i + j*len) < aseq1.length) {\r
321           if ( intToStr[aseq1[count+i+j*len]].equals(intToStr[aseq2[count+i+j*len]]) && !intToStr[aseq1[count+i+j*len]].equals("-")) {\r
322             pid++;\r
323             output = output + ("|");\r
324           } else if (type.equals("pep")) {\r
325             if (ResidueProperties.getPAM250(intToStr[aseq1[count+i+j*len]],intToStr[aseq2[count+i+j*len]]) > 0) {\r
326               output  = output + (".");\r
327             } else {\r
328               output = output + (" ");\r
329             }\r
330           } else {\r
331             output = output + (" ");\r
332           }\r
333 \r
334         }\r
335       }\r
336       // Now print the second aligned sequence\r
337       output = output + ("\n");\r
338       output = output + (new Format("%" + (maxid) + "s").form(s2.getName()) + " " );\r
339       for (int i = 0; i < len ; i++) {\r
340         if ((count + i + j*len) < aseq1.length) {\r
341           output = output + (new Format("%s").form(intToStr[aseq2[count + i + j*len]]));\r
342         }\r
343       }\r
344       output = output + ("\n\n");\r
345     }\r
346     pid = pid/(float)(aseq1.length-count)*100;\r
347     output = output + (new Format("Percentage ID = %2.2f\n\n").form(pid));\r
348 \r
349   }\r
350 \r
351   public void printScoreMatrix(int[][] mat) {\r
352     int n = seq1.length;\r
353     int m = seq2.length;\r
354 \r
355     for (int i = 0; i < n;i++) {\r
356       // Print the top sequence\r
357       if (i == 0) {\r
358         Format.print(System.out,"%8s",s2str.substring(0,1));\r
359         for (int jj = 1;jj < m; jj++) {\r
360           Format.print(System.out,"%5s",s2str.substring(jj,jj+1));\r
361         }\r
362         System.out.println();\r
363       }\r
364 \r
365       for (int j = 0;j < m; j++) {\r
366         if (j == 0) {\r
367           Format.print(System.out,"%3s",s1str.substring(i,i+1));\r
368         }\r
369         Format.print(System.out,"%3d ",mat[i][j]/10);\r
370       }\r
371       System.out.println();\r
372     }\r
373   }\r
374 \r
375   public int findTrace(int i,int j) {\r
376     int t = 0;\r
377     int max = score[i-1][j-1] + lookup[seq1[i]][seq2[j]] * 10;\r
378 \r
379     if (F[i][j] > max) {\r
380       max = F[i][j];\r
381       t = -1;\r
382     } else if (F[i][j] == max) {\r
383       if (prev == -1) {\r
384         max = F[i][j];\r
385         t = -1;\r
386       }\r
387     }\r
388     if (E[i][j] >= max) {\r
389       max = E[i][j];\r
390       t = 1;\r
391     } else if (E[i][j] == max) {\r
392       if (prev == 1) {\r
393         max = E[i][j];\r
394         t = 1;\r
395       }\r
396     }\r
397     prev = t;\r
398     return t;\r
399   }\r
400 \r
401   public void calcScoreMatrix() {\r
402 \r
403 \r
404     int n = seq1.length;\r
405     int m = seq2.length;\r
406 \r
407 \r
408     // top left hand element\r
409     score[0][0] =  lookup[seq1[0]][seq2[0]] * 10;\r
410     E[0][0] = -gapExtend;\r
411     F[0][0] = 0;\r
412 \r
413     // Calculate the top row first\r
414     for (int j=1; j < m; j++) {\r
415       // What should these values be? 0 maybe\r
416       E[0][j] = max(score[0][j-1] - gapOpen,E[0][j-1] - gapExtend);\r
417       F[0][j] = -gapExtend;\r
418 \r
419       score[0][j] = max( lookup[seq1[0]][seq2[j]] * 10 ,-gapOpen,-gapExtend);\r
420 \r
421       traceback[0][j] = 1;\r
422     }\r
423 \r
424     // Now do the left hand column\r
425     for (int i=1; i < n; i++) {\r
426       E[i][0] =  -gapOpen;\r
427       F[i][0] =  max(score[i-1][0]-gapOpen,F[i-1][0]-gapExtend);\r
428 \r
429       score[i][0] = max( lookup[seq1[i]][seq2[0]] * 10 ,E[i][0],F[i][0]);\r
430       traceback[i][0] = -1;\r
431     }\r
432 \r
433     // Now do all the other rows\r
434     for (int i = 1; i < n; i++) {\r
435       for (int j = 1; j < m; j++) {\r
436 \r
437         E[i][j] = max(score[i][j-1] - gapOpen,  E[i][j-1] - gapExtend);\r
438         F[i][j] = max(score[i-1][j] - gapOpen,  F[i-1][j] - gapExtend);\r
439 \r
440         score[i][j] = max(score[i-1][j-1] + lookup[seq1[i]][seq2[j]]*10,\r
441                           E[i][j],\r
442                           F[i][j]);\r
443         traceback[i][j] = findTrace(i,j);\r
444       }\r
445     }\r
446 \r
447   }\r
448   public static String extractChars(String chars, String seq) {\r
449     String out = seq;\r
450     for (int i=0; i < chars.length(); i++) {\r
451       String gap = chars.substring(i,i+1);\r
452       out = extractGaps(gap,out);\r
453     }\r
454     return out;\r
455   }\r
456   public static String extractGaps(String gapChar, String seq) {\r
457     StringTokenizer str = new StringTokenizer(seq,gapChar);\r
458     String newString = "";\r
459 \r
460     while (str.hasMoreTokens()) {\r
461       newString  = newString + str.nextToken();\r
462     }\r
463     return newString;\r
464   }\r
465 \r
466 \r
467   public int max(int i1, int i2, int i3) {\r
468     int max = i1;\r
469     if (i2 > i1)  {\r
470       max = i2;\r
471     }\r
472     if (i3 > max) {\r
473       max = i3;\r
474     }\r
475     return max;\r
476   }\r
477 \r
478   public int max(int i1, int i2) {\r
479     int max = i1;\r
480     if (i2 > i1)  {\r
481       max = i2;\r
482     }\r
483     return max;\r
484   }\r
485 \r
486   public int[] stringToInt(String s,String type) {\r
487     int[] seq1 = new int[s.length()];\r
488 \r
489     for (int i = 0;i < s.length(); i++) {\r
490       String ss = s.substring(i,i+1).toUpperCase();\r
491       try {\r
492         if (type.equals("pep")) {\r
493           seq1[i] = ((Integer)ResidueProperties.aaHash.get(ss)).intValue();\r
494         } else if (type.equals("dna")) {\r
495           seq1[i] = ((Integer)dnaHash.get(ss)).intValue();\r
496         }\r
497         if (seq1[i] > 23) {\r
498           seq1[i] = 23;\r
499         }\r
500       } catch (Exception e) {\r
501         if (type.equals("dna")) {\r
502           seq1[i] = 4;\r
503         } else {\r
504           seq1[i] = 23;\r
505         }\r
506       }\r
507     }\r
508     return seq1;\r
509   }\r
510 \r
511   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,int psize) {\r
512 \r
513     int max = -1000;\r
514     int min = 1000;\r
515 \r
516     for (int i=0; i < n; i++) {\r
517       for (int j=0; j < m; j++) {\r
518         if (mat[i][j] >= max) {\r
519           max = mat[i][j];\r
520         }\r
521         if (mat[i][j] <= min) {\r
522           min = mat[i][j];\r
523         }\r
524       }\r
525     }\r
526     System.out.println(max + " " + min);\r
527     for (int i=0; i < n; i++) {\r
528       for (int j=0; j < m; j++) {\r
529         int x = psize*i;\r
530         int y = psize*j;\r
531         //      System.out.println(mat[i][j]);\r
532         float score = (float)(mat[i][j] - min)/(float)(max-min);\r
533         g.setColor(new Color(score,0,0));\r
534         g.fillRect(x,y,psize,psize);\r
535         //      System.out.println(x + " " + y + " " + score);\r
536       }\r
537 \r
538     }\r
539   }\r
540 \r
541 }\r