/**************************************************************************************************/
-Alignment::Alignment() { /* do nothing */ }
+Alignment::Alignment() { m = MothurOut::getInstance(); /* do nothing */ }
/**************************************************************************************************/
Alignment::Alignment(int A) : nCols(A), nRows(A) {
-
- 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
- }
-
+ 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) {
+ 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
- // to fill the values of seqAaln and seqBaln
- seqAaln = "";
- seqBaln = "";
- int row = lB-1;
- int column = lA-1;
-// seqAstart = 1;
-// seqAend = column;
-
- AlignmentCell currentCell = alignment[row][column]; // Start the traceback from the bottom-right corner of the
- // 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
- 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
- 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
- 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
- currentCell = alignment[--row][--column];
+ try {
+ BBaseMap.clear();
+ ABaseMap.clear(); // to fill the values of seqAaln and seqBaln
+ seqAaln = "";
+ seqBaln = "";
+ int row = lB-1;
+ int column = lA-1;
+ // seqAstart = 1;
+ // seqAend = column;
+
+ AlignmentCell currentCell = alignment[row][column]; // Start the traceback from the bottom-right corner of the
+ // 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
+ 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();
+ 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++; }
+ else { break; }
+ }
+
+ pairwiseLength -= (seqAstart + seqBstart - 2);
+
+ for(int i=seqAaln.length()-1; i>=0;i--){
+ if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAend++; }
+ else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBend++; }
+ else { break; }
+ }
+ pairwiseLength -= (seqAend + seqBend);
+
+ seqAend = seqA.length() - seqAend - 1;
+ seqBend = seqB.length() - seqBend - 1;
}
-
- pairwiseLength = seqAaln.length();
- seqAstart = 1; seqAend = 0;
- seqBstart = 1; seqBend = 0;
-
- for(int i=0;i<seqAaln.length();i++){
- if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAstart++; }
- else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBstart++; }
- else { break; }
- }
-
- pairwiseLength -= (seqAstart + seqBstart - 2);
-
- for(int i=seqAaln.length()-1; i>=0;i--){
- if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAend++; }
- else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBend++; }
- else { break; }
+ catch(exception& e) {
+ m->errorOut(e, "Alignment", "traceBack");
+ exit(1);
}
- pairwiseLength -= (seqAend + seqBend);
-
- seqAend = seqA.length() - seqAend - 1;
- seqBend = seqB.length() - seqBend - 1;
-
}
/**************************************************************************************************/
}
}
catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the Alignment class Function ~Alignment. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
- exit(1);
- }
- catch(...) {
- cout << "An unknown error has occurred in the Alignment class function ~Alignment. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ 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(){
projects / mothur.git / blobdiff