18#if !defined(_SMITH_WATERMAN_)
19#define _SMITH_WATERMAN_
88#include "CigarRoller.h"
100#define MAX(x,y) ((x)>(y) ? (x) : (y))
103#define MIN(x,y) ((x)<(y) ? (x) : (y))
126template <
int maxReadLengthH,
int maxReferenceLengthH,
typename HCellType,
typename Atype,
typename Btype,
typename QualityType,
typename readIndexType,
typename referenceIndexType>
155 HCellType H[maxReadLengthH][maxReferenceLengthH];
158 QualityType *qualities;
161 readIndexType MOffset;
162 referenceIndexType NOffset;
164 int allowedInsertDelete;
169 vector<pair<int,int> > alignment;
170 void clearAlignment()
175 HCellType maxCostValue;
176 pair<int,int> maxCostPosition;
187 maxCostPosition.first = 0;
188 maxCostPosition.second = 0;
196 allowedInsertDelete = 0;
213 QualityType *qualities,
217 int allowedInsertDelete = INT_MAX,
219 readIndexType MOffset = 0,
220 referenceIndexType NOffset = 0):
222 qualities(qualities),
226 allowedInsertDelete(allowedInsertDelete),
227 direction(direction),
230 maxCostValue((HCellType) 0)
234 void setRead(Atype *A)
238 void setReadQuality(QualityType *qualities)
240 this->qualities = qualities;
242 void setReference(Btype *B)
250 void setReadLength(
int m)
254 void setReadOffset(readIndexType MOffset)
256 this->MOffset = MOffset;
263 void setReferenceLength(
int n)
267 void setReferenceOffset(referenceIndexType NOffset)
269 this->NOffset = NOffset;
281 void setAllowedInsertDelete(
int allowedInsertDelete = INT_MAX)
283 this->allowedInsertDelete = allowedInsertDelete;
290 void setDirection(
int direction)
292 this->direction = direction;
297 memset(H, 0,
sizeof(H));
305 for (
int i=1; i<=m ; i++)
309 int low = MAX(1, i - allowedInsertDelete);
310 int high = MIN(n, i + allowedInsertDelete);
312 for (
int j=low; j<=high ; j++)
316 if (direction>0) c = MAX(c, H[i-1][j-1] + (((*A)[MOffset + i-1]==(*B)[NOffset + j-1]) ? w.match : w.misMatch));
317 else c = MAX(c, H[i-1][j-1] + (((*A)[MOffset + m-i+0]==(*B)[NOffset + n-j+0]) ? w.match : w.misMatch));
318 c = MAX(c, H[i-1][j] + w.del);
319 c = MAX(c, H[i][j-1] + w.insert);
324 maxCostPosition.first = i;
325 maxCostPosition.second = j;
334 void printH(
bool prettyPrint =
true)
337 for (
int i=-1; i<=m ; i++)
339 for (
int j=-1; j<=n ; j++)
341 if (prettyPrint) cout << setw(3);
344 if (prettyPrint) cout <<
" ";
349 if (!prettyPrint) cout <<
"\t";
350 if (i==0) cout <<
"-";
351 else cout << (*A)[MOffset + direction>0 ? i-1 : m - i];
355 if (!prettyPrint) cout <<
"\t";
356 if (j==0) cout <<
"-";
357 else cout << (*B)[NOffset + direction>0 ? j-1 : n - j];
361 if (!prettyPrint) cout <<
"\t";
369 void debugPrint(
bool doPrintH =
true)
371 if (doPrintH) printH();
372 cout <<
"maxCostPosition = " << maxCostPosition << std::endl;
373 if (alignment.empty()) cout <<
"alignment vector is empty.\n";
376 cout <<
"alignment vector:\n";
377 for (vector<pair<int,int> >::iterator i=alignment.begin(); i < alignment.end(); i++)
379 cout << (i - alignment.begin()) <<
": " << *i <<
"\n";
389 void populateAlignment()
394 i = maxCostPosition.first;
395 j = maxCostPosition.second;
403 while (H[i][j] > 0 || (i>0 && j>0))
406#if defined(DEBUG_ALIGNMENT_VECTOR)
407 cout <<
"alignment.push_back(" << i <<
", " << j <<
")" << endl;
409 alignment.push_back(pair<int,int>(i,j));
410 if (H[i-1][j-1]>=H[i-1][j] && H[i-1][j-1]>=H[i][j-1])
416 else if (H[i-1][j] < H[i][j-1])
427 alignment.push_back(pair<int,int>(i,j));
428#if defined(DEBUG_ALIGNMENT_VECTOR)
429 cout <<
"alignment.push_back(" << i <<
", " << j <<
")" << endl;
430 cout <<
"alignment.size(): " << alignment.size() << endl;
450 if (direction>0)
return getSumQForward();
451 else return getSumQBackward();
457 vector<pair<int,int> >::reverse_iterator i;
459 for (i=alignment.rbegin(); i < alignment.rend() - 1; i++)
462#if defined(DEBUG_GETSUMQ)
465 if ((*(i+1)).first == ((*i).first+1) && (*(i+1)).second == ((*i).second + 1))
468#if defined(DEBUG_GETSUMQ)
469 cout <<
"Match/Mismatch";
471 if ((*A)[MOffset + (*i).first] != (*B)[NOffset + (*i).second])
472 sumQ += (*qualities)[MOffset + (*i).first] -
'!';
474 else if ((*(i+1)).first == ((*i).first+1) && (*(i+1)).second == ((*i).second))
477#if defined(DEBUG_GETSUMQ)
482 else if ((*(i+1)).first == ((*i).first) && (*(i+1)).second == ((*i).second + 1))
485#if defined(DEBUG_GETSUMQ)
491#if defined(DEBUG_GETSUMQ)
497 int getSumQBackward()
500 vector<pair<int,int> >::iterator i;
502 for (i=alignment.begin(); i < alignment.end() - 1; i++)
504#if defined(DEBUG_GETSUMQ)
507 if ((*(i+1)).first == ((*i).first-1) && (*(i+1)).second == ((*i).second - 1))
510#if defined(DEBUG_GETSUMQ)
511 cout <<
"Match/Mismatch";
513 if ((*A)[MOffset + m - (*i).first] != (*B)[NOffset + n - (*i).second])
514 sumQ += (*qualities)[MOffset + m - (*i).first] -
'!';
516 else if ((*(i+1)).first == ((*i).first-1) && (*(i+1)).second == ((*i).second))
519#if defined(DEBUG_GETSUMQ)
524 else if ((*(i+1)).first == ((*i).first) && (*(i+1)).second == ((*i).second - 1))
527#if defined(DEBUG_GETSUMQ)
533#if defined(DEBUG_GETSUMQ)
542 vector<pair<int,int> >::reverse_iterator i;
544 for (i=alignment.rbegin(); i < alignment.rend() - 1; i++)
546#if defined(DEBUG_ALIGNMENT_VECTOR)
547 cout <<
"i: " << i - alignment.rbegin() << *i << endl;
552 if ((*A)[MOffset + (*i).first] != (*B)[NOffset + (*i).second])
553 sumQ += (*qualities)[MOffset + (*i).first] -
'!';
558 if ((*A)[MOffset + m - (*i).first] != (*B)[NOffset + n - (*i).second])
559 sumQ += (*qualities)[MOffset + m - (*i).first] -
'!';
576 vector<pair<int,int> >::reverse_iterator i;
578 for (i=alignment.rbegin(); i < alignment.rend() - 1; i++)
581#if defined(DEBUG_CIGAR)
584 if ((*(i+1)).first == ((*i).first+1) && (*(i+1)).second == ((*i).second + 1))
587#if defined(DEBUG_CIGAR)
588 cout <<
"Match/Mismatch";
592 else if ((*(i+1)).first == ((*i).first+1) && (*(i+1)).second == ((*i).second))
595#if defined(DEBUG_CIGAR)
600 else if ((*(i+1)).first == ((*i).first) && (*(i+1)).second == ((*i).second + 1))
603#if defined(DEBUG_CIGAR)
612#if defined(DEBUG_CIGAR)
627 vector<pair<int,int> >::iterator i;
633 i = alignment.begin();
635 for (i=alignment.begin();
636 i < alignment.end() - 1;
639#if defined(DEBUG_CIGAR)
642 if ((*(i+1)).first == ((*i).first-1) && (*(i+1)).second == ((*i).second - 1))
645#if defined(DEBUG_CIGAR)
646 cout <<
"Match/Mismatch";
650 else if ((*(i+1)).first == ((*i).first-1) && (*(i+1)).second == ((*i).second))
653#if defined(DEBUG_CIGAR)
658 else if ((*(i+1)).first == ((*i).first) && (*(i+1)).second == ((*i).second - 1))
661#if defined(DEBUG_CIGAR)
667#if defined(DEBUG_CIGAR)
680 int getSoftClipCount()
685 return m - maxCostPosition.first;
692 return m - maxCostPosition.first;
698 if (direction>0) rollCigarForward(cigar);
699 else rollCigarBackward(cigar);
716 readIndexType readLength,
717 QualityType &quality,
719 referenceIndexType referenceLength,
720 referenceIndexType referenceOffset,
722 uint32_t &softClipCount,
723 referenceIndexType &cigarStartingPoint,
733 setAllowedInsertDelete(bandSize);
737 setReadLength(readLength);
739 setReadQuality(&quality);
741 setReference(&reference);
742 setReferenceOffset(referenceOffset);
743 setReferenceLength(referenceLength);
747 softClipCount = getSoftClipCount();
751 rollCigar(cigarRoller);
The purpose of this class is to provide accessors for setting, updating, modifying the CIGAR object....
void Add(Operation operation, int count)
Append the specified operation with the specified count to this object.
void clear()
Clear this object so that it has no Cigar Operations.
@ del
deletion from the reference (the reference contains bases that have no corresponding base in the quer...
@ match
match/mismatch operation. Associated with CIGAR Operation "M"
@ insert
insertion to the reference (the query sequence contains bases that have no corresponding base in the ...
@ softClip
Soft clip on the read (clipped sequence present in the query sequence, but not in reference)....