4 * Created by Pat Schloss on 12/15/08.
5 * Copyright 2008 Patrick D. Schloss. All rights reserved.
7 * This is a class for an abstract datatype for classes that implement various types of alignment algorithms.
8 * As of 12/18/08 these included alignments based on blastn, needleman-wunsch, and the Gotoh algorithms
17 #include "alignmentcell.hpp"
18 #include "alignment.hpp"
22 /**************************************************************************************************/
24 Alignment::Alignment() { /* do nothing */ }
26 /**************************************************************************************************/
28 Alignment::Alignment(int A) : nCols(A), nRows(A) {
30 alignment.resize(nRows); // For the Gotoh and Needleman-Wunsch we initialize the dynamic programming
31 for(int i=0;i<nRows;i++){ // matrix by initializing a matrix that is A x A. By default we will set A
32 alignment[i].resize(nCols); // at 2000 for 16S rRNA gene sequences
37 /**************************************************************************************************/
39 void Alignment::traceBack(){ // This traceback routine is used by the dynamic programming algorithms
40 // to fill the values of seqAaln and seqBaln
48 AlignmentCell currentCell = alignment[row][column]; // Start the traceback from the bottom-right corner of the
51 if(currentCell.prevCell == 'x'){ seqAaln = seqBaln = "NOALIGNMENT"; }//If there's an 'x' in the bottom-
52 else{ // right corner bail out because it means nothing got aligned
53 while(currentCell.prevCell != 'x'){ // while the previous cell isn't an 'x', keep going...
55 if(currentCell.prevCell == 'u'){ // if the pointer to the previous cell is 'u', go up in the
56 seqAaln = '-' + seqAaln; // matrix. this indicates that we need to insert a gap in
57 seqBaln = seqB[row] + seqBaln; // seqA and a base in seqB
58 currentCell = alignment[--row][column];
60 else if(currentCell.prevCell == 'l'){ // if the pointer to the previous cell is 'l', go to the left
61 seqBaln = '-' + seqBaln; // in the matrix. this indicates that we need to insert a gap
62 seqAaln = seqA[column] + seqAaln; // in seqB and a base in seqA
63 currentCell = alignment[row][--column];
66 seqAaln = seqA[column] + seqAaln; // otherwise we need to go diagonally up and to the left,
67 seqBaln = seqB[row] + seqBaln; // here we add a base to both alignments
68 currentCell = alignment[--row][--column];
73 pairwiseLength = seqAaln.length();
74 seqAstart = 1; seqAend = 0;
75 seqBstart = 1; seqBend = 0;
77 for(int i=0;i<seqAaln.length();i++){
78 if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAstart++; }
79 else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBstart++; }
83 pairwiseLength -= (seqAstart + seqBstart - 2);
85 for(int i=seqAaln.length()-1; i>=0;i--){
86 if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAend++; }
87 else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBend++; }
90 pairwiseLength -= (seqAend + seqBend);
92 seqAend = seqA.length() - seqAend - 1;
93 seqBend = seqB.length() - seqBend - 1;
97 /**************************************************************************************************/
99 string Alignment::getSeqAAln(){
100 return seqAaln; // this is called to get the alignment of seqA
103 /**************************************************************************************************/
105 string Alignment::getSeqBAln(){
106 return seqBaln; // this is called to get the alignment of seqB
109 /**************************************************************************************************/
111 int Alignment::getCandidateStartPos(){
112 return seqAstart; // this is called to report the quality of the alignment
115 /**************************************************************************************************/
117 int Alignment::getCandidateEndPos(){
118 return seqAend; // this is called to report the quality of the alignment
121 /**************************************************************************************************/
123 int Alignment::getTemplateStartPos(){
124 return seqBstart; // this is called to report the quality of the alignment
127 /**************************************************************************************************/
129 int Alignment::getTemplateEndPos(){
130 return seqBend; // this is called to report the quality of the alignment
133 /**************************************************************************************************/
135 int Alignment::getPairwiseLength(){
136 return pairwiseLength; // this is the pairwise alignment length
139 /**************************************************************************************************/
141 //int Alignment::getLongestTemplateGap(){
143 // int length = seqBaln.length();
145 // int gapLength = 0;
147 // int start = seqAstart;
148 // if(seqAstart < seqBstart){ start = seqBstart; }
149 // for(int i=seqAstart;i<length;i++){
150 // if(seqBaln[i] == '-'){
154 // if(gapLength > 0){
155 // if(gapLength > longGap){ longGap = gapLength; }
163 /**************************************************************************************************/