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

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

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