working on pam

[mothur.git] / needlemanoverlap.cpp

diff --git a/needlemanoverlap.cpp b/needlemanoverlap.cpp
index 1239beb1a2a9f61d3f8c2cec178161fb36a7b865..1c2ef223539c65a95ca1dac66a823d226b6108b8 100644 (file)
--- a/needlemanoverlap.cpp
+++ b/needlemanoverlap.cpp
@@ -103,6 +103,98 @@ void NeedlemanOverlap::align(string A, string B){
        }
 
 }
+/**************************************************************************************************/
 
+void NeedlemanOverlap::alignPrimer(string A, string B){
+       try {
+        
+               seqA = ' ' + A; lA = seqA.length();             //      algorithm requires a dummy space at the beginning of each string
+               seqB = ' ' + B; lB = seqB.length();             //      algorithm requires a dummy space at the beginning of each string
+        
+               if (lA > nRows) { m->mothurOut("One of your candidate sequences is longer than you longest template sequence. Your longest template sequence is " + toString(nRows) + ". Your candidate is " + toString(lA) + "."); m->mothurOutEndLine();  }
+               
+               for(int i=1;i<lB;i++){                                  //      This code was largely translated from Perl code provided in Ex 3.1
+            
+                       for(int j=1;j<lA;j++){                          //      of the O'Reilly BLAST book.  I found that the example output had a
+                
+                               //      number of errors
+                               float diagonal;
+                               if(isEquivalent(seqB[i],seqA[j]))       {       diagonal = alignment[i-1][j-1].cValue + match;          }
+                               else                                {   diagonal = alignment[i-1][j-1].cValue + mismatch;       }
+                
+                               float up        = alignment[i-1][j].cValue + gap;
+                               float left      = alignment[i][j-1].cValue + gap;
+                               
+                               if(diagonal >= up){
+                                       if(diagonal >= left){
+                                               alignment[i][j].cValue = diagonal;
+                                               alignment[i][j].prevCell = 'd';
+                                       }
+                                       else{
+                                               alignment[i][j].cValue = left;
+                                               alignment[i][j].prevCell = 'l';
+                                       }
+                               }
+                               else{
+                                       if(up >= left){
+                                               alignment[i][j].cValue = up;
+                                               alignment[i][j].prevCell = 'u';
+                                       }
+                                       else{
+                                               alignment[i][j].cValue = left;
+                                               alignment[i][j].prevCell = 'l';
+                                       }
+                               }
+                       }
+               }
+        
+               Overlap over;
+               over.setOverlap(alignment, lA, lB, 0);          //      Fix gaps at the beginning and end of the sequences
+               traceBack();                                                            //      Traceback the alignment to populate seqAaln and seqBaln
+        
+       }
+       catch(exception& e) {
+               m->errorOut(e, "NeedlemanOverlap", "alignPrimer");
+               exit(1);
+       }
+    
+}
+//********************************************************************/
+bool NeedlemanOverlap::isEquivalent(char oligo, char seq){
+       try {
+               
+        bool same = true;
+                                       
+        oligo = toupper(oligo);
+        seq = toupper(seq);
+        
+        if(oligo != seq){
+            if(oligo == 'A' && (seq != 'A' && seq != 'M' && seq != 'R' && seq != 'W' && seq != 'D' && seq != 'H' && seq != 'V'))       {       same = false;   }
+            else if(oligo == 'C' && (seq != 'C' && seq != 'Y' && seq != 'M' && seq != 'S' && seq != 'B' && seq != 'H' && seq != 'V'))       {  same = false;   }
+            else if(oligo == 'G' && (seq != 'G' && seq != 'R' && seq != 'K' && seq != 'S' && seq != 'B' && seq != 'D' && seq != 'V'))       {  same = false;   }
+            else if(oligo == 'T' && (seq != 'T' && seq != 'Y' && seq != 'K' && seq != 'W' && seq != 'B' && seq != 'D' && seq != 'H'))       {  same = false;   }
+            else if((oligo == '.' || oligo == '-'))           {        same = false;   }
+            else if((oligo == 'N' || oligo == 'I') && (seq == 'N'))                         {  same = false;   }
+            else if(oligo == 'R' && (seq != 'A' && seq != 'G'))                        {       same = false;   }
+            else if(oligo == 'Y' && (seq != 'C' && seq != 'T'))                        {       same = false;   }
+            else if(oligo == 'M' && (seq != 'C' && seq != 'A'))                        {       same = false;   }
+            else if(oligo == 'K' && (seq != 'T' && seq != 'G'))                        {       same = false;   }
+            else if(oligo == 'W' && (seq != 'T' && seq != 'A'))                        {       same = false;   }
+            else if(oligo == 'S' && (seq != 'C' && seq != 'G'))                        {       same = false;   }
+            else if(oligo == 'B' && (seq != 'C' && seq != 'T' && seq != 'G'))       {  same = false;   }
+            else if(oligo == 'D' && (seq != 'A' && seq != 'T' && seq != 'G'))       {  same = false;   }
+            else if(oligo == 'H' && (seq != 'A' && seq != 'T' && seq != 'C'))       {  same = false;   }
+            else if(oligo == 'V' && (seq != 'A' && seq != 'C' && seq != 'G'))       {  same = false;   }
+        }
+
+               
+               
+               return same;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "TrimOligos", "countDiffs");
+               exit(1);
+       }
+}
 /**************************************************************************************************/