added tree node property for sequence placeholders.
[jalview.git] / src / jalview / analysis / NJTree.java
1 package jalview.analysis;\r
2 \r
3 import jalview.datamodel.*;\r
4 import jalview.util.*;\r
5 import jalview.schemes.ResidueProperties;\r
6 import java.util.*;\r
7 \r
8 import jalview.io.NewickFile;\r
9 \r
10 public class NJTree {\r
11 \r
12   Vector cluster;\r
13   SequenceI[] sequence;\r
14 \r
15   int done[];\r
16   int noseqs;\r
17   int noClus;\r
18 \r
19   float distance[][];\r
20 \r
21   int mini;\r
22   int minj;\r
23   float ri;\r
24   float rj;\r
25 \r
26   Vector groups = new Vector();\r
27   SequenceNode maxdist;\r
28   SequenceNode top;\r
29 \r
30   float maxDistValue;\r
31   float maxheight;\r
32 \r
33   int ycount;\r
34 \r
35   Vector node;\r
36 \r
37   String type;\r
38   String pwtype;\r
39 \r
40   Object found = null;\r
41   Object leaves = null;\r
42 \r
43   int start;\r
44   int end;\r
45 \r
46   public NJTree(SequenceNode node) {\r
47     top = node;\r
48     maxheight = findHeight(top);\r
49   }\r
50 \r
51   public String toString()\r
52   {\r
53     jalview.io.NewickFile fout = new jalview.io.NewickFile(getTopNode());\r
54     return fout.print(false,true); // distances only\r
55   }\r
56 \r
57   public NJTree(SequenceI[] seqs, NewickFile treefile) {\r
58     top = treefile.getTree();\r
59     maxheight = findHeight(top);\r
60     SequenceIdMatcher algnIds = new SequenceIdMatcher(seqs);\r
61 \r
62     Vector leaves = new Vector();\r
63     findLeaves(top, leaves);\r
64 \r
65     int i = 0;\r
66     int namesleft = seqs.length;\r
67 \r
68     SequenceNode j;\r
69     SequenceI nam;\r
70     String realnam;\r
71     while (i < leaves.size())\r
72     {\r
73       j = (SequenceNode) leaves.elementAt(i++);\r
74       realnam = j.getName();\r
75       nam = null;\r
76       if (namesleft>-1)\r
77         nam = algnIds.findIdMatch(realnam);\r
78       if (nam != null) {\r
79         j.setElement(nam);\r
80         namesleft--;\r
81       } else {\r
82         j.setElement(new Sequence(realnam, "THISISAPLACEHLDER"));\r
83         j.setPlaceholder(true);\r
84 \r
85       }\r
86     }\r
87   }\r
88 \r
89   public NJTree(SequenceI[] sequence,int start, int end) {\r
90     this(sequence,"NJ","BL",start,end);\r
91   }\r
92 \r
93   public NJTree(SequenceI[] sequence,String type,String pwtype,int start, int end ) {\r
94 \r
95     this.sequence = sequence;\r
96     this.node     = new Vector();\r
97     this.type     = type;\r
98     this.pwtype   = pwtype;\r
99     this.start    = start;\r
100     this.end      = end;\r
101 \r
102     if (!(type.equals("NJ"))) {\r
103       type = "AV";\r
104     }\r
105 \r
106     if (!(pwtype.equals("PID"))) {\r
107       type = "BL";\r
108     }\r
109 \r
110     int i=0;\r
111 \r
112     done = new int[sequence.length];\r
113 \r
114 \r
115     while (i < sequence.length  && sequence[i] != null) {\r
116       done[i] = 0;\r
117       i++;\r
118     }\r
119 \r
120     noseqs = i++;\r
121 \r
122     distance = findDistances();\r
123 \r
124     makeLeaves();\r
125 \r
126     noClus = cluster.size();\r
127 \r
128     cluster();\r
129 \r
130   }\r
131 \r
132 \r
133   public void cluster() {\r
134 \r
135     while (noClus > 2) {\r
136       if (type.equals("NJ")) {\r
137         float mind = findMinNJDistance();\r
138       } else {\r
139         float mind = findMinDistance();\r
140       }\r
141 \r
142       Cluster c = joinClusters(mini,minj);\r
143 \r
144 \r
145       done[minj] = 1;\r
146 \r
147       cluster.setElementAt(null,minj);\r
148       cluster.setElementAt(c,mini);\r
149 \r
150       noClus--;\r
151     }\r
152 \r
153     boolean onefound = false;\r
154 \r
155     int one = -1;\r
156     int two = -1;\r
157 \r
158     for (int i=0; i < noseqs; i++) {\r
159       if (done[i] != 1) {\r
160         if (onefound == false) {\r
161           two = i;\r
162           onefound = true;\r
163         } else {\r
164           one = i;\r
165         }\r
166       }\r
167     }\r
168 \r
169     Cluster c = joinClusters(one,two);\r
170     top = (SequenceNode)(node.elementAt(one));\r
171 \r
172     reCount(top);\r
173     findHeight(top);\r
174     findMaxDist(top);\r
175 \r
176   }\r
177 \r
178   public Cluster joinClusters(int i, int j) {\r
179 \r
180     float dist = distance[i][j];\r
181 \r
182     int noi = ((Cluster)cluster.elementAt(i)).value.length;\r
183     int noj = ((Cluster)cluster.elementAt(j)).value.length;\r
184 \r
185     int[] value = new int[noi + noj];\r
186 \r
187     for (int ii = 0; ii < noi;ii++) {\r
188       value[ii] =  ((Cluster)cluster.elementAt(i)).value[ii];\r
189     }\r
190 \r
191     for (int ii = noi; ii < noi+ noj;ii++) {\r
192       value[ii] =  ((Cluster)cluster.elementAt(j)).value[ii-noi];\r
193     }\r
194 \r
195     Cluster c = new Cluster(value);\r
196 \r
197     ri = findr(i,j);\r
198     rj = findr(j,i);\r
199 \r
200     if (type.equals("NJ")) {\r
201       findClusterNJDistance(i,j);\r
202     } else {\r
203       findClusterDistance(i,j);\r
204     }\r
205 \r
206     SequenceNode sn = new SequenceNode();\r
207 \r
208     sn.setLeft((SequenceNode)(node.elementAt(i)));\r
209     sn.setRight((SequenceNode)(node.elementAt(j)));\r
210 \r
211     SequenceNode tmpi = (SequenceNode)(node.elementAt(i));\r
212     SequenceNode tmpj = (SequenceNode)(node.elementAt(j));\r
213 \r
214     if (type.equals("NJ")) {\r
215       findNewNJDistances(tmpi,tmpj,dist);\r
216     } else {\r
217       findNewDistances(tmpi,tmpj,dist);\r
218     }\r
219 \r
220     tmpi.setParent(sn);\r
221     tmpj.setParent(sn);\r
222 \r
223     node.setElementAt(sn,i);\r
224     return c;\r
225   }\r
226 \r
227   public void findNewNJDistances(SequenceNode tmpi, SequenceNode tmpj, float dist) {\r
228 \r
229     float ih = 0;\r
230     float jh = 0;\r
231 \r
232     SequenceNode sni = tmpi;\r
233     SequenceNode snj = tmpj;\r
234 \r
235     tmpi.dist = (dist + ri - rj)/2;\r
236     tmpj.dist = (dist - tmpi.dist);\r
237 \r
238     if (tmpi.dist < 0) {\r
239       tmpi.dist = 0;\r
240     }\r
241     if (tmpj.dist < 0) {\r
242       tmpj.dist = 0;\r
243     }\r
244   }\r
245 \r
246   public void findNewDistances(SequenceNode tmpi,SequenceNode tmpj,float dist) {\r
247 \r
248     float ih = 0;\r
249     float jh = 0;\r
250 \r
251     SequenceNode sni = tmpi;\r
252     SequenceNode snj = tmpj;\r
253 \r
254     while (sni != null) {\r
255       ih = ih + sni.dist;\r
256       sni = (SequenceNode)sni.left();\r
257     }\r
258 \r
259     while (snj != null) {\r
260       jh = jh + snj.dist;\r
261       snj = (SequenceNode)snj.left();\r
262     }\r
263 \r
264     tmpi.dist = (dist/2 - ih);\r
265     tmpj.dist = (dist/2 - jh);\r
266   }\r
267 \r
268 \r
269 \r
270   public void findClusterDistance(int i, int j) {\r
271 \r
272     int noi = ((Cluster)cluster.elementAt(i)).value.length;\r
273     int noj = ((Cluster)cluster.elementAt(j)).value.length;\r
274 \r
275     // New distances from cluster to others\r
276     float[] newdist = new float[noseqs];\r
277 \r
278     for (int l = 0; l < noseqs; l++) {\r
279       if ( l != i && l != j) {\r
280         newdist[l] = (distance[i][l] * noi + distance[j][l] * noj)/(noi + noj);\r
281       } else {\r
282         newdist[l] = 0;\r
283       }\r
284     }\r
285 \r
286     for (int ii=0; ii < noseqs;ii++) {\r
287       distance[i][ii] = newdist[ii];\r
288       distance[ii][i] = newdist[ii];\r
289     }\r
290   }\r
291 \r
292   public void findClusterNJDistance(int i, int j) {\r
293 \r
294     int noi = ((Cluster)cluster.elementAt(i)).value.length;\r
295     int noj = ((Cluster)cluster.elementAt(j)).value.length;\r
296 \r
297     // New distances from cluster to others\r
298     float[] newdist = new float[noseqs];\r
299 \r
300     for (int l = 0; l < noseqs; l++) {\r
301       if ( l != i && l != j) {\r
302         newdist[l] = (distance[i][l] + distance[j][l] - distance[i][j])/2;\r
303       } else {\r
304         newdist[l] = 0;\r
305       }\r
306     }\r
307 \r
308     for (int ii=0; ii < noseqs;ii++) {\r
309       distance[i][ii] = newdist[ii];\r
310       distance[ii][i] = newdist[ii];\r
311     }\r
312   }\r
313 \r
314   public float findr(int i, int j) {\r
315 \r
316     float tmp = 1;\r
317     for (int k=0; k < noseqs;k++) {\r
318       if (k!= i && k!= j && done[k] != 1) {\r
319         tmp = tmp + distance[i][k];\r
320       }\r
321     }\r
322 \r
323     if (noClus > 2) {\r
324       tmp = tmp/(noClus - 2);\r
325     }\r
326 \r
327     return tmp;\r
328   }\r
329 \r
330   public float findMinNJDistance() {\r
331 \r
332     float min = 100000;\r
333 \r
334     for (int i=0; i < noseqs-1; i++) {\r
335       for (int j=i+1;j < noseqs;j++) {\r
336         if (done[i] != 1 && done[j] != 1) {\r
337           float tmp = distance[i][j] - (findr(i,j) + findr(j,i));\r
338           if (tmp < min) {\r
339 \r
340             mini = i;\r
341             minj = j;\r
342 \r
343             min = tmp;\r
344 \r
345           }\r
346         }\r
347       }\r
348     }\r
349     return min;\r
350   }\r
351 \r
352   public float findMinDistance() {\r
353 \r
354     float min = 100000;\r
355 \r
356     for (int i=0; i < noseqs-1;i++) {\r
357       for (int j = i+1; j < noseqs;j++) {\r
358         if (done[i] != 1 && done[j] != 1) {\r
359           if (distance[i][j] < min) {\r
360             mini = i;\r
361             minj = j;\r
362 \r
363             min = distance[i][j];\r
364           }\r
365         }\r
366       }\r
367     }\r
368     return min;\r
369   }\r
370 \r
371   public float[][] findDistances() {\r
372 \r
373     float[][] distance = new float[noseqs][noseqs];\r
374     if (pwtype.equals("PID")) {\r
375       for (int i = 0; i < noseqs-1; i++) {\r
376         for (int j = i; j < noseqs; j++) {\r
377           if (j==i) {\r
378             distance[i][i] = 0;\r
379           } else {\r
380             distance[i][j] = 100-Comparison.PID(sequence[i], sequence[j]);\r
381             distance[j][i] = distance[i][j];\r
382           }\r
383         }\r
384       }\r
385     } else if (pwtype.equals("BL")) {\r
386       int   maxscore = 0;\r
387 \r
388       for (int i = 0; i < noseqs-1; i++) {\r
389         for (int j = i; j < noseqs; j++) {\r
390           int score = 0;\r
391           for (int k=0; k < sequence[i].getLength(); k++) {\r
392             try{\r
393               score +=\r
394                   ResidueProperties.getBLOSUM62(sequence[i].getSequence(k,\r
395                   k + 1),\r
396                                                 sequence[j].getSequence(k,\r
397                   k + 1));\r
398             }catch(Exception ex){System.out.println("err creating BLOSUM62 tree");}\r
399           }\r
400           distance[i][j] = (float)score;\r
401           if (score > maxscore) {\r
402             maxscore = score;\r
403           }\r
404         }\r
405       }\r
406       for (int i = 0; i < noseqs-1; i++) {\r
407         for (int j = i; j < noseqs; j++) {\r
408           distance[i][j] =  (float)maxscore - distance[i][j];\r
409           distance[j][i] = distance[i][j];\r
410         }\r
411       }\r
412     } else if (pwtype.equals("SW")) {\r
413       float max = -1;\r
414       for (int i = 0; i < noseqs-1; i++) {\r
415         for (int j = i; j < noseqs; j++) {\r
416           AlignSeq as = new AlignSeq(sequence[i],sequence[j],"pep");\r
417           as.calcScoreMatrix();\r
418           as.traceAlignment();\r
419           as.printAlignment();\r
420           distance[i][j] = (float)as.maxscore;\r
421           if (max < distance[i][j]) {\r
422             max = distance[i][j];\r
423           }\r
424         }\r
425       }\r
426       for (int i = 0; i < noseqs-1; i++) {\r
427         for (int j = i; j < noseqs; j++) {\r
428           distance[i][j] =  max - distance[i][j];\r
429           distance[j][i] = distance[i][j];\r
430         }\r
431       }\r
432     }\r
433 \r
434     return distance;\r
435   }\r
436 \r
437   public void makeLeaves() {\r
438     cluster = new Vector();\r
439 \r
440     for (int i=0; i < noseqs; i++) {\r
441       SequenceNode sn = new SequenceNode();\r
442 \r
443       sn.setElement(sequence[i]);\r
444       sn.setName(sequence[i].getName());\r
445       node.addElement(sn);\r
446 \r
447       int[] value = new int[1];\r
448       value[0] = i;\r
449 \r
450       Cluster c = new Cluster(value);\r
451       cluster.addElement(c);\r
452     }\r
453   }\r
454 \r
455   public Vector findLeaves(SequenceNode node, Vector leaves) {\r
456     if (node == null) {\r
457       return leaves;\r
458     }\r
459 \r
460     if (node.left() == null && node.right() == null) {\r
461       leaves.addElement(node);\r
462       return leaves;\r
463     } else {\r
464       findLeaves((SequenceNode)node.left(),leaves);\r
465       findLeaves((SequenceNode)node.right(),leaves);\r
466     }\r
467     return leaves;\r
468   }\r
469 \r
470   public Object findLeaf(SequenceNode node, int count) {\r
471     found = _findLeaf(node,count);\r
472 \r
473     return found;\r
474   }\r
475   public Object _findLeaf(SequenceNode node,int count) {\r
476     if (node == null) {\r
477       return null;\r
478     }\r
479     if (node.ycount == count) {\r
480       found = node.element();\r
481       return found;\r
482     } else {\r
483       _findLeaf((SequenceNode)node.left(),count);\r
484       _findLeaf((SequenceNode)node.right(),count);\r
485     }\r
486 \r
487     return found;\r
488   }\r
489 \r
490   public void printNode(SequenceNode node) {\r
491     if (node == null) {\r
492       return;\r
493     }\r
494     if (node.left() == null && node.right() == null) {\r
495       System.out.println("Leaf = " + ((SequenceI)node.element()).getName());\r
496       System.out.println("Dist " + ((SequenceNode)node).dist);\r
497       System.out.println("Boot " + node.getBootstrap());\r
498     } else {\r
499       System.out.println("Dist " + ((SequenceNode)node).dist);\r
500       printNode((SequenceNode)node.left());\r
501       printNode((SequenceNode)node.right());\r
502     }\r
503   }\r
504   public void findMaxDist(SequenceNode node) {\r
505     if (node == null) {\r
506       return;\r
507     }\r
508     if (node.left() == null && node.right() == null) {\r
509 \r
510       float dist = ((SequenceNode)node).dist;\r
511       if (dist > maxDistValue) {\r
512           maxdist      = (SequenceNode)node;\r
513           maxDistValue = dist;\r
514       }\r
515     } else {\r
516       findMaxDist((SequenceNode)node.left());\r
517       findMaxDist((SequenceNode)node.right());\r
518     }\r
519   }\r
520     public Vector getGroups() {\r
521         return groups;\r
522     }\r
523     public float getMaxHeight() {\r
524         return maxheight;\r
525     }\r
526   public void  groupNodes(SequenceNode node, float threshold) {\r
527     if (node == null) {\r
528       return;\r
529     }\r
530 \r
531     if (node.height/maxheight > threshold) {\r
532       groups.addElement(node);\r
533     } else {\r
534       groupNodes((SequenceNode)node.left(),threshold);\r
535       groupNodes((SequenceNode)node.right(),threshold);\r
536     }\r
537   }\r
538 \r
539   public float findHeight(SequenceNode node) {\r
540 \r
541     if (node == null) {\r
542       return maxheight;\r
543     }\r
544 \r
545     if (node.left() == null && node.right() == null) {\r
546       node.height = ((SequenceNode)node.parent()).height + node.dist;\r
547 \r
548       if (node.height > maxheight) {\r
549         return node.height;\r
550       } else {\r
551         return maxheight;\r
552       }\r
553     } else {\r
554       if (node.parent() != null) {\r
555         node.height = ((SequenceNode)node.parent()).height + node.dist;\r
556       } else {\r
557         maxheight = 0;\r
558         node.height = (float)0.0;\r
559       }\r
560 \r
561       maxheight = findHeight((SequenceNode)(node.left()));\r
562       maxheight = findHeight((SequenceNode)(node.right()));\r
563     }\r
564     return maxheight;\r
565   }\r
566   public SequenceNode reRoot() {\r
567     if (maxdist != null) {\r
568       ycount = 0;\r
569       float tmpdist = maxdist.dist;\r
570 \r
571       // New top\r
572       SequenceNode sn = new SequenceNode();\r
573       sn.setParent(null);\r
574 \r
575       // New right hand of top\r
576       SequenceNode snr = (SequenceNode)maxdist.parent();\r
577       changeDirection(snr,maxdist);\r
578       System.out.println("Printing reversed tree");\r
579       printN(snr);\r
580       snr.dist = tmpdist/2;\r
581       maxdist.dist = tmpdist/2;\r
582 \r
583       snr.setParent(sn);\r
584       maxdist.setParent(sn);\r
585 \r
586       sn.setRight(snr);\r
587       sn.setLeft(maxdist);\r
588 \r
589       top = sn;\r
590 \r
591       ycount = 0;\r
592       reCount(top);\r
593       findHeight(top);\r
594 \r
595     }\r
596     return top;\r
597   }\r
598   public static void printN(SequenceNode node) {\r
599     if (node == null) {\r
600       return;\r
601     }\r
602 \r
603     if (node.left() != null && node.right() != null) {\r
604       printN((SequenceNode)node.left());\r
605       printN((SequenceNode)node.right());\r
606     } else {\r
607       System.out.println(" name = " + ((SequenceI)node.element()).getName());\r
608     }\r
609     System.out.println(" dist = " + ((SequenceNode)node).dist + " " + ((SequenceNode)node).count + " " + ((SequenceNode)node).height);\r
610   }\r
611 \r
612     public void reCount(SequenceNode node) {\r
613         ycount = 0;\r
614         _reCount(node);\r
615     }\r
616   public void _reCount(SequenceNode node) {\r
617     if (node == null) {\r
618       return;\r
619     }\r
620 \r
621     if (node.left() != null && node.right() != null) {\r
622       _reCount((SequenceNode)node.left());\r
623       _reCount((SequenceNode)node.right());\r
624 \r
625       SequenceNode l = (SequenceNode)node.left();\r
626       SequenceNode r = (SequenceNode)node.right();\r
627 \r
628       ((SequenceNode)node).count  = l.count + r.count;\r
629       ((SequenceNode)node).ycount = (l.ycount + r.ycount)/2;\r
630 \r
631     } else {\r
632       ((SequenceNode)node).count = 1;\r
633       ((SequenceNode)node).ycount = ycount++;\r
634     }\r
635 \r
636   }\r
637     public void swapNodes(SequenceNode node) {\r
638         if (node == null) {\r
639             return;\r
640         }\r
641         SequenceNode tmp = (SequenceNode)node.left();\r
642 \r
643         node.setLeft(node.right());\r
644         node.setRight(tmp);\r
645     }\r
646   public void changeDirection(SequenceNode node, SequenceNode dir) {\r
647     if (node == null) {\r
648       return;\r
649     }\r
650     if (node.parent() != top) {\r
651       changeDirection((SequenceNode)node.parent(), node);\r
652 \r
653       SequenceNode tmp = (SequenceNode)node.parent();\r
654 \r
655       if (dir == node.left()) {\r
656         node.setParent(dir);\r
657         node.setLeft(tmp);\r
658       } else if (dir == node.right()) {\r
659         node.setParent(dir);\r
660         node.setRight(tmp);\r
661       }\r
662 \r
663     } else {\r
664       if (dir == node.left()) {\r
665         node.setParent(node.left());\r
666 \r
667         if (top.left() == node) {\r
668           node.setRight(top.right());\r
669         } else {\r
670           node.setRight(top.left());\r
671         }\r
672       } else {\r
673         node.setParent(node.right());\r
674 \r
675         if (top.left() == node) {\r
676           node.setLeft(top.right());\r
677         } else {\r
678           node.setLeft(top.left());\r
679         }\r
680       }\r
681     }\r
682   }\r
683     public void setMaxDist(SequenceNode node) {\r
684         this.maxdist = maxdist;\r
685     }\r
686     public SequenceNode getMaxDist() {\r
687         return maxdist;\r
688     }\r
689     public SequenceNode getTopNode() {\r
690         return top;\r
691     }\r
692 \r
693 }\r
694 \r
695 \r
696 \r
697 class Cluster {\r
698 \r
699   int[] value;\r
700 \r
701   public Cluster(int[] value) {\r
702     this.value = value;\r
703   }\r
704 \r
705 }\r
706 \r