1 // ==========================================================================
2 //                 SeqAn - The Library for Sequence Analysis
3 // ==========================================================================
4 // Copyright (c) 2006-2010, Knut Reinert, FU Berlin
5 // All rights reserved.
6 //
7 // Redistribution and use in source and binary forms, with or without
8 // modification, are permitted provided that the following conditions are met:
9 //
10 //     * Redistributions of source code must retain the above copyright
11 //       notice, this list of conditions and the following disclaimer.
12 //     * Redistributions in binary form must reproduce the above copyright
13 //       notice, this list of conditions and the following disclaimer in the
14 //       documentation and/or other materials provided with the distribution.
15 //     * Neither the name of Knut Reinert or the FU Berlin nor the names of
16 //       its contributors may be used to endorse or promote products derived
17 //       from this software without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
23 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
25 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
26 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
29 // DAMAGE.
30 //
31 // ==========================================================================
32 
33 #ifndef SEQAN_HEADER_GRAPH_ALIGN_TCOFFEE_PROGRESSIVE_H
34 #define SEQAN_HEADER_GRAPH_ALIGN_TCOFFEE_PROGRESSIVE_H
35 
36 namespace SEQAN_NAMESPACE_MAIN
37 {
38 
39 //////////////////////////////////////////////////////////////////////////////
40 // Progressive Alignment
41 //////////////////////////////////////////////////////////////////////////////
42 
43 //////////////////////////////////////////////////////////////////////////////
44 
45 template<typename TStringSet, typename TCargo, typename TSpec, typename TPosition, typename TSequence>
46 inline void
_buildLeafString(Graph<Alignment<TStringSet,TCargo,TSpec>> const & g,TPosition const pos,TSequence & alignSeq)47 _buildLeafString(Graph<Alignment<TStringSet, TCargo, TSpec> > const& g,
48 				 TPosition const pos,
49 				 TSequence& alignSeq)
50 {
51 	SEQAN_CHECKPOINT
52 	typedef Graph<Alignment<TStringSet, TCargo, TSpec> > TGraph;
53 	typedef typename Size<TGraph>::Type TSize;
54 	typedef typename Id<TGraph>::Type TId;
55 	typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
56 	typedef typename Value<TSequence>::Type TVertexString;
57 
58 	//TVertexDescriptor nilVertex = getNil<TVertexDescriptor>();
59 	TStringSet& str = stringSet(g);
60 	TId seqId = positionToId(str, pos);
61 	TSize lenRoot = length(str[pos]);
62 	TSize i = 0;
63 	while(i<lenRoot) {
64 		TVertexDescriptor nextVertex = findVertex(const_cast<TGraph&>(g), seqId, i);
65 		//SEQAN_ASSERT_TRUE(nextVertex != nilVertex);
66 		//if (nextVertex == nilVertex) {
67 		//	std::cout << "Warning: Nil Vertex" << std::endl;
68 		//	TSize j = i + 1;
69 		//	while ((j < lenRoot) && (findVertex(const_cast<TGraph&>(g), seqId, j) == nilVertex)) ++j;
70 		//	nextVertex = addVertex(const_cast<TGraph&>(g), seqId, i, j-i);
71 		//}
72 		appendValue(alignSeq, TVertexString(nextVertex), Generous());
73 		i += fragmentLength(g, nextVertex);
74 	}
75 }
76 
77 //////////////////////////////////////////////////////////////////////////////
78 
79 template<typename TStringSet, typename TCargo, typename TSpec, typename TSegmentString, typename TOutGraph>
80 inline void
_createAlignmentGraph(Graph<Alignment<TStringSet,TCargo,TSpec>> const & g,TSegmentString & alignSeq,TOutGraph & gOut)81 _createAlignmentGraph(Graph<Alignment<TStringSet, TCargo, TSpec> > const& g,
82 					  TSegmentString& alignSeq,
83 					  TOutGraph& gOut)
84 {
85 	SEQAN_CHECKPOINT
86 	typedef Graph<Alignment<TStringSet, TCargo, TSpec> > TGraph;
87 	typedef typename Size<TGraph>::Type TSize;
88 	typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
89 	typedef String<TVertexDescriptor> TVertexString;
90 
91 	// Create the alignment graph
92 	TSize alignSeqLen = length(alignSeq);
93 	for(TSize i = 0; i<alignSeqLen;++i) {
94 		TVertexString& alignSeq_i = alignSeq[i];
95 		TSize len_i = length(alignSeq_i);
96 		for(TSize j=0; j<len_i; ++j) {
97 			TVertexDescriptor v = alignSeq_i[j];
98 			SEQAN_ASSERT_TRUE(fragmentBegin(g,v) < length(getValueById(stringSet(g), sequenceId(g,v))));
99 			SEQAN_ASSERT_TRUE(fragmentLength(g,v) > 0);
100 			SEQAN_ASSERT_TRUE(fragmentBegin(g,v) + fragmentLength(g,v) <= length(getValueById(stringSet(g), sequenceId(g,v))));
101 			TVertexDescriptor l = addVertex(gOut, sequenceId(g, v), fragmentBegin(g,v), fragmentLength(g,v));
102 			//std::cout << l << label(gOut, l) << ',';
103 			TSize count = 1;
104 			for(TSize k = j; k>0; --k) {
105 				SEQAN_ASSERT_TRUE(fragmentLength(gOut,l) == fragmentLength(gOut,l - count));
106 				addEdge(gOut, (TVertexDescriptor) (l - count), (TVertexDescriptor) l);
107 				++count;
108 			}
109 		}
110 		//std::cout << std::endl;
111 	}
112 }
113 
114 //////////////////////////////////////////////////////////////////////////////
115 
116 /**
117 .Function.progressiveAlignment:
118 ..summary:Performs a progressive alignment.
119 ..cat:Graph
120 ..signature:
121 progressiveAlignment(inputGraph, guideTree, outputGraph)
122 ..param.inputGraph:The alignment graph with multiple sequence information.
123 ...type:Spec.Alignment Graph
124 ..param.guideTree:A guide tree.
125 ...type:Spec.Tree
126 ..param.outputGraph:An alignment graph for the final MSA.
127 ...type:Spec.Alignment Graph
128 ..returns:void
129 ..include:seqan/graph_msa.h
130 */
131 template<typename TStringSet, typename TCargo, typename TSpec, typename TGuideTree, typename TOutGraph>
132 inline void
progressiveAlignment(Graph<Alignment<TStringSet,TCargo,TSpec>> & g,TGuideTree & tree,TOutGraph & gOut)133 progressiveAlignment(Graph<Alignment<TStringSet, TCargo, TSpec> >& g,
134 					 TGuideTree& tree,
135 					 TOutGraph& gOut)
136 {
137 	SEQAN_CHECKPOINT
138 	typedef Graph<Alignment<TStringSet, TCargo, TSpec> > TGraph;
139 	typedef typename Size<TGraph>::Type TSize;
140 	typedef typename VertexDescriptor<TGuideTree>::Type TVertexDescriptor;
141 	typedef typename Iterator<TGuideTree, BfsIterator>::Type TBfsIterator;
142 	typedef typename Iterator<TGuideTree, AdjacencyIterator>::Type TAdjacencyIterator;
143 	typedef String<TVertexDescriptor> TVertexString;
144 	typedef String<TVertexString> TSegmentString;
145 
146 	// Initialization
147 	TVertexDescriptor rootVertex = getRoot(tree);
148 	TSize nVertices = numVertices(tree);
149 
150 	// Vertices in reversed bfs order
151 	TVertexString vertices;
152 	resize(vertices, nVertices);
153 
154 	// All Strings of Strings of vertices for each node of the guide tree
155 	String<TSegmentString> segString;
156 	resize(segString, nVertices);
157 
158 	// Walk through the tree in bfs order
159 	typedef typename Iterator<TVertexString, Standard>::Type TVertexIter;
160 	TVertexIter itVert = begin(vertices, Standard());
161 	TVertexIter itVertEnd = end(vertices, Standard());
162 	--itVertEnd;
163 	TBfsIterator bfsIt(tree, rootVertex);
164 	for(;!atEnd(bfsIt);goNext(bfsIt), --itVertEnd)
165 		*itVertEnd = *bfsIt;
166 
167 	// Progressive alignment
168 	itVert = begin(vertices, Standard());
169 	itVertEnd = end(vertices, Standard());
170 	for(;itVert != itVertEnd; ++itVert) {
171 		if(isLeaf(tree, *itVert)) _buildLeafString(g, *itVert, segString[*itVert]);
172 		else {
173 			// Align the two children (Binary tree)
174 			TAdjacencyIterator adjIt(tree, *itVert);
175 			TVertexDescriptor child1 = *adjIt; goNext(adjIt);
176 			heaviestCommonSubsequence(g, segString[child1], segString[*adjIt], segString[*itVert]);
177 			clear(segString[child1]);
178 			clear(segString[*adjIt]);
179 		}
180 	}
181 
182 	// Create the alignment graph
183 	_createAlignmentGraph(g, segString[rootVertex], gOut);
184 }
185 
186 
187 
188 //////////////////////////////////////////////////////////////////////////////
189 // Progressive Matching, just for testing purposes
190 //////////////////////////////////////////////////////////////////////////////
191 
192 
193 //////////////////////////////////////////////////////////////////////////////
194 
195 template<typename TStringSet, typename TCargo, typename TSpec, typename TSegmentString, typename TEdgeMap, typename TOutGraph>
196 inline void
_createMatchingGraph(Graph<Alignment<TStringSet,TCargo,TSpec>> const & g,TSegmentString & alignSeq,TEdgeMap & edgeMap,TOutGraph & gOut,TEdgeMap & edgeMapOut)197 _createMatchingGraph(Graph<Alignment<TStringSet, TCargo, TSpec> > const& g,
198 					 TSegmentString& alignSeq,
199 					 TEdgeMap& edgeMap,
200 					 TOutGraph& gOut,
201 					 TEdgeMap& edgeMapOut)
202 {
203 	typedef Graph<Alignment<TStringSet, TCargo, TSpec> > TGraph;
204 	typedef typename Size<TGraph>::Type TSize;
205 	typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
206 	typedef typename EdgeDescriptor<TGraph>::Type TEdgeDescriptor;
207 	typedef String<TVertexDescriptor> TVertexString;
208 
209 	// Initialization
210 	TVertexDescriptor nilVertex = getNil<TVertexDescriptor>();
211 
212 	// Create the matching graph
213 	clear(edgeMapOut);
214 	TSize alignSeqLen = length(alignSeq);
215 	for(TSize i = 0; i<alignSeqLen;++i) {
216 		TVertexString& alignSeq_i = alignSeq[i];
217 		TSize len_i = length(alignSeq_i);
218 		for(TSize j=0; j<len_i - 1; ++j) {
219 			for(TSize k=j+1; k<len_i; ++k) {
220 				TVertexDescriptor v1 = getValue(alignSeq_i, j);
221 				TVertexDescriptor v2 = getValue(alignSeq_i, k);
222 				if ((v1 == nilVertex) || (v2 == nilVertex)) continue;
223 
224 				TVertexDescriptor v1New = findVertex(gOut, sequenceId(g, v1), fragmentBegin(g,v1));
225 				if (v1New == nilVertex) v1New = addVertex(gOut, sequenceId(g, v1), fragmentBegin(g,v1), fragmentLength(g,v1));
226 				TVertexDescriptor v2New = findVertex(gOut, sequenceId(g, v2), fragmentBegin(g,v2));
227 				if (v2New == nilVertex) v2New = addVertex(gOut, sequenceId(g, v2), fragmentBegin(g,v2), fragmentLength(g,v2));
228 
229 				TEdgeDescriptor e = findEdge(g, v1, v2);
230 				addEdge(gOut, v1New, v2New, cargo(e));
231 				appendValue(edgeMapOut, property(edgeMap, e));
232 			}
233 		}
234 	}
235 }
236 
237 //////////////////////////////////////////////////////////////////////////////
238 
239 template<typename TStringSet, typename TCargo, typename TSpec, typename TString, typename TOutString>
240 inline TCargo
heaviestMatching(Graph<Alignment<TStringSet,TCargo,TSpec>> const & g,TString const & str1,TString const & str2,TOutString & align)241 heaviestMatching(Graph<Alignment<TStringSet, TCargo, TSpec> > const& g,
242 				 TString const& str1,
243 				 TString const& str2,
244 				 TOutString& align)
245 {
246 	typedef Graph<Alignment<TStringSet, TCargo, TSpec> > TGraph;
247 	typedef __int64 TLargeSize;
248 	typedef typename Size<TStringSet>::Type TSize;
249 	typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
250 	typedef typename Iterator<TGraph, OutEdgeIterator>::Type TOutEdgeIterator;
251 	typedef typename Value<TString>::Type TVertexSet;
252 
253 	TSize m = length(str1);  // How many sets of vertex descriptors in seq1
254 	TSize n = length(str2);  // How many sets of vertex descriptors in seq1
255 
256 	// Size of the sequences
257 	// Note for profile alignments every member of the sequence is a String!!! of vertex descriptors
258 	TCargo divider = (TCargo) length(str1[0]) * (TCargo) length(str2[0]);
259 	TVertexDescriptor nilVertex = getNil<TVertexDescriptor>();
260 
261 	// Fill the vertex to position map for str1
262 	// Remember for each vertex descriptor the position in the sequence
263 	typedef std::map<TVertexDescriptor, TSize> TVertexToPosMap;
264 	typedef typename TVertexToPosMap::const_iterator TVertexToPosMapIter;
265 	TVertexToPosMap map;
266 	typedef typename Iterator<TString const>::Type TStringIterConst;
267 	typedef typename Iterator<TVertexSet const>::Type TVertexSetIterConst;
268 	TStringIterConst itStrEnd1 = end(str1);
269 	TSize pos = 0;
270 	for(TStringIterConst itStr1 = begin(str1);itStr1 != itStrEnd1;++itStr1, ++pos) {
271 		TVertexSetIterConst itVEnd = end(getValue(itStr1));
272 		for(TVertexSetIterConst itV = begin(getValue(itStr1));itV != itVEnd;++itV) {
273 			if (*itV != nilVertex) map.insert(std::make_pair(*itV, pos));
274 		}
275 	}
276 
277 	// Given a full graph, what positions are occupied?
278 	typedef std::set<TLargeSize> TOccupiedPositions;
279 	TOccupiedPositions occupiedPositions;
280 	TStringIterConst itStrEnd2 = end(str2);
281 	TSize posItStr2 = 0;
282 	for(TStringIterConst itStr2 = begin(str2);itStr2 != itStrEnd2;++itStr2, ++posItStr2) {
283 		TVertexSetIterConst itVEnd = end(getValue(itStr2));
284 		for(TVertexSetIterConst itV = begin(getValue(itStr2));itV != itVEnd;++itV) {
285 			if (*itV != nilVertex) {
286 				TOutEdgeIterator itOut(g, *itV);
287 				for(;!atEnd(itOut); ++itOut) {
288 					// Target vertex must be in the map
289 					TVertexToPosMapIter pPos = map.find(targetVertex(itOut));
290 					if (pPos != map.end()) occupiedPositions.insert( (TLargeSize) (pPos->second) * (TLargeSize) n + (TLargeSize) (n - posItStr2 - 1) );
291 				}
292 			}
293 		}
294 	}
295 	// Map the occupied positions to slots
296 	typedef std::map<TLargeSize, TSize> TPositionToSlotMap;
297 	TPositionToSlotMap posToSlotMap;
298 	TSize counter = 0;
299 	for(typename TOccupiedPositions::const_iterator setIt = occupiedPositions.begin();setIt != occupiedPositions.end(); ++setIt, ++counter) {
300 		posToSlotMap.insert(std::make_pair(*setIt, counter));
301 	}
302 	occupiedPositions.clear();
303 
304 	// Walk through str2 and fill in the weights of the actual edges
305 	typedef String<TCargo> TWeights;
306 	typedef typename Iterator<TWeights>::Type TWeightsIter;
307 	TWeights weights;
308 	resize(weights, posToSlotMap.size(),0);
309 	posItStr2 = 0;
310 	for(TStringIterConst itStr2 = begin(str2);itStr2 != itStrEnd2;++itStr2, ++posItStr2) {
311 		TVertexSetIterConst itVEnd = end(getValue(itStr2));
312 		for(TVertexSetIterConst itV = begin(getValue(itStr2));itV != itVEnd;++itV) {
313 			if (*itV != nilVertex) {
314 				TOutEdgeIterator itOut(g, *itV);
315 				for(;!atEnd(itOut); ++itOut) {
316 					// Target vertex must be in the map
317 					TVertexToPosMapIter pPos = map.find(targetVertex(itOut));
318 					if (pPos != map.end()) weights[posToSlotMap[ (TLargeSize) (pPos->second) * (TLargeSize) n + (TLargeSize) (n - posItStr2 - 1) ]] += (TCargo) cargo(*itOut);
319 				}
320 			}
321 		}
322 	}
323 	map.clear();
324 	// Average weights
325 	TWeightsIter itWeights = begin(weights);
326 	TWeightsIter itWeightsEnd = begin(weights);
327 	for(;itWeights != itWeightsEnd; ++itWeights) *itWeights /= divider;
328 
329 	// Create the match graph
330 	Graph<Undirected<void> > matchGraph;
331 	typedef typename VertexDescriptor<Graph<Undirected<void> > >::Type TVD;
332 	String<bool> vertexMap;
333 	resize(vertexMap, m + n, false);
334 	for(TSize i = 0; i < m; ++i) addVertex(matchGraph);
335 	for(TSize i = 0; i < n; ++i) value(vertexMap, addVertex(matchGraph)) = true;
336 	for(typename TPositionToSlotMap::const_iterator mapIt = posToSlotMap.begin();mapIt != posToSlotMap.end(); ++mapIt) {
337 		TLargeSize pos = mapIt->first;
338 		TSize i = (TSize) (pos / (TLargeSize) n);   // Get the index in str1
339 		TSize j = n - 1 - (TSize) (pos % (TLargeSize) n); // Get the index in str2
340 		addEdge(matchGraph, i, m+j);
341 	}
342 
343 	// Compute the best matching
344 	typedef String<Pair<TVD, TVD> > TEdges;
345 	TEdges edges;
346 
347 	//std::fstream strm;
348 	//strm.open("Z:\\my_graph.dot", std::ios_base::out | std::ios_base::trunc);
349 	//write(strm,matchGraph,DotDrawing());
350 	//strm.close();
351 
352 	TCargo val = weightedBipartiteMatching(matchGraph, vertexMap, weights, edges);
353 
354 	// Retrieve the aligned segments
355 	TSize seqsInStr1 = length(str1[0]);
356 	TSize seqsInStr2 = length(str2[0]);
357 	typedef typename Value<TString>::Type TVertexSet;
358 	typedef typename Iterator<TString const, Rooted>::Type TStringIter;
359 	typedef typename Iterator<TString, Rooted>::Type TSIter;
360 	typedef typename Iterator<TVertexSet const, Rooted>::Type TVertexSetIter;
361 	typedef typename Iterator<TVertexSet, Rooted>::Type TIter;
362 	clear(align);
363 	// Retrieve all matches
364 	String<bool> matchedVertices;
365 	resize(matchedVertices, n + m, false);
366 	typedef typename Iterator<TEdges>::Type TEdgesIter;
367 	TEdgesIter itEdges = begin(edges);
368 	TEdgesIter itEdgesEnd = end(edges);
369 	for(;itEdges != itEdgesEnd; ++itEdges) {
370 		TSize i = (value(itEdges)).i1;
371 		TSize j = (value(itEdges)).i2 - m;
372 		value(matchedVertices, i) = true;
373 		value(matchedVertices, m+j) = true;
374 		TVertexSet tmp;
375 		resize(tmp, (seqsInStr1 + seqsInStr2), nilVertex);
376 		TVertexSetIter itVEnd = end(value(str1,i));
377 		TSize count = 0;
378 		for(TVertexSetIter itV = begin(value(str1,i));itV != itVEnd;++itV) {
379 			tmp[count] = *itV;
380 			++count;
381 		}
382 		TVertexSetIter itVEnd2 = end(value(str2,j));
383 		for(TVertexSetIter itV2 = begin(value(str2,j));itV2 != itVEnd2;++itV2) {
384 			tmp[count] = *itV2;
385 			++count;
386 		}
387 		appendValue(align, tmp);
388 	}
389 	typedef typename Iterator<Graph<Undirected<void> >, VertexIterator>::Type TVertexIterator;
390 	TVertexIterator itVertex(matchGraph);
391 	for(;!atEnd(itVertex);++itVertex) {
392 		if (getProperty(matchedVertices, *itVertex) == true) continue;
393 		if (*itVertex < m) {
394 			TSize i = *itVertex;
395 			TVertexSet tmp;
396 			resize(tmp, (seqsInStr1 + seqsInStr2), nilVertex);
397 			TVertexSetIter itVEnd = end(value(str1, i));
398 			TSize count = 0;
399 			for(TVertexSetIter itV = begin(value(str1, i));itV != itVEnd;++itV) {
400 				tmp[count] = *itV;
401 				++count;
402 			}
403 			appendValue(align, tmp);
404 		} else {
405 			TSize j = *itVertex - m;
406 			TVertexSet tmp;
407 			resize(tmp, (seqsInStr1 + seqsInStr2), nilVertex);
408 			TVertexSetIter itVEnd = end(value(str2, j));
409 			TSize count = 0;
410 			for(TVertexSetIter itV = begin(value(str2, j));itV != itVEnd;++itV) {
411 				tmp[count] = *itV;
412 				++count;
413 			}
414 			appendValue(align, tmp);
415 		}
416 	}
417 
418 	return val;
419 }
420 
421 //////////////////////////////////////////////////////////////////////////////
422 
423 template<typename TStringSet, typename TCargo, typename TSpec, typename TGuideTree, typename TVertexDescriptor, typename TSequence>
424 inline void
_recursiveProgressiveMatching(Graph<Alignment<TStringSet,TCargo,TSpec>> & g,TGuideTree & tree,TVertexDescriptor const root,TSequence & alignSeq)425 _recursiveProgressiveMatching(Graph<Alignment<TStringSet, TCargo, TSpec> >& g,
426 							  TGuideTree& tree,
427 							  TVertexDescriptor const root,
428 							  TSequence& alignSeq)
429 {
430 	typedef Graph<Alignment<TStringSet, TCargo, TSpec> > TGraph;
431 	typedef typename Size<TGraph>::Type TSize;
432 	typedef typename Id<TGraph>::Type TId;
433 	typedef typename Iterator<TGuideTree, AdjacencyIterator>::Type TAdjacencyIterator;
434 	typedef typename Iterator<TGuideTree, DfsPreorder>::Type TDfsPreorderIterator;
435 
436 	if(isLeaf(tree, root)) {
437 		_buildLeafString(g, root, alignSeq);
438 	} else {
439 		// Align the two children (Binary tree)
440 		typedef String<String<TVertexDescriptor> > TSegmentString;
441 		TSegmentString seq1;
442 		TSegmentString seq2;
443 		TAdjacencyIterator adjIt(tree, root);
444 		_recursiveProgressiveMatching(g,tree, *adjIt, seq1);
445 		goNext(adjIt);
446 		_recursiveProgressiveMatching(g,tree, *adjIt, seq2);
447 
448 		heaviestMatching(g, seq1, seq2, alignSeq);
449 
450 		//// Debug Code
451 		//for(TSize i = 0; i<length(alignSeq);++i) {
452 		//	std::cout << '(';
453 		//	for(TSize j=0; j<length(alignSeq[i]);++j) {
454 		//		std::cout << getValue(alignSeq[i], j) << ',';
455 		//	}
456 		//	std::cout << ')' << ',';
457 		//}
458 		//std::cout << std::endl;
459 	}
460 }
461 
462 //////////////////////////////////////////////////////////////////////////////
463 
464 template<typename TStringSet, typename TCargo, typename TSpec, typename TGuideTree, typename TEdgeMap, typename TOutGraph>
465 inline void
progressiveMatching(Graph<Alignment<TStringSet,TCargo,TSpec>> & g,TGuideTree & tree,TEdgeMap & edgeMap,TOutGraph & gOut,TEdgeMap & edgeMapOut)466 progressiveMatching(Graph<Alignment<TStringSet, TCargo, TSpec> >& g,
467 					TGuideTree& tree,
468 					TEdgeMap& edgeMap,
469 					TOutGraph& gOut,
470 					TEdgeMap& edgeMapOut)
471 {
472 	typedef Graph<Alignment<TStringSet, TCargo, TSpec> > TGraph;
473 	typedef typename Size<TGraph>::Type TSize;
474 	typedef typename VertexDescriptor<TGraph>::Type TVertexDescriptor;
475 	typedef String<TVertexDescriptor> TVertexString;
476 	typedef String<TVertexString> TSegmentString;
477 
478 	// Perform progressive alignment
479 	TSegmentString alignSeq;
480 	_recursiveProgressiveMatching(g,tree,getRoot(tree),alignSeq);
481 
482 	// Create the alignment graph
483 	_createMatchingGraph(g, alignSeq, edgeMap, gOut, edgeMapOut);
484 }
485 
486 
487 }// namespace SEQAN_NAMESPACE_MAIN
488 
489 #endif //#ifndef SEQAN_HEADER_...
490