fixed bug with realigner for chimera slayer

[mothur.git] / chimerarealigner.cpp

diff --git a/chimerarealigner.cpp b/chimerarealigner.cpp
index 57c3641d14f4f2664dd12c3e1a7f89d13722ab91..bfd283db0d1b66db7df4606d980cc9d54dad58a7 100644 (file)
--- a/chimerarealigner.cpp
+++ b/chimerarealigner.cpp
@@ -12,7 +12,7 @@
 #include "nast.hpp"
 
 //***************************************************************************************************************
-ChimeraReAligner::ChimeraReAligner(vector<Sequence*> t, int m, int mm) : match(m), misMatch(mm) {  templateSeqs = t;  }
+ChimeraReAligner::ChimeraReAligner(vector<Sequence*> t, int ms, int mm) : match(ms), misMatch(mm) {  templateSeqs = t;   m = MothurOut::getInstance(); }
 //***************************************************************************************************************
 ChimeraReAligner::~ChimeraReAligner() {}       
 //***************************************************************************************************************
@@ -24,50 +24,66 @@ void ChimeraReAligner::reAlign(Sequence* query, vector<results> parents) {
                        
                        string qAligned = query->getAligned();
                        string newQuery = "";
-                       
+               
                        //sort parents by region start
                        sort(parents.begin(), parents.end(), compareRegionStart);
 
                        //make sure you don't cutoff beginning of query 
-                       if (parents[0].nastRegionStart > 0) {  newQuery += qAligned.substr(0, parents[0].nastRegionStart+1);  }
+                       if (parents[0].nastRegionStart > 0) {  newQuery += qAligned.substr(0, parents[0].nastRegionStart);  }
                        int longest = 0;
+
                        //take query and break apart into pieces using breakpoints given by results of parents
                        for (int i = 0; i < parents.size(); i++) {
                                int length = parents[i].nastRegionEnd - parents[i].nastRegionStart+1;
                                string q = qAligned.substr(parents[i].nastRegionStart, length);
+       
                                Sequence* queryFrag = new Sequence(query->getName(), q);
-                               
+
                                queryParts.push_back(queryFrag);
-                               
+                       
                                Sequence* parent = getSequence(parents[i].parent);
                                string p = parent->getAligned();
+               
                                p = p.substr(parents[i].nastRegionStart, length);
                                parent->setAligned(p);
                                
                                parentParts.push_back(parent);
-                               
-                               if (q.length() > longest)       { longest = q.length(); }
-                               if (p.length() > longest)       { longest = p.length(); }
+
+                               if (queryFrag->getUnaligned().length() > longest)       { longest = queryFrag->getUnaligned().length(); }
+                               if (parent->getUnaligned().length() > longest)  { longest = parent->getUnaligned().length();    }
                        }
-                                                                                       
+
                        //align each peice to correct parent from results
                        for (int i = 0; i < queryParts.size(); i++) {
-                               alignment = new NeedlemanOverlap(-2.0, match, misMatch, longest+1); //default gapopen, match, mismatch, longestbase
-                               Nast nast(alignment, queryParts[i], parentParts[i]);
-                               delete alignment;
+                               if ((queryParts[i]->getUnaligned() == "") || (parentParts[i]->getUnaligned() == "")) {;}
+                               else {
+                                       Alignment* alignment = new NeedlemanOverlap(-2.0, 1.0, -1.0, longest+1); //default gapopen, match, mismatch, longestbase
+                               
+                                       Nast nast(alignment, queryParts[i], parentParts[i]);
+                                       delete alignment;
+                               }
                        }
-                                                                                               
+
                        //recombine pieces to form new query sequence
                        for (int i = 0; i < queryParts.size(); i++) {
+                               //sometimes the parent regions do not meet, for example region 1 may end at 1000 and region 2 starts at 1100.  
+                               //we don't want to loose length so in this case we will leave query alone
+                               if (i != 0) {
+                                       int space = parents[i].nastRegionStart - parents[i-1].nastRegionEnd - 1;
+                                       if (space > 0) { //they don't meet and we need to add query piece
+                                               string q = qAligned.substr(parents[i-1].nastRegionEnd+1, space);
+                                               newQuery += q;
+                                       }
+                               }
+
                                newQuery += queryParts[i]->getAligned();
                        }
                        
                        //make sure you don't cutoff end of query 
-                       if (parents[parents.size()-1].nastRegionEnd < qAligned.length()) {  newQuery += qAligned.substr(parents[parents.size()-1].nastRegionEnd-1);  }
+                       if (parents[parents.size()-1].nastRegionEnd < (qAligned.length()-1)) {  newQuery += qAligned.substr(parents[parents.size()-1].nastRegionEnd+1);  }
                        
-                       //set query to new aligned string
                        query->setAligned(newQuery);
-                       
+
                        //free memory
                        for (int i = 0; i < queryParts.size(); i++) { delete queryParts[i];  }
                        for (int i = 0; i < parentParts.size(); i++) { delete parentParts[i];  }
@@ -76,7 +92,7 @@ void ChimeraReAligner::reAlign(Sequence* query, vector<results> parents) {
                
        }
        catch(exception& e) {
-               errorOut(e, "ChimeraReAligner", "reAlign");
+               m->errorOut(e, "ChimeraReAligner", "reAlign");
                exit(1);
        }
 }
@@ -94,14 +110,14 @@ Sequence* ChimeraReAligner::getSequence(string name) {
                        }
                }
                
-               if(spot == -1) { mothurOut("Error: Could not find sequence."); mothurOutEndLine(); return NULL; }
+               if(spot == -1) { m->mothurOut("Error: Could not find sequence."); m->mothurOutEndLine(); return NULL; }
                
                temp = new Sequence(templateSeqs[spot]->getName(), templateSeqs[spot]->getAligned());
                
                return temp;
        }
        catch(exception& e) {
-               errorOut(e, "ChimeraReAligner", "getSequence");
+               m->errorOut(e, "ChimeraReAligner", "getSequence");
                exit(1);
        }
 }