seqan-dev 1.3.1-1 / usr / include / seqan / consensus / consensus_library.h

// ==========================================================================
//                 SeqAn - The Library for Sequence Analysis
// ==========================================================================
// Copyright (c) 2006-2010, Knut Reinert, FU Berlin
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * Redistributions in binary form must reproduce the above copyright
//       notice, this list of conditions and the following disclaimer in the
//       documentation and/or other materials provided with the distribution.
//     * Neither the name of Knut Reinert or the FU Berlin nor the names of
//       its contributors may be used to endorse or promote products derived
//       from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
// DAMAGE.
//
// ==========================================================================
// Author: Tobias Rausch <rausch@embl.de>
// ==========================================================================

#ifndef SEQAN_HEADER_CONSENSUS_LIBRARY_H
#define SEQAN_HEADER_CONSENSUS_LIBRARY_H

namespace SEQAN_NAMESPACE_MAIN
{

//////////////////////////////////////////////////////////////////////////////

template<typename TSize, typename TCargo1, typename TCargo2>
inline void 
_getAlignmentStatistics(Nothing&,
						 TSize,
						 TSize,
						 TSize,
						 TCargo1,
						 TCargo2)
{
	SEQAN_CHECKPOINT
}

//////////////////////////////////////////////////////////////////////////////

template<typename TValue, typename TSpec, typename TSize, typename TCargo1, typename TCargo2>
inline void 
_getAlignmentStatistics(String<TValue, TSpec>& dist,
						 TSize i,
						 TSize j,
						 TSize nseq,
						 TCargo1,
						 TCargo2 quality)
{
	SEQAN_CHECKPOINT
	dist[i*nseq + j] = (TValue) (100 - quality);
}

//////////////////////////////////////////////////////////////////////////////

template<typename TCargo, typename TSpec, typename TSize, typename TCargo1, typename TCargo2>
inline void 
_getAlignmentStatistics(Graph<Undirected<TCargo, TSpec> >& dist,
						 TSize i,
						 TSize j,
						 TSize,
						 TCargo1,
						 TCargo2 quality)
{
	SEQAN_CHECKPOINT
	addEdge(dist, i, j, (TCargo) (100 - quality));
}

//////////////////////////////////////////////////////////////////////////////
// Layout-based pair selection
//////////////////////////////////////////////////////////////////////////////

template<typename TSize>
struct LessPair_ :
	public ::std::unary_function<Pair<TSize, TSize>, bool>
{
	inline bool 
	operator() (Pair<TSize, TSize> const& a1, Pair<TSize, TSize> const& a2) const {
		if (a1.i1 == a2.i1) return (a1.i2 < a2.i2);
		else return (a1.i1 < a2.i1);
	}
};

template<typename TSize>
struct _LessTripel :
	public ::std::unary_function<Pair<TSize, Triple<TSize, TSize, TSize> >, bool>
{
	inline bool 
	operator() (Pair<TSize, Triple<TSize, TSize, TSize> > const& a1, Pair<TSize, Triple<TSize, TSize, TSize> > const& a2) {
		return (a1.i1 < a2.i1);
	}
};


//////////////////////////////////////////////////////////////////////////////

template<typename TString, typename TSpec, typename TBegEndPos, typename TSize, typename TPairList, typename TPos, typename TSpec2>
inline void 
selectPairs(StringSet<TString, TSpec> const& str,
			TBegEndPos const& begEndPos,
			TSize bandwidth,							
			TPairList& pList,
			String<Pair<TPos, TPos>, TSpec2>& dList)
{
	SEQAN_CHECKPOINT
	typedef String<Pair<TPos, TPos>, TSpec2>  TDistanceList;
	typedef StringSet<TString, TSpec> TStringSet;
	typedef Pair<TPos, TPos> TDiagPair;
	typedef typename Value<TPairList>::Type TPair;
	typedef typename Iterator<TPairList, Standard>::Type TPairIter;
	typedef typename Iterator<TBegEndPos, Standard>::Type TBegEndIter;

	// Initialization
	TSize nseq = length(str);

	// Workaround for strange celera behaviour (just for contained reads)
#ifdef CELERA_OFFSET
	TSize contained_offset=200;
#else
	TSize contained_offset=0;
#endif
	
	// Sort the reads by their first index position
	TBegEndIter begEndIt = begin(begEndPos, Standard());
	TBegEndIter begEndItEnd = end(begEndPos, Standard());
	typedef Triple<TSize, TSize, TSize> TInfo;
	typedef String<Pair<TSize, TInfo> > TPosIndexList;
	typedef typename Iterator<TPosIndexList, Standard>::Type TPosIter;
	TPosIndexList posIndex;
	resize(posIndex, length(begEndPos));
	TPosIter posIndexIt = begin(posIndex, Standard());
	TPosIter posIndexItEnd = end(posIndex, Standard());
	for(TSize index = 0;begEndIt != begEndItEnd; ++begEndIt, ++posIndexIt, ++index) 
		*posIndexIt = ((*begEndIt).i1 < (*begEndIt).i2) ? Pair<TSize, TInfo>((*begEndIt).i1,TInfo(index, (*begEndIt).i1, (*begEndIt).i2)) : Pair<TSize, TInfo>((*begEndIt).i2,TInfo(index, (*begEndIt).i1, (*begEndIt).i2));
	std::sort(begin(posIndex, Standard() ), end(posIndex, Standard() ), _LessTripel<TSize>() );

	// The expected overlap by a pair of reads (represented by its index)
	typedef String<Pair<TSize, TSize> > TOverlapIndexList;
	TOverlapIndexList ovlIndex;
	TSize pairLen = 0;  // Pair Counter
	TPos const initialRadius = (bandwidth + 1) / 2;
	TPos const lengthDivider = 5;	// Overlap / 2^lengthDivider is added to the radius

	// Find all overlapping reads
	TDistanceList preDList;
	TPairList prePList;
	reserve(preDList, nseq * 40);
	reserve(prePList, nseq * 40);
	reserve(ovlIndex, nseq * 40);
	posIndexIt = begin(posIndex, Standard());
	posIndexItEnd = end(posIndex, Standard());
	for(;posIndexIt != posIndexItEnd; ++posIndexIt) {
		TSize index1 = ((*posIndexIt).i2).i1;
		TPos posIi1 = ((*posIndexIt).i2).i2;
		TPos posIi2 = ((*posIndexIt).i2).i3;
		bool forwardI = (posIi1 < posIi2) ? true : false;
		TSize lenI = (posIi1 < posIi2) ? posIi2 - posIi1 : posIi1 - posIi2;
		TPosIter posIndexIt2 = posIndexIt;
		++posIndexIt2;
		for(;posIndexIt2 != posIndexItEnd; ++posIndexIt2) {
			if ((*posIndexIt).i1 + lenI <= (*posIndexIt2).i1) break;
			TSize index2 = ((*posIndexIt2).i2).i1;
			TPos posJi1 = ((*posIndexIt2).i2).i2;
			TPos posJi2 = ((*posIndexIt2).i2).i3;

			// Diagonal boundaries of the band
			// Initialization values are used if one read is contained in the other 
			TSize lenJ = (posJi1 < posJi2) ? posJi2 - posJi1 : posJi1 - posJi2;
			bool forwardJ = (posJi1 < posJi2) ? true : false;
			TPos diagLow = -1 * (TPos) lenJ;
			TPos diagHigh = (TPos) lenI;
			TPos radius = initialRadius;		// Increased by overlap length

			// Read orientations
			if (forwardI) {
				// 1) Forward - Forward
				if (forwardJ) {
					if ((posJi2 < posIi2) && (posJi1 < posIi1)) {
						TPos offset = (posIi1 - posJi1);
						radius += (posJi2 - posJi1 - offset) >> lengthDivider;
						if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
						if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
					} else if ((posJi1 > posIi1) && (posJi2 > posIi2)) {
						TPos offset = (posJi1 - posIi1);
						radius += (posIi2 - posIi1 - offset) >> lengthDivider;
						if (offset + radius < diagHigh) diagHigh = offset + radius;
						if (offset - radius > diagLow) diagLow = offset - radius;
					} else {
						radius += contained_offset;
						if (posIi1 < posJi1) {
							TPos offset = (posJi1 - posIi1);
							radius += (posJi2 - posJi1) >> lengthDivider;
							if (offset + radius < diagHigh) diagHigh = offset + radius;
							if (offset - radius > diagLow) diagLow = offset - radius;
						} else {
							TPos offset = (posIi1 - posJi1);
							radius += (posIi2 - posIi1) >> lengthDivider;
							if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
							if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
						}
					}
				} else { // 2) Forward - Reverse
					if ((posJi1 < posIi2) && (posJi2 < posIi1)) {
						TPos offset = (posIi1 - posJi2);
						radius += (posJi1 - posJi2 - offset) >> lengthDivider;
						if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
						if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
					} else if ((posJi2 > posIi1) && (posJi1 > posIi2)) {
						TPos offset = (posJi2 - posIi1);
						radius += (posIi2 - posIi1 - offset) >> lengthDivider;
						if (offset + radius < diagHigh) diagHigh = offset + radius;
						if (offset - radius > diagLow) diagLow = offset - radius;
					} else {
						radius += contained_offset;  
						if (posIi1 < posJi2) {
							TPos offset = (posJi2 - posIi1);
							radius += (posJi1 - posJi2) >> lengthDivider;
							if (offset + radius < diagHigh) diagHigh = offset + radius;
							if (offset - radius > diagLow) diagLow = offset - radius;
						} else {
							TPos offset = (posIi1 - posJi2);
							radius += (posIi2 - posIi1) >> lengthDivider;
							if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
							if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
						}	
					}
				}
			} else { 
				// 3) Reverse - Forward
				if (forwardJ) {
					if ((posIi1 > posJi2) && (posIi2 > posJi1)) {
						TPos offset = (posIi2 - posJi1);
						radius += (posJi2 - posJi1 - offset) >> lengthDivider;
						if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
						if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
					} else if ((posJi1 > posIi2) && (posJi2 > posIi1)) {
						TPos offset = (posJi1 - posIi2);
						radius += (posIi1 - posIi2 - offset) >> lengthDivider;
						if (offset + radius < diagHigh) diagHigh = offset + radius;
						if (offset - radius > diagLow) diagLow = offset - radius;
					} else {
						radius += contained_offset;  
						if (posIi2 < posJi1) {
							TPos offset = (posJi1 - posIi2);
							radius += (posJi2 - posJi1) >> lengthDivider;
							if (offset + radius < diagHigh) diagHigh = offset + radius;
							if (offset - radius > diagLow) diagLow = offset - radius;
						} else {
							TPos offset = (posIi2 - posJi1);
							radius += (posIi1 - posIi2) >> lengthDivider;
							if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
							if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
						}
					}
				} else { // 4) Reverse - Reverse
					if ((posJi1 < posIi1) && (posJi2 < posIi2)) {
						TPos offset = (posIi2 - posJi2);
						radius += (posJi1 - posJi2 - offset) >> lengthDivider;
						if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
						if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
					} else if ((posJi2 > posIi2) && (posJi1 > posIi1)) {
						TPos offset = (posJi2 - posIi2);
						radius += (posIi1 - posIi2 - offset) >> lengthDivider;
						if (offset + radius < diagHigh) diagHigh = offset + radius;
						if (offset - radius > diagLow) diagLow = offset - radius;
					} else {
						radius += contained_offset;  
						if (posIi2 < posJi2) {
							TPos offset = (posJi2 - posIi2);
							radius += (posJi1 - posJi2) >> lengthDivider;
							if (offset + radius < diagHigh) diagHigh = offset + radius;
							if (offset - radius > diagLow) diagLow = offset - radius;
						} else {
							TPos offset = (posIi2 - posJi2);
							radius += (posIi1 - posIi2) >> lengthDivider;
							if (-1 * offset + radius < diagHigh) diagHigh = -1 * offset + radius;
							if (-1 * offset - radius > diagLow) diagLow = -1 * offset - radius;
						}
					}
				}
			}

			// Append this pair of reads
			if (index1 < index2) {
				appendValue(prePList, TPair(positionToId(str, index1), positionToId(str, index2)), Generous());
				appendValue(preDList, TDiagPair(diagLow, diagHigh), Generous());
			} else {
				appendValue(prePList, TPair(positionToId(str, index2), positionToId(str, index1)), Generous());
				appendValue(preDList, TDiagPair(-1 * diagHigh, -1 * diagLow), Generous());
			}
			// Estimate the overlap quality
			TPos avgDiag = (diagLow + diagHigh) / 2;
			if (avgDiag < 0) avgDiag *= -1;
			appendValue(ovlIndex, Pair<TSize, TSize>((TSize) (avgDiag), pairLen), Generous()); 
			++pairLen;
		}
	}

	// Sort the pairs, better expected overlaps come first
	std::sort(begin(ovlIndex, Standard() ), end(ovlIndex, Standard() ), LessPair_<TSize>() );
	typedef typename Iterator<TOverlapIndexList, Standard>::Type TOVLIter;
	TOVLIter itOvl = begin(ovlIndex, Standard());
	TOVLIter itOvlEnd = end(ovlIndex, Standard());
	reserve(dList, pairLen);
	reserve(pList, pairLen);
	for(;itOvl != itOvlEnd; ++itOvl) {
		TSize count = (*itOvl).i2;
		appendValue(dList, preDList[count], Generous());
		appendValue(pList, prePList[count], Generous());
	}
}



//////////////////////////////////////////////////////////////////////////////

template<typename TString, typename TSpec, typename TBegEndPos, typename TSize, typename TPairList, typename TPos, typename TSpec2>
inline void 
selectPairsIndel(StringSet<TString, TSpec> const& str,
				 TBegEndPos const& begEndPos,
				 TSize lookAround,
				 TPairList& pList,
				 String<Pair<TPos, TPos>, TSpec2>& dList)
{	
	SEQAN_CHECKPOINT	
	typedef String<Pair<TPos, TPos>, TSpec2>  TDistanceList;
	typedef StringSet<TString, TSpec> TStringSet;
	typedef Pair<TPos, TPos> TDiagPair;
	typedef typename Value<TPairList>::Type TPair;
	typedef typename Iterator<TBegEndPos, Standard>::Type TBegEndIter;
	
	TBegEndIter beIt = begin(begEndPos, Standard());
	TBegEndIter beItEnd = end(begEndPos, Standard());
	TSize index1 = 0;
	for(;beIt != beItEnd; ++beIt, ++index1) {
		TPos beg2 = (beIt->i2 < beIt->i1) ? beIt->i2 : beIt->i1;
		TPos diagHigh = (beIt->i2 < beIt->i1) ? beIt->i1 - beIt->i2 : beIt->i2 - beIt->i1;
		TBegEndIter beIt2 = beIt;
		TSize index2 = index1 + 1;
		for(;++beIt2 != beItEnd; ++index2) {
			TPos beg1 = (beIt2->i2 < beIt2->i1) ? beIt2->i2 : beIt2->i1;
			TPos diagLow = (beIt2->i2 < beIt2->i1) ? beIt2->i2 - beIt2->i1 : beIt2->i1 - beIt2->i2;
			TPos diff = (beg1 > beg2) ? beg1 - beg2 : beg2 - beg1;
			if (diff < (TPos) lookAround) {
				appendValue(pList, TPair(positionToId(str, index1), positionToId(str, index2)), Generous());
				appendValue(dList, TDiagPair(diagLow, diagHigh), Generous());
			}
		}
	}
}


//////////////////////////////////////////////////////////////////////////////

template<typename TString, typename TSpec, typename TId, typename TDiagList, typename TBegEndPos, typename TScore, typename TSize, typename TSegmentMatches, typename TScoreValues, typename TDistance>
inline void 
appendSegmentMatches(StringSet<TString, TSpec> const& str,
					 String<Pair<TId, TId> > const& pList,
					 TDiagList const& dList,
					 TBegEndPos const& begEndPos,
					 TScore const& score_type,
					 TSize thresholdMatchlength,
					 TSize thresholdQuality,
					 TSize maxOvl,
					 TSegmentMatches& matches,
					 TScoreValues& scores,
					 TDistance& dist,
					 OverlapLibrary)
{
	SEQAN_CHECKPOINT
	typedef StringSet<TString, Dependent<> > TStringSet;
	typedef String<Pair<TId, TId> > TPairList;
	typedef typename Value<TScoreValues>::Type TScoreValue;
	typedef typename Iterator<TPairList, Standard>::Type TPairIter;
	typedef typename Iterator<TDiagList, Standard>::Type TDiagIter;

	// Initialization
	TSize nseq = length(str);
	_resizeWithRespectToDistance(dist, nseq);

	// "Front" and "Back"-overlap counter for each read
	String<TSize> frontOvl;
	String<TSize> backOvl;
	resize(frontOvl, nseq, 0);
	resize(backOvl, nseq, 0);
	
	// Pairwise alignments
	String<bool> aligned;
	resize(aligned, length(pList), true);
	typedef Iterator<String<bool>, Standard>::Type TBoolIter;
	TBoolIter itAligned = begin(aligned, Standard());
	TPairIter itPair = begin(pList, Standard());
	TDiagIter itDiag = begin(dList, Standard());
	TPairIter itPairEnd = end(pList, Standard());
	TSize dropCount = 0;
	for(;itPair != itPairEnd; ++itPair, ++itDiag, ++itAligned) {
		TId id1 = itPair->i1;
		TId id2 = itPair->i2;
		TSize seq1 = idToPosition(str, id1);
		TSize seq2 = idToPosition(str, id2);
		if ((frontOvl[seq1] > maxOvl) && (backOvl[seq1] > maxOvl) &&
			(frontOvl[seq2] > maxOvl) && (backOvl[seq2] > maxOvl)) {
				++dropCount;
				continue;
		}


		// Make a pairwise string-set
		TStringSet pairSet;
		assignValueById(pairSet, const_cast<StringSet<TString, TSpec>&>(str), id1);
		assignValueById(pairSet, const_cast<StringSet<TString, TSpec>&>(str), id2);
		

		// Overlap alignment
		TSize from = length(matches);
		TScoreValue myScore = globalAlignment(matches, pairSet, score_type, AlignConfig<true,true,true,true>(), itDiag->i1, itDiag->i2, BandedGotoh() );
		TSize to = length(matches);

		// Determine a sequence weight
		TSize matchLen = 0;
		TSize overlapLen = 0;
		TSize alignLen = 0;
		if (from != to) getAlignmentStatistics(matches, pairSet, from, to, matchLen, overlapLen, alignLen);


		// Get only the good overlap alignments
		if ((overlapLen) && ((matchLen * 100) / overlapLen >= thresholdQuality) && (matchLen >= thresholdMatchlength)) {

			//// Debug Code
			//Graph<Alignment<TStringSet, TSize> > tmp(pairSet);
			//globalAlignment(tmp, pairSet, score_type, AlignConfig<true,true,true,true>(), (value(itDiag)).i1, (value(itDiag)).i2, BandedGotoh() );
			////globalAlignment(tmp, pairSet, score_type, Gotoh() );
			//std::cout << "Match length: " << matchLen << std::endl;
			//std::cout << "Overlap length: " << overlapLen << std::endl;
			//std::cout << "Align length: " << alignLen << std::endl;
			//std::cout << "Quality: " << quality << std::endl;
			//std::cout << tmp << std::endl;

			// Create a corresponding edge
			if (seq1<seq2) _getAlignmentStatistics(dist, seq1, seq2, nseq, matchLen, (matchLen * 100) / overlapLen);
			else _getAlignmentStatistics(dist, seq2, seq1, nseq, matchLen, (matchLen * 100) / overlapLen);
			
			// Record the scores
			resize(scores, to);
			typedef typename Iterator<TScoreValues, Standard>::Type TScoreIter;
			TScoreIter itScore = begin(scores, Standard());
			TScoreIter itScoreEnd = end(scores, Standard());
			itScore += from;
			for(;itScore != itScoreEnd; ++itScore) 
				*itScore = myScore;

			// Update the overlap counter
			TSize lenLast = matches[from].len; 
			if (matches[to - 1].begin1 == 0) ++frontOvl[seq1];
			if (matches[to - 1].begin2 == 0) ++frontOvl[seq2];
			if (matches[from].begin1 + lenLast == length(pairSet[0])) ++backOvl[seq1];
			if (matches[from].begin2 + lenLast == length(pairSet[1])) ++backOvl[seq2];
		} else {
			resize(matches, from);
			*itAligned = false;
		}
	}
	//std::cout << "Filtration ration: " << (double) dropCount / (double) length(pList) << std::endl;

	// Find sequences that have no overlap in the front or back
	String<TSize> noFront;
	String<TSize> noBack;
	for(TSize seqI = 0; seqI < nseq; ++seqI) {
		if (frontOvl[seqI] == 0) appendValue(noFront, seqI, Generous());
		else if (backOvl[seqI] == 0) appendValue(noBack, seqI, Generous());
	}
	// Drop the first and the last sequence
	typedef typename Iterator<TBegEndPos, Standard>::Type TBegEndIter;
	TBegEndIter begEndIt = begin(begEndPos, Standard());
	TBegEndIter begEndItEnd = end(begEndPos, Standard());
	TSize minVal = maxValue<TSize>();
	TSize maxVal = 0;
	for(;begEndIt != begEndItEnd; ++begEndIt) {
		TSize pos1 = begEndIt->i1;
		TSize pos2 = begEndIt->i2;
		if (pos1 > pos2) { TSize tmp = pos1; pos1 = pos2; pos2 = tmp;}
		if (pos1 < minVal) minVal = pos1;
		if (pos2 > maxVal) maxVal = pos2;
	}
	// Insert all remaining sequences into a string
	String<TSize> unalignedReads;
	for(TSize i = 0; i < (TSize) length(noFront); ++i) {
		TSize p1 = begEndPos[noFront[i]].i1;
		TSize p2 = begEndPos[noFront[i]].i2;
		if (p1 > p2) {TSize tmp = p1; p1 = p2; p2 = tmp; }
		if (p1 != minVal) appendValue(unalignedReads, noFront[i], Generous());
	}
	for(TSize i = 0; i < (TSize) length(noBack); ++i) {
		TSize p1 = begEndPos[noBack[i]].i1;
		TSize p2 = begEndPos[noBack[i]].i2;
		if (p1 > p2) {TSize tmp = p1; p1 = p2; p2 = tmp; }
		if (p2 != maxVal) appendValue(unalignedReads, noBack[i], Generous());
	}
	TSize countUnalignedReads = length(unalignedReads);
	//std::cout << "Unaligned reads: " << countUnalignedReads << std::endl;
	if (countUnalignedReads > 0) {
		// Sort unaligned reads
		::std::sort(begin(unalignedReads, Standard()), end(unalignedReads, Standard()));
		
		// Realign all unaligned sequences
		itPair = begin(pList, Standard());
		itDiag = begin(dList, Standard());
		itAligned = begin(aligned, Standard());
		for(;itPair != itPairEnd; ++itPair, ++itDiag, ++itAligned) {
			if (*itAligned == true) continue;
			TId id1 = itPair->i1;
			TId id2 = itPair->i2;
			TSize seq1 = idToPosition(str, id1);
			TSize seq2 = idToPosition(str, id2);
			if ((!::std::binary_search(begin(unalignedReads, Standard()), end(unalignedReads, Standard()), seq1)) &&
				(!::std::binary_search(begin(unalignedReads, Standard()), end(unalignedReads, Standard()), seq2))) 
				continue;
			if ((frontOvl[seq1] > maxOvl) && (backOvl[seq1] > maxOvl) &&
				(frontOvl[seq2] > maxOvl) && (backOvl[seq2] > maxOvl)) 
				continue;

			// Make a pairwise string-set
			TStringSet pairSet;
			assignValueById(pairSet, const_cast<StringSet<TString, TSpec>&>(str), id1);
			assignValueById(pairSet, const_cast<StringSet<TString, TSpec>&>(str), id2);
		
			// Overlap alignment
			TSize from = length(matches);
			TScoreValue myScore = globalAlignment(matches, pairSet, score_type, AlignConfig<true,true,true,true>(), Gotoh() );
			TSize to = length(matches);

			// Determine a sequence weight
			TSize matchLen = 0;
			TSize overlapLen = 0;
			TSize alignLen = 0;
			getAlignmentStatistics(matches, pairSet, from, to, matchLen, overlapLen, alignLen);

			if (((matchLen * 100) / overlapLen >= 80) && (matchLen >= 5)) {
				// Create a corresponding edge
				if (seq1<seq2) _getAlignmentStatistics(dist, seq1, seq2, nseq, matchLen, (matchLen * 100) / overlapLen);
				else _getAlignmentStatistics(dist, seq2, seq1, nseq, matchLen, (matchLen * 100) / overlapLen);

				// Record the scores
				resize(scores, to);
				typedef typename Iterator<TScoreValues, Standard>::Type TScoreIter;
				TScoreIter itScore = begin(scores, Standard());
				TScoreIter itScoreEnd = end(scores, Standard());
				itScore+=from;
				for(;itScore != itScoreEnd; ++itScore) *itScore = myScore;

				// Update the overlap counter
				TSize lenLast = matches[from].len; 
				if (matches[to - 1].begin1 == 0) ++frontOvl[seq1];
				if (matches[to - 1].begin2 == 0) ++frontOvl[seq2];
				if (matches[from].begin1 + lenLast == length(pairSet[0])) ++backOvl[seq1];
				if (matches[from].begin2 + lenLast == length(pairSet[1])) ++backOvl[seq2];
			} else resize(matches, from);
		}
	}
}

}// namespace SEQAN_NAMESPACE_MAIN

#endif //#ifndef SEQAN_HEADER_...