X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=chimerarealigner.cpp;h=a07e43392a1d96727d4870154bccef452119755b;hb=92de7f976371d41441ad41f02ca83af8b43cef5c;hp=fa9dd566c83014124fa302be6d5505102acdddac;hpb=e72551c9cc5542e6a354f0f3e415fea261421d72;p=mothur.git diff --git a/chimerarealigner.cpp b/chimerarealigner.cpp index fa9dd56..a07e433 100644 --- a/chimerarealigner.cpp +++ b/chimerarealigner.cpp @@ -12,7 +12,7 @@ #include "nast.hpp" //*************************************************************************************************************** -ChimeraReAligner::ChimeraReAligner(vector t, int m, int mm) : match(m), misMatch(mm) { templateSeqs = t; } +ChimeraReAligner::ChimeraReAligner(vector t, int ms, int mm) : match(ms), misMatch(mm) { templateSeqs = t; m = MothurOut::getInstance(); } //*************************************************************************************************************** ChimeraReAligner::~ChimeraReAligner() {} //*************************************************************************************************************** @@ -31,12 +31,14 @@ void ChimeraReAligner::reAlign(Sequence* query, vector parents) { //make sure you don't cutoff beginning of query 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); @@ -46,30 +48,43 @@ void ChimeraReAligner::reAlign(Sequence* query, vector parents) { parent->setAligned(p); parentParts.push_back(parent); - + if (q.length() > longest) { longest = q.length(); } if (p.length() > longest) { longest = p.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 = new NeedlemanOverlap(-2.0, match, misMatch, 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); } - + 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]; } @@ -78,7 +93,7 @@ void ChimeraReAligner::reAlign(Sequence* query, vector parents) { } catch(exception& e) { - errorOut(e, "ChimeraReAligner", "reAlign"); + m->errorOut(e, "ChimeraReAligner", "reAlign"); exit(1); } } @@ -96,14 +111,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); } }