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);