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
12 #include "alignmentcell.hpp"
13 #include "alignment.hpp"
16 /**************************************************************************************************/
18 Alignment::Alignment() { m = MothurOut::getInstance(); /* do nothing */ }
20 /**************************************************************************************************/
22 Alignment::Alignment(int A) : nCols(A), nRows(A) {
25 m = MothurOut::getInstance();
26 alignment.resize(nRows); // For the Gotoh and Needleman-Wunsch we initialize the dynamic programming
27 for(int i=0;i<nRows;i++){ // matrix by initializing a matrix that is A x A. By default we will set A
28 alignment[i].resize(nCols); // at 2000 for 16S rRNA gene sequences
32 m->errorOut(e, "Alignment", "Alignment");
36 /**************************************************************************************************/
37 void Alignment::resize(int A) {
42 alignment.resize(nRows);
43 for(int i=0;i<nRows;i++){
44 alignment[i].resize(nCols);
48 m->errorOut(e, "Alignment", "resize");
52 /**************************************************************************************************/
54 void Alignment::traceBack(){ // This traceback routine is used by the dynamic programming algorithms
57 ABaseMap.clear(); // to fill the values of seqAaln and seqBaln
65 AlignmentCell currentCell = alignment[row][column]; // Start the traceback from the bottom-right corner of the
68 if(currentCell.prevCell == 'x'){ seqAaln = seqBaln = "NOALIGNMENT"; }//If there's an 'x' in the bottom-
69 else{ // right corner bail out because it means nothing got aligned
71 while(currentCell.prevCell != 'x'){ // while the previous cell isn't an 'x', keep going...
73 if(currentCell.prevCell == 'u'){ // if the pointer to the previous cell is 'u', go up in the
74 seqAaln = '-' + seqAaln; // matrix. this indicates that we need to insert a gap in
75 seqBaln = seqB[row] + seqBaln; // seqA and a base in seqB
76 BBaseMap[row] = count;
77 currentCell = alignment[--row][column];
79 else if(currentCell.prevCell == 'l'){ // if the pointer to the previous cell is 'l', go to the left
80 seqBaln = '-' + seqBaln; // in the matrix. this indicates that we need to insert a gap
81 seqAaln = seqA[column] + seqAaln; // in seqB and a base in seqA
82 ABaseMap[column] = count;
83 currentCell = alignment[row][--column];
86 seqAaln = seqA[column] + seqAaln; // otherwise we need to go diagonally up and to the left,
87 seqBaln = seqB[row] + seqBaln; // here we add a base to both alignments
88 BBaseMap[row] = count;
89 ABaseMap[column] = count;
90 currentCell = alignment[--row][--column];
97 pairwiseLength = seqAaln.length();
98 seqAstart = 1; seqAend = 0;
99 seqBstart = 1; seqBend = 0;
100 //flip maps since we now know the total length
101 map<int, int> newAMap;
102 for (map<int, int>::iterator it = ABaseMap.begin(); it != ABaseMap.end(); it++) {
103 int spot = it->second;
104 newAMap[pairwiseLength-spot-1] = it->first-1;
107 map<int, int> newBMap;
108 for (map<int, int>::iterator it = BBaseMap.begin(); it != BBaseMap.end(); it++) {
109 int spot = it->second;
110 newBMap[pairwiseLength-spot-1] = it->first-1;
114 for(int i=0;i<seqAaln.length();i++){
115 if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAstart++; }
116 else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBstart++; }
120 pairwiseLength -= (seqAstart + seqBstart - 2);
122 for(int i=seqAaln.length()-1; i>=0;i--){
123 if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAend++; }
124 else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBend++; }
127 pairwiseLength -= (seqAend + seqBend);
129 seqAend = seqA.length() - seqAend - 1;
130 seqBend = seqB.length() - seqBend - 1;
132 catch(exception& e) {
133 m->errorOut(e, "Alignment", "traceBack");
137 /**************************************************************************************************/
139 Alignment::~Alignment(){
141 for (int i = 0; i < alignment.size(); i++) {
142 for (int j = (alignment[i].size()-1); j >= 0; j--) { alignment[i].pop_back(); }
145 catch(exception& e) {
146 m->errorOut(e, "Alignment", "~Alignment");
151 /**************************************************************************************************/
153 string Alignment::getSeqAAln(){
154 return seqAaln; // this is called to get the alignment of seqA
157 /**************************************************************************************************/
159 string Alignment::getSeqBAln(){
160 return seqBaln; // this is called to get the alignment of seqB
163 /**************************************************************************************************/
165 int Alignment::getCandidateStartPos(){
166 return seqAstart; // this is called to report the quality of the alignment
169 /**************************************************************************************************/
171 int Alignment::getCandidateEndPos(){
172 return seqAend; // this is called to report the quality of the alignment
175 /**************************************************************************************************/
177 int Alignment::getTemplateStartPos(){
178 return seqBstart; // this is called to report the quality of the alignment
180 /**************************************************************************************************/
182 map<int, int> Alignment::getSeqAAlnBaseMap(){
185 /**************************************************************************************************/
187 map<int, int> Alignment::getSeqBAlnBaseMap(){
190 /**************************************************************************************************/
192 int Alignment::getTemplateEndPos(){
193 return seqBend; // this is called to report the quality of the alignment
196 /**************************************************************************************************/
198 int Alignment::getPairwiseLength(){
199 return pairwiseLength; // this is the pairwise alignment length
202 /**************************************************************************************************/
204 //int Alignment::getLongestTemplateGap(){
206 // int length = seqBaln.length();
208 // int gapLength = 0;
210 // int start = seqAstart;
211 // if(seqAstart < seqBstart){ start = seqBstart; }
212 // for(int i=seqAstart;i<length;i++){
213 // if(seqBaln[i] == '-'){
217 // if(gapLength > 0){
218 // if(gapLength > longGap){ longGap = gapLength; }
226 /**************************************************************************************************/