/**************************************************************************************************/
-Alignment::Alignment() { /* do nothing */ }
+Alignment::Alignment() { m = MothurOut::getInstance(); /* do nothing */ }
/**************************************************************************************************/
Alignment::Alignment(int A) : nCols(A), nRows(A) {
try {
+
+ m = MothurOut::getInstance();
alignment.resize(nRows); // For the Gotoh and Needleman-Wunsch we initialize the dynamic programming
for(int i=0;i<nRows;i++){ // matrix by initializing a matrix that is A x A. By default we will set A
alignment[i].resize(nCols); // at 2000 for 16S rRNA gene sequences
}
}
catch(exception& e) {
- errorOut(e, "Alignment", "Alignment");
+ m->errorOut(e, "Alignment", "Alignment");
exit(1);
}
}
+/**************************************************************************************************/
+void Alignment::resize(int A) {
+ try {
+ nCols = A;
+ nRows = A;
+ alignment.resize(nRows);
+ for(int i=0;i<nRows;i++){
+ alignment[i].resize(nCols);
+ }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "Alignment", "resize");
+ exit(1);
+ }
+}
/**************************************************************************************************/
void Alignment::traceBack(){ // This traceback routine is used by the dynamic programming algorithms
try {
- // to fill the values of seqAaln and seqBaln
+ BBaseMap.clear();
+ ABaseMap.clear(); // to fill the values of seqAaln and seqBaln
seqAaln = "";
seqBaln = "";
int row = lB-1;
// matrix
if(currentCell.prevCell == 'x'){ seqAaln = seqBaln = "NOALIGNMENT"; }//If there's an 'x' in the bottom-
- else{ // right corner bail out because it means nothing got aligned
+ else{ // right corner bail out because it means nothing got aligned
+ int count = 0;
while(currentCell.prevCell != 'x'){ // while the previous cell isn't an 'x', keep going...
if(currentCell.prevCell == 'u'){ // if the pointer to the previous cell is 'u', go up in the
seqAaln = '-' + seqAaln; // matrix. this indicates that we need to insert a gap in
seqBaln = seqB[row] + seqBaln; // seqA and a base in seqB
+ BBaseMap[row] = count;
currentCell = alignment[--row][column];
}
else if(currentCell.prevCell == 'l'){ // if the pointer to the previous cell is 'l', go to the left
seqBaln = '-' + seqBaln; // in the matrix. this indicates that we need to insert a gap
seqAaln = seqA[column] + seqAaln; // in seqB and a base in seqA
+ ABaseMap[column] = count;
currentCell = alignment[row][--column];
}
else{
seqAaln = seqA[column] + seqAaln; // otherwise we need to go diagonally up and to the left,
seqBaln = seqB[row] + seqBaln; // here we add a base to both alignments
+ BBaseMap[row] = count;
+ ABaseMap[column] = count;
currentCell = alignment[--row][--column];
}
+ count++;
}
}
- pairwiseLength = seqAaln.length();
+
+ pairwiseLength = seqAaln.length();
seqAstart = 1; seqAend = 0;
seqBstart = 1; seqBend = 0;
-
+ //flip maps since we now know the total length
+ map<int, int> newAMap;
+ for (map<int, int>::iterator it = ABaseMap.begin(); it != ABaseMap.end(); it++) {
+ int spot = it->second;
+ newAMap[pairwiseLength-spot-1] = it->first-1;
+ }
+ ABaseMap = newAMap;
+ map<int, int> newBMap;
+ for (map<int, int>::iterator it = BBaseMap.begin(); it != BBaseMap.end(); it++) {
+ int spot = it->second;
+ newBMap[pairwiseLength-spot-1] = it->first-1;
+ }
+ BBaseMap = newBMap;
+
for(int i=0;i<seqAaln.length();i++){
if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAstart++; }
else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBstart++; }
seqBend = seqB.length() - seqBend - 1;
}
catch(exception& e) {
- errorOut(e, "Alignment", "traceBack");
+ m->errorOut(e, "Alignment", "traceBack");
exit(1);
}
}
}
}
catch(exception& e) {
- errorOut(e, "Alignment", "~Alignment");
+ m->errorOut(e, "Alignment", "~Alignment");
exit(1);
}
}
int Alignment::getTemplateStartPos(){
return seqBstart; // this is called to report the quality of the alignment
}
+/**************************************************************************************************/
+map<int, int> Alignment::getSeqAAlnBaseMap(){
+ return ABaseMap;
+}
+/**************************************************************************************************/
+
+map<int, int> Alignment::getSeqBAlnBaseMap(){
+ return BBaseMap;
+}
/**************************************************************************************************/
int Alignment::getTemplateEndPos(){