X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=maligner.cpp;h=1d5258964650cad8c2e90d39414022e7dc2f0455;hb=b20d7d9b1b0eeeaedb78d8fdf26833fd212d2e81;hp=b393194bbc3a4c89fd2ae8ea3f7f7ba67f84b927;hpb=9ada98592a54c82d08f3d46c9b1d8c3e472a922d;p=mothur.git diff --git a/maligner.cpp b/maligner.cpp index b393194..1d52589 100644 --- a/maligner.cpp +++ b/maligner.cpp @@ -8,12 +8,12 @@ */ #include "maligner.h" -#include "database.hpp" +#include "kmerdb.hpp" #include "blastdb.hpp" /***********************************************************************/ -Maligner::Maligner(vector temp, int num, int match, int misMatch, float div, int ms, int minCov, string mode) : - db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov), searchMethod(mode) {} +Maligner::Maligner(vector temp, int num, int match, int misMatch, float div, int ms, int minCov, string mode, Database* dataLeft, Database* dataRight) : + db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov), searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight) { m = MothurOut::getInstance(); } /***********************************************************************/ string Maligner::getResults(Sequence* q, DeCalculator* decalc) { try { @@ -25,15 +25,19 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) { string chimera; - if (searchMethod != "blast") { + if (searchMethod == "distance") { //find closest seqs to query in template - returns copies of seqs so trim does not destroy - remember to deallocate - refSeqs = decalc->findClosest(query, db, numWanted); - }else{ - refSeqs = getBlastSeqs(query, numWanted); - } + refSeqs = decalc->findClosest(query, db, numWanted, indexes); + }else if (searchMethod == "blast") { + refSeqs = getBlastSeqs(query, numWanted); //fills indexes + }else if (searchMethod == "kmer") { + refSeqs = getKmerSeqs(query, numWanted); //fills indexes + }else { m->mothurOut("not valid search."); exit(1); } //should never get here - refSeqs = minCoverageFilter(refSeqs); + if (m->control_pressed) { return chimera; } + refSeqs = minCoverageFilter(refSeqs); + if (refSeqs.size() < 2) { for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } percentIdenticalQueryChimera = 0.0; @@ -41,19 +45,21 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) { } int chimeraPenalty = computeChimeraPenalty(); - + //fills outputResults chimera = chimeraMaligner(chimeraPenalty, decalc); + if (m->control_pressed) { return chimera; } //free memory delete query; + for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } return chimera; } catch(exception& e) { - errorOut(e, "Maligner", "getResults"); + m->errorOut(e, "Maligner", "getResults"); exit(1); } } @@ -65,33 +71,39 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { //trims seqs to first non gap char in all seqs and last non gap char in all seqs spotMap = decalc->trimSeqs(query, refSeqs); - + vector temp = refSeqs; temp.push_back(query); verticalFilter(temp); +//for (int i = 0; i < temp.size(); i++) { cout << temp[i]->getName() << '\n' << temp[i]->getAligned() << endl; } return "no"; vector< vector > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes - fillScoreMatrix(matrix, refSeqs, chimeraPenalty); + if (m->control_pressed) { return chimera; } + fillScoreMatrix(matrix, refSeqs, chimeraPenalty); + vector path = extractHighestPath(matrix); - vector trace = mapTraceRegionsToAlignment(path, refSeqs); + if (m->control_pressed) { return chimera; } + vector trace = mapTraceRegionsToAlignment(path, refSeqs); + if (trace.size() > 1) { chimera = "yes"; } else { chimera = "no"; } int traceStart = path[0].col; - int traceEnd = path[path.size()-1].col; - + int traceEnd = path[path.size()-1].col; string queryInRange = query->getAligned(); queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1)); string chimeraSeq = constructChimericSeq(trace, refSeqs); - + percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq); + if (m->control_pressed) { return chimera; } + //save output results for (int i = 0; i < trace.size(); i++) { int regionStart = trace[i].col; @@ -101,6 +113,7 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { results temp; temp.parent = refSeqs[seqIndex]->getName(); + temp.parentAligned = db[indexes[seqIndex]]->getAligned(); temp.nastRegionStart = spotMap[regionStart]; temp.nastRegionEnd = spotMap[regionEnd]; temp.regionStart = regionStart; @@ -126,7 +139,7 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { return chimera; } catch(exception& e) { - errorOut(e, "Maligner", "chimeraMaligner"); + m->errorOut(e, "Maligner", "chimeraMaligner"); exit(1); } } @@ -168,7 +181,7 @@ vector Maligner::minCoverageFilter(vector ref){ return newRefs; } catch(exception& e) { - errorOut(e, "Maligner", "minCoverageFilter"); + m->errorOut(e, "Maligner", "minCoverageFilter"); exit(1); } } @@ -185,7 +198,7 @@ int Maligner::computeChimeraPenalty() { } catch(exception& e) { - errorOut(e, "Maligner", "computeChimeraPenalty"); + m->errorOut(e, "Maligner", "computeChimeraPenalty"); exit(1); } } @@ -237,7 +250,7 @@ void Maligner::verticalFilter(vector seqs) { spotMap = newMap; } catch(exception& e) { - errorOut(e, "Maligner", "verticalFilter"); + m->errorOut(e, "Maligner", "verticalFilter"); exit(1); } } @@ -264,12 +277,12 @@ vector< vector > Maligner::buildScoreMatrix(int cols, int rows) { return m; } catch(exception& e) { - errorOut(e, "Maligner", "buildScoreMatrix"); + m->errorOut(e, "Maligner", "buildScoreMatrix"); exit(1); } } //*************************************************************************************************************** -void Maligner::fillScoreMatrix(vector >& m, vector seqs, int penalty) { +void Maligner::fillScoreMatrix(vector >& ms, vector seqs, int penalty) { try{ //get matrix dimensions @@ -283,11 +296,11 @@ void Maligner::fillScoreMatrix(vector >& m, vector >& m, vector m[i][j].score) { - m[i][j].score = sumScore; - m[i][j].prev = prevIndex; + if (sumScore > ms[i][j].score) { + ms[i][j].score = sumScore; + ms[i][j].prev = prevIndex; } } } @@ -329,17 +342,17 @@ void Maligner::fillScoreMatrix(vector >& m, vectorerrorOut(e, "Maligner", "fillScoreMatrix"); exit(1); } } //*************************************************************************************************************** -vector Maligner::extractHighestPath(vector > m) { +vector Maligner::extractHighestPath(vector > ms) { try { //get matrix dimensions int numCols = query->getAligned().length(); - int numRows = m.size(); + int numRows = ms.size(); //find highest score scoring matrix @@ -348,9 +361,9 @@ vector Maligner::extractHighestPath(vector > for (int i = 0; i < numRows; i++) { for (int j = 0; j < numCols; j++) { - if (m[i][j].score > highestScore) { - highestScore = m[i][j].score; - highestStruct = m[i][j]; + if (ms[i][j].score > highestScore) { + highestScore = ms[i][j].score; + highestStruct = ms[i][j]; } } } @@ -362,7 +375,7 @@ vector Maligner::extractHighestPath(vector > int score = highestStruct.score; while (pos >= 0 && score > 0) { - score_struct temp = m[rowIndex][pos]; + score_struct temp = ms[rowIndex][pos]; score = temp.score; if (score > 0) { path.push_back(temp); } @@ -377,7 +390,7 @@ vector Maligner::extractHighestPath(vector > } catch(exception& e) { - errorOut(e, "Maligner", "extractHighestPath"); + m->errorOut(e, "Maligner", "extractHighestPath"); exit(1); } } @@ -420,7 +433,7 @@ vector Maligner::mapTraceRegionsToAlignment(vector p } catch(exception& e) { - errorOut(e, "Maligner", "mapTraceRegionsToAlignment"); + m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment"); exit(1); } } @@ -438,7 +451,7 @@ string Maligner::constructChimericSeq(vector trace, vectorerrorOut(e, "Maligner", "constructChimericSeq"); exit(1); } } @@ -447,9 +460,9 @@ float Maligner::computePercentID(string queryAlign, string chimera) { try { if (queryAlign.length() != chimera.length()) { - mothurOut("Error, alignment strings are of different lengths: "); mothurOutEndLine(); - mothurOut(toString(queryAlign.length())); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine(); - mothurOut(toString(chimera.length())); mothurOutEndLine(); + m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine(); + m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine(); + m->mothurOut(toString(chimera.length())); m->mothurOutEndLine(); return -1.0; } @@ -474,18 +487,90 @@ float Maligner::computePercentID(string queryAlign, string chimera) { } catch(exception& e) { - errorOut(e, "Maligner", "computePercentID"); + m->errorOut(e, "Maligner", "computePercentID"); exit(1); } } //*************************************************************************************************************** vector Maligner::getBlastSeqs(Sequence* q, int num) { try { - //generate blastdb - Database* database = new BlastDB(-2.0, -1.0, matchScore, misMatchPenalty); - for (int i = 0; i < db.size(); i++) { database->addSequence(*db[i]); } - database->generateDB(); - database->setNumSeqs(db.size()); + indexes.clear(); + vector refResults; + + //get parts of query + string queryUnAligned = q->getUnaligned(); + string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence + string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence + + Sequence* queryLeft = new Sequence(q->getName()+"left", leftQuery); + Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery); + + vector tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1); + vector tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1); + + //if ((tempIndexesRight.size() == 0) && (tempIndexesLeft.size() == 0)) { m->mothurOut("megablast returned " + toString(tempIndexesRight.size()) + " results for the right end, and " + toString(tempIndexesLeft.size()) + " for the left end. Needed " + toString(num+1) + ". Unable to process sequence " + q->getName()); m->mothurOutEndLine(); return refResults; } + + vector smaller; + vector larger; + + if (tempIndexesRight.size() < tempIndexesLeft.size()) { smaller = tempIndexesRight; larger = tempIndexesLeft; } + else { smaller = tempIndexesLeft; larger = tempIndexesRight; } + + //merge results + map seen; + map::iterator it; + + vector mergedResults; + for (int i = 0; i < smaller.size(); i++) { + //add left if you havent already + it = seen.find(smaller[i]); + if (it == seen.end()) { + mergedResults.push_back(smaller[i]); + seen[smaller[i]] = smaller[i]; + } + + //add right if you havent already + it = seen.find(larger[i]); + if (it == seen.end()) { + mergedResults.push_back(larger[i]); + seen[larger[i]] = larger[i]; + } + } + + for (int i = smaller.size(); i < larger.size(); i++) { + it = seen.find(larger[i]); + if (it == seen.end()) { + mergedResults.push_back(larger[i]); + seen[larger[i]] = larger[i]; + } + } + + if (mergedResults.size() < numWanted) { numWanted = mergedResults.size(); } +//cout << q->getName() << endl; + for (int i = 0; i < numWanted; i++) { +//cout << db[mergedResults[i]]->getName() << endl; + if (db[mergedResults[i]]->getName() != q->getName()) { + Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); + refResults.push_back(temp); + indexes.push_back(mergedResults[i]); + } +//cout << mergedResults[i] << endl; + } +//cout << endl; + delete queryRight; + delete queryLeft; + + return refResults; + } + catch(exception& e) { + m->errorOut(e, "Maligner", "getBlastSeqs"); + exit(1); + } +} +//*************************************************************************************************************** +vector Maligner::getKmerSeqs(Sequence* q, int num) { + try { + indexes.clear(); //get parts of query string queryUnAligned = q->getUnaligned(); @@ -495,8 +580,8 @@ vector Maligner::getBlastSeqs(Sequence* q, int num) { Sequence* queryLeft = new Sequence(q->getName(), leftQuery); Sequence* queryRight = new Sequence(q->getName(), rightQuery); - vector tempIndexesRight = database->findClosestMegaBlast(queryRight, num+1); - vector tempIndexesLeft = database->findClosestMegaBlast(queryLeft, num+1); + vector tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted); + vector tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted); //merge results map seen; @@ -519,24 +604,24 @@ vector Maligner::getBlastSeqs(Sequence* q, int num) { } } - +//cout << q->getName() << endl; vector refResults; for (int i = 0; i < numWanted; i++) { +//cout << db[mergedResults[i]]->getName() << endl; Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); refResults.push_back(temp); + indexes.push_back(mergedResults[i]); } - +//cout << endl; delete queryRight; delete queryLeft; - delete database; return refResults; } catch(exception& e) { - errorOut(e, "Maligner", "getBlastSeqs"); + m->errorOut(e, "Maligner", "getBlastSeqs"); exit(1); } } - //***************************************************************************************************************