--- /dev/null
+/*
+ * alignment.cpp
+ *
+ * Created by Pat Schloss on 12/15/08.
+ * Copyright 2008 Patrick D. Schloss. All rights reserved.
+ *
+ * This is a class for an abstract datatype for classes that implement various types of alignment algorithms.
+ * As of 12/18/08 these included alignments based on blastn, needleman-wunsch, and the Gotoh algorithms
+ *
+ */
+
+#include "alignmentcell.hpp"
+#include "alignment.hpp"
+
+
+/**************************************************************************************************/
+
+Alignment::Alignment() { /* 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) {
+ 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
+ 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];
+ }
+ }
+ }
+
+ 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; }
+ }
+ pairwiseLength -= (seqAend + seqBend);
+
+ seqAend = seqA.length() - seqAend - 1;
+ seqBend = seqB.length() - seqBend - 1;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "Alignment", "traceBack");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+
+Alignment::~Alignment(){
+ try {
+ for (int i = 0; i < alignment.size(); i++) {
+ for (int j = (alignment[i].size()-1); j >= 0; j--) { alignment[i].pop_back(); }
+ }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "Alignment", "~Alignment");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+string Alignment::getSeqAAln(){
+ return seqAaln; // this is called to get the alignment of seqA
+}
+
+/**************************************************************************************************/
+
+string Alignment::getSeqBAln(){
+ return seqBaln; // this is called to get the alignment of seqB
+}
+
+/**************************************************************************************************/
+
+int Alignment::getCandidateStartPos(){
+ return seqAstart; // this is called to report the quality of the alignment
+}
+
+/**************************************************************************************************/
+
+int Alignment::getCandidateEndPos(){
+ return seqAend; // this is called to report the quality of the alignment
+}
+
+/**************************************************************************************************/
+
+int Alignment::getTemplateStartPos(){
+ return seqBstart; // this is called to report the quality of the alignment
+}
+
+/**************************************************************************************************/
+
+int Alignment::getTemplateEndPos(){
+ return seqBend; // this is called to report the quality of the alignment
+}
+
+/**************************************************************************************************/
+
+int Alignment::getPairwiseLength(){
+ return pairwiseLength; // this is the pairwise alignment length
+}
+
+/**************************************************************************************************/
+
+//int Alignment::getLongestTemplateGap(){
+//
+// int length = seqBaln.length();
+// int longGap = 0;
+// int gapLength = 0;
+//
+// int start = seqAstart;
+// if(seqAstart < seqBstart){ start = seqBstart; }
+// for(int i=seqAstart;i<length;i++){
+// if(seqBaln[i] == '-'){
+// gapLength++;
+// }
+// else{
+// if(gapLength > 0){
+// if(gapLength > longGap){ longGap = gapLength; }
+// }
+// gapLength = 0;
+// }
+// }
+// return longGap;
+//}
+
+/**************************************************************************************************/
projects / mothur.git / blobdiff