Context Navigation

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

Last change on this file since 1588 was 1588, checked in by rkrimmel, 10 years ago
Fixed Needleman Wunsch Algorithm
File size: 9.0 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
14	package de.ugoe.cs.autoquest.tasktrees.alignment.pal.tree;
15
16	import java.util.ArrayList;
17	import java.util.logging.Level;
18
19	import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm;
20	import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithmFactory;
21	import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.NumberSequence;
22	import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.SmithWatermanRepeated;
23	import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.BinaryAlignmentStorage;
24	import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.ObjectDistanceSubstitionMatrix;
25	import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.UPGMAMatrix;
26	import de.ugoe.cs.autoquest.tasktrees.alignment.pal.misc.Identifier;
27	import 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	*/
38	public 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	this.getRoot().setSequences(alignSequences(this.getRoot()));
140	distance = null;
141	}
142
143	private void findNextPair()
144	{
145	besti = 0;
146	bestj = 1;
147	double dmin = getDist(0, 1);
148	for (int i = 0; i < numClusters-1; i++)
149	{
150	for (int j = i+1; j < numClusters; j++)
151	{
152	if (getDist(i, j) < dmin)
153	{
154	dmin = getDist(i, j);
155	besti = i;
156	bestj = j;
157	}
158	}
159	}
160	abi = alias[besti];
161	abj = alias[bestj];
162	//System.out.println("Found best pair: " + abi + "/" +abj + " - "+ besti+ "/"+bestj +" with distance " + dmin);
163
164	}
165
166	private void newBranchLengths()
167	{
168	double dij = getDist(besti, bestj);
169
170	getRoot().getChild(besti).setBranchLength(dij/2.0-height[abi]);
171	getRoot().getChild(bestj).setBranchLength(dij/2.0-height[abj]);
172	}
173
174	private void newCluster()
175	{
176	// Update distances
177	for (int k = 0; k < numClusters; k++)
178	{
179	if (k != besti && k != bestj)
180	{
181	int ak = alias[k];
182	double updated = updatedDistance(besti,bestj,k);
183	distance[ak][abi] = distance[abi][ak] = updated;
184	}
185	}
186	distance[abi][abi] = 0.0;
187
188	// Update UPGMA variables
189	height[abi] = getDist(besti, bestj)/2.0;
190	oc[abi] += oc[abj];
191
192	// Index besti now represent the new cluster
193	Node newNode = getRoot().joinChildren(besti, bestj);
194
195	if(newNode instanceof FengDoolittleNode) {
196	newNode.setSequences(alignSequences(newNode));
197	}
198
199	// Update alias
200	for (int i = bestj; i < numClusters-1; i++)
201	{
202	alias[i] = alias[i+1];
203	}
204
205	numClusters--;
206	}
207
208
209	public ArrayList<NumberSequence> alignSequences(Node parent) {
210	ArrayList<NumberSequence> alignment = new ArrayList<NumberSequence>();
211	if(parent.getChildCount()<3) {
212
213	Node node1 = parent.getChild(0);
214	Node node2 = parent.getChild(1);
215
216	int seqCount1 = node1.getSequences().size();
217	int seqCount2 = node2.getSequences().size();
218
219	/*
220	for(int i = 0; i < seqCount1; i++) {
221	for(int j = 0; j < seqCount2; j++) {
222	node1.getSequence(i).printSequence();
223	node2.getSequence(j).printSequence();
224	}
225	}
226	*/
227
228	Console.traceln(Level.INFO,"Merging node " + node1.getIdentifier() + " with " + node2.getIdentifier());
229	//Console.println("SeqCount1: " + seqCount1 + " seqCount2 " + seqCount2);
230	//Align 2 sequences
231	if(seqCount1 == 1 && seqCount2 == 1) {
232	alignment = (alignments.get(node1.getNumber(), node2.getNumber())).getAlignment();
233
234	}
235	//Align a sequence to a group
236	else if( seqCount1 > 1 && seqCount2 == 1) {
237	alignment.addAll(node1.getSequences());
238
239	BinaryAlignmentStorage tempStorage = new BinaryAlignmentStorage(seqCount1,seqCount2);
240	double maxScore = 0.0;
241	int maxIndex = 0;
242	for(int i=0;i<seqCount1;i++) {
243	tempStorage.set(i, 1, AlignmentAlgorithmFactory.create(node1.getSequence(i).getSequence(), node2.getSequence(0).getSequence() , submat, 5));
244	if(maxScore < tempStorage.get(i, 1).getAlignmentScore()) {
245	maxScore = tempStorage.get(i, 1).getAlignmentScore();
246	maxIndex = i;
247	}
248	}
249	//if(maxScore > 0)
250
251	alignment.add(tempStorage.get(maxIndex, 1).getAlignment().get(1));
252	}
253	//Align a sequence to a group
254	else if(seqCount1 == 1 && seqCount2 > 1) {
255	alignment.addAll(node2.getSequences());
256	BinaryAlignmentStorage tempStorage = new BinaryAlignmentStorage(seqCount1,seqCount2);
257	double maxScore = 0.0;
258	int maxIndex = 0;
259	for(int i=0;i<seqCount2;i++) {
260	tempStorage.set(1, i, AlignmentAlgorithmFactory.create(node2.getSequence(i).getSequence(), node1.getSequence(0).getSequence() , submat, 5));
261	if(maxScore < tempStorage.get(1, i).getAlignmentScore()) {
262	maxScore = tempStorage.get(1, i).getAlignmentScore();
263	maxIndex = i;
264	}
265	}
266	//if(maxScore > 0)
267
268	alignment.add(tempStorage.get(1,maxIndex).getAlignment().get(1));
269	}
270	//Align 2 groups
271	else if((seqCount1 > 1) && (seqCount2 > 1)){
272	BinaryAlignmentStorage tempStorage1 = new BinaryAlignmentStorage(seqCount2,1);
273	BinaryAlignmentStorage tempStorage2 = new BinaryAlignmentStorage(seqCount1,1);
274	double maxScore1 = 0.0;
275	double maxScore2 = 0.0;
276	int maxIndex1 = 0;
277	int maxIndex2 = 0;
278	for(int i=0;i<seqCount1;i++) {
279	for(int j=0;j<seqCount2;j++) {
280	tempStorage1.set(j, 0, AlignmentAlgorithmFactory.create(node1.getSequence(i).getSequence(), node2.getSequence(j).getSequence() , submat, 5));
281	if(maxScore1 < tempStorage1.get(j, 0).getAlignmentScore()) {
282	maxScore1 = tempStorage1.get(j, 0).getAlignmentScore();
283	maxIndex1 = j;
284	}
285	}
286	//if(maxScore1 > 0)
287	alignment.add(tempStorage1.get(maxIndex1,0).getAlignment().get(0));
288	}
289	for(int i=0; i<seqCount2;i++) {
290	for (int j=0;j<seqCount1;j++) {
291	tempStorage2.set(j, 0, AlignmentAlgorithmFactory.create(node2.getSequence(i).getSequence(),node1.getSequence(j).getSequence(),submat,5));
292	if(maxScore2 < tempStorage2.get(j, 0).getAlignmentScore()) {
293	maxScore2 = tempStorage2.get(j, 0).getAlignmentScore();
294	maxIndex2 = j;
295	}
296	}
297	//if(maxScore2 > 0)
298	alignment.add(tempStorage2.get(maxIndex2,0).getAlignment().get(0));
299	}
300
301	}
302	else {
303	Console.traceln(Level.WARNING,"No sequences to align while merging " + node1.getIdentifier() + " with " + node2.getIdentifier());
304	}
305	}
306	else {
307	Console.traceln(Level.WARNING,"More than 2 children! This should never happen, it's a binary tree.");
308	}
309	return alignment;
310	}
311
312
313	/**
314	* compute updated distance between the new cluster (i,j)
315	* to any other cluster k
316	*/
317	private double updatedDistance(int i, int j, int k)
318	{
319	int ai = alias[i];
320	int aj = alias[j];
321
322	double ocsum = (double) (oc[ai]+oc[aj]);
323	double idist = getDist(k,i);
324	double jdist = getDist(k,j);
325	//TODO: Dirty hack to deal with infinity, insert proper solution here
326	if(Double.isInfinite(idist)) {
327	idist = 100;
328	}
329	if(Double.isInfinite(jdist)) {
330	jdist = 100;
331	}
332
333	return (oc[ai]/ocsum)*idist+
334	(oc[aj]/ocsum)*jdist;
335	}
336	}
337
338

Note: See TracBrowser for help on using the repository browser.

Download in other formats: