source: branches/ralph/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/pal/tree/UPGMAAligningTree.java @ 1585

Last change on this file since 1585 was 1585, checked in by rkrimmel, 10 years ago

Multiple Alignment first version
File size: 8.5 KB
Line 
1// UPGMATree.java
2//
3// (c) 1999-2001 PAL Development Core Team
4//
5// This package may be distributed under the
6// terms of the Lesser GNU General Public License (LGPL)
7
8// Known bugs and limitations:
9// - computational complexity O(numSeqs^3)
10//   (this could be brought down to O(numSeqs^2)
11//   but this needs more clever programming ...)
12
13
14package de.ugoe.cs.autoquest.tasktrees.alignment.pal.tree;
15
16import java.util.ArrayList;
17import java.util.Iterator;
18import java.util.logging.Level;
19
20import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.NumberSequence;
21import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.SmithWatermanRepeated;
22import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.BinaryAlignmentStorage;
23import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.ObjectDistanceSubstitionMatrix;
24import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.UPGMAMatrix;
25import de.ugoe.cs.autoquest.tasktrees.alignment.pal.misc.Identifier;
26import de.ugoe.cs.util.console.Console;
27
28
29/**
30 * constructs a UPGMA tree from pairwise distances
31 *
32 * @version $Id: UPGMATree.java,v 1.9 2001/07/13 14:39:13 korbinian Exp $
33 *
34 * @author Korbinian Strimmer
35 * @author Alexei Drummond
36 */
37public class UPGMAAligningTree extends SimpleTree
38{
39        //
40        // Public stuff
41        //     
42
43        /**
44         * constructor UPGMA tree
45         * @param numberseqs
46         *
47         * @param m distance matrix
48         */
49        public UPGMAAligningTree(ArrayList<NumberSequence> numberseqs, BinaryAlignmentStorage alignments, ObjectDistanceSubstitionMatrix submat)
50        {
51                if (alignments.getDistanceMatrix().size() < 2)
52                {
53                        new IllegalArgumentException("LESS THAN 2 TAXA IN DISTANCE MATRIX");
54                }
55       
56                this.numberseqs = numberseqs;
57                this.alignments = alignments;
58                this.submat = submat;
59                init(alignments.getDistanceMatrix());
60
61                while (true)
62                {
63                        findNextPair();
64                        newBranchLengths();
65                       
66                        if (numClusters == 2)
67                        {
68                                break;
69                        }
70                       
71                        newCluster();
72                }
73               
74                finish();
75                createNodeList();
76        }
77
78
79        //
80        // Private stuff
81        //
82        private ArrayList<NumberSequence> numberseqs;
83        private BinaryAlignmentStorage alignments;
84        private ObjectDistanceSubstitionMatrix submat;
85        private int numClusters;
86        private int besti, abi;
87        private int bestj, abj;
88        private int[] alias;
89        private double[][] distance;
90
91        private double[] height;
92        private int[] oc;
93
94        private double getDist(int a, int b)
95        {
96                return distance[alias[a]][alias[b]];
97        }
98       
99        private void init(UPGMAMatrix m)
100        {
101                numClusters = m.size();
102
103                distance = new double[numClusters][numClusters];
104                for (int i = 0; i < numClusters; i++)
105                {
106                        for (int j = 0; j < numClusters; j++)
107                        {
108                                distance[i][j] = m.get(i,j);
109                        }
110                }
111
112                for (int i = 0; i < numClusters; i++)
113                {
114                        Node tmp = NodeFactory.createNode();
115                        tmp.setIdentifier(new Identifier(Integer.toString(i)));
116                        tmp.setNumber(i);
117                        tmp.addSequence(numberseqs.get(i));
118                        getRoot().addChild(tmp);
119                }
120               
121                alias = new int[numClusters];
122                for (int i = 0; i < numClusters; i++)
123                {
124                        alias[i] = i;
125                }
126                               
127                height = new double[numClusters];
128                oc = new int[numClusters];
129                for (int i = 0; i < numClusters; i++)
130                {
131                        height[i] = 0.0;
132                        oc[i] = 1;
133                }
134        }
135
136        private void finish()
137        {
138                distance = null;               
139        }
140
141        private void findNextPair()
142        {
143                besti = 0;
144                bestj = 1;
145                double dmin = getDist(0, 1);
146                for (int i = 0; i < numClusters-1; i++)
147                {
148                        for (int j = i+1; j < numClusters; j++)
149                        {
150                                if (getDist(i, j) < dmin)
151                                {
152                                        dmin = getDist(i, j);
153                                        besti = i;
154                                        bestj = j;
155                                }
156                        }
157                }
158                abi = alias[besti];
159                abj = alias[bestj];
160                //System.out.println("Found best pair: " + abi + "/" +abj + " - "+ besti+ "/"+bestj +" with distance " + dmin);
161               
162        }
163
164        private void newBranchLengths()
165        {
166                double dij = getDist(besti, bestj);
167               
168                getRoot().getChild(besti).setBranchLength(dij/2.0-height[abi]);
169                getRoot().getChild(bestj).setBranchLength(dij/2.0-height[abj]);
170        }
171
172        private void newCluster()
173        {
174                // Update distances
175                for (int k = 0; k < numClusters; k++)
176                {
177                        if (k != besti && k != bestj)
178                        {
179                                int ak = alias[k];
180                                double updated = updatedDistance(besti,bestj,k);
181                                distance[ak][abi] = distance[abi][ak] = updated;
182                        }
183                }
184                distance[abi][abi] = 0.0;
185
186                // Update UPGMA variables
187                height[abi] = getDist(besti, bestj)/2.0;
188                oc[abi] += oc[abj];
189               
190                // Index besti now represent the new cluster
191                Node newNode = getRoot().joinChildren(besti, bestj);
192               
193                if(newNode instanceof FengDoolittleNode) {
194                        newNode.setSequences(alignSequences(newNode));
195                }
196               
197                // Update alias
198                for (int i = bestj; i < numClusters-1; i++)
199                {
200                        alias[i] = alias[i+1];
201                }
202               
203                numClusters--;
204        }
205       
206       
207        public ArrayList<NumberSequence> alignSequences(Node parent) {
208                ArrayList<NumberSequence> alignment = new ArrayList<NumberSequence>();
209                if(parent.getChildCount()<3) {
210                       
211                        Node node1 = parent.getChild(0);
212                        Node node2 = parent.getChild(1);
213                       
214                        int seqCount1 = node1.getSequences().size();
215                        int seqCount2 = node2.getSequences().size();
216                        Console.traceln(Level.INFO,"Merging node " + node1.getIdentifier() + " with " + node2.getIdentifier());
217                        //Console.println("SeqCount1: " + seqCount1 + " seqCount2 " + seqCount2);
218                        //Align 2 sequences
219                        if(seqCount1 == 1 && seqCount2 == 1) {
220                                alignment = (alignments.get(node1.getNumber(), node2.getNumber())).getAlignment();
221                        }
222                        //Align a sequence to a group
223                        else if( seqCount1 > 1 && seqCount2 == 1) {
224                                alignment.addAll(node1.getSequences());
225                               
226                                BinaryAlignmentStorage tempStorage = new BinaryAlignmentStorage(seqCount1,seqCount2);
227                                double maxScore = 0.0;
228                                int maxIndex = 0;
229                                for(int i=0;i<seqCount1;i++) {
230                                        tempStorage.set(i, 1, new SmithWatermanRepeated(node1.getSequence(i).getSequence(), node2.getSequence(0).getSequence() , submat, 5));
231                                        if(maxScore < tempStorage.get(i, 1).getAlignmentScore()) {
232                                                maxScore = tempStorage.get(i, 1).getAlignmentScore();
233                                                maxIndex = i;
234                                        }
235                                }
236                                alignment.add(tempStorage.get(maxIndex, 1).getAlignment().get(1));
237                        }
238                        //Align a sequence to a group
239                        else if(seqCount1 == 1 && seqCount2 > 1) {
240                                alignment.addAll(node2.getSequences());
241                                BinaryAlignmentStorage tempStorage = new BinaryAlignmentStorage(seqCount1,seqCount2);
242                                double maxScore = 0.0;
243                                int maxIndex = 0;
244                                for(int i=0;i<seqCount2;i++) {
245                                        tempStorage.set(1, i, new SmithWatermanRepeated(node2.getSequence(i).getSequence(), node1.getSequence(0).getSequence() , submat, 5));
246                                        if(maxScore < tempStorage.get(1, i).getAlignmentScore()) {
247                                                maxScore = tempStorage.get(1, i).getAlignmentScore();
248                                                maxIndex = i;
249                                        }
250                                }
251                                alignment.add(tempStorage.get(1,maxIndex).getAlignment().get(1));
252                        }
253                        //Align 2 groups
254                        else if((seqCount1 > 1) && (seqCount2 > 1)){
255                                        BinaryAlignmentStorage tempStorage1 = new BinaryAlignmentStorage(seqCount2,1);
256                                        BinaryAlignmentStorage tempStorage2 = new BinaryAlignmentStorage(seqCount1,1);
257                                        double maxScore1 = 0.0;
258                                        double maxScore2 = 0.0;
259                                        int maxIndex1 = 0;
260                                        int maxIndex2 = 0;
261                                        for(int i=0;i<seqCount1;i++) {
262                                                for(int j=0;j<seqCount2;j++) {
263                                                        tempStorage1.set(j, 0, new SmithWatermanRepeated(node1.getSequence(i).getSequence(), node2.getSequence(j).getSequence() , submat, 5));
264                                                        if(maxScore1 < tempStorage1.get(j, 0).getAlignmentScore()) {
265                                                                maxScore1 = tempStorage1.get(j, 0).getAlignmentScore();
266                                                                maxIndex1 = j;
267                                                        }
268                                                }
269                                                alignment.add(tempStorage1.get(maxIndex1,0).getAlignment().get(0));
270                                        }
271                                        for(int i=0; i<seqCount2;i++) {
272                                                for (int j=0;j<seqCount1;j++) {
273                                                        tempStorage2.set(j, 0, new SmithWatermanRepeated(node2.getSequence(i).getSequence(),node1.getSequence(j).getSequence(),submat,5));
274                                                        if(maxScore2 < tempStorage2.get(j, 0).getAlignmentScore()) {
275                                                                maxScore2 = tempStorage2.get(j, 0).getAlignmentScore();
276                                                                maxIndex2 = j;
277                                                        }
278                                                }
279                                                alignment.add(tempStorage2.get(maxIndex2,0).getAlignment().get(0));
280                                        }
281                                       
282                        }
283                        else {
284                                Console.traceln(Level.WARNING,"No sequences to align while merging " + node1.getIdentifier() + " with " + node2.getIdentifier());
285                        }
286                }
287                else {
288                        Console.traceln(Level.WARNING,"More than 2 children! This should never happen, it's a binary tree.");
289                }
290                return alignment;
291        }
292
293       
294        /**
295         * compute updated distance between the new cluster (i,j)
296         * to any other cluster k
297         */
298        private double updatedDistance(int i, int j, int k)
299        {
300                int ai = alias[i];
301                int aj = alias[j];
302               
303                double ocsum = (double) (oc[ai]+oc[aj]);
304                double idist = getDist(k,i);
305                double jdist = getDist(k,j);
306                //TODO: Dirty hack to deal with infinity, insert proper solution here
307                if(Double.isInfinite(idist)) {
308                        idist = 100;
309                }
310                if(Double.isInfinite(jdist)) {
311                        jdist = 100;
312                }
313               
314                return  (oc[ai]/ocsum)*idist+
315                        (oc[aj]/ocsum)*jdist;
316        }
317}
318
319
Note: See TracBrowser for help on using the repository browser.