X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=maligner.cpp;h=76689e9f11bf1ab3547ef9a33053830a7b436d14;hb=74844a60d80c6dd06e3fb02ee9b928424f9019b0;hp=de6004d7e387d0e241dfbc205695d359dbcf6dda;hpb=5a1e62397b91f57d0d3aff635891df04b8999a88;p=mothur.git diff --git a/maligner.cpp b/maligner.cpp index de6004d..76689e9 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,35 +25,17 @@ 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 -//ofstream out; -//string name = toString(numi+1); -//cout << name << endl; -//name = name.substr(name.find_first_of("|")+1); -//cout << name << endl; -//name = name.substr(name.find_first_of("|")+1); -//cout << name << endl; -//name = name.substr(0, name.find_last_of("|")); -//cout << name << endl; -//string filename = name + ".seqsformaligner"; -//openOutputFile(filename, out); -//for (int u = 0; u < refSeqs.size(); u++) { refSeqs[u]->printSequence(out); } -//out.close(); -//filename = name + ".fasta"; -//openOutputFile(filename, out); -//query->printSequence(out); -//out.close(); - -//for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl; } -//cout << "before = " << refSeqs.size() << endl; refSeqs = minCoverageFilter(refSeqs); -//cout << "after = " << refSeqs.size() << endl; + if (refSeqs.size() < 2) { for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } percentIdenticalQueryChimera = 0.0; @@ -61,19 +43,20 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) { } int chimeraPenalty = computeChimeraPenalty(); -//cout << chimeraPenalty << endl; + //fills outputResults chimera = chimeraMaligner(chimeraPenalty, decalc); - + //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); } } @@ -85,31 +68,31 @@ 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); - + vector path = extractHighestPath(matrix); 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); //save output results @@ -121,6 +104,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; @@ -146,7 +130,7 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { return chimera; } catch(exception& e) { - errorOut(e, "Maligner", "chimeraMaligner"); + m->errorOut(e, "Maligner", "chimeraMaligner"); exit(1); } } @@ -188,7 +172,7 @@ vector Maligner::minCoverageFilter(vector ref){ return newRefs; } catch(exception& e) { - errorOut(e, "Maligner", "minCoverageFilter"); + m->errorOut(e, "Maligner", "minCoverageFilter"); exit(1); } } @@ -205,7 +189,7 @@ int Maligner::computeChimeraPenalty() { } catch(exception& e) { - errorOut(e, "Maligner", "computeChimeraPenalty"); + m->errorOut(e, "Maligner", "computeChimeraPenalty"); exit(1); } } @@ -257,7 +241,7 @@ void Maligner::verticalFilter(vector seqs) { spotMap = newMap; } catch(exception& e) { - errorOut(e, "Maligner", "verticalFilter"); + m->errorOut(e, "Maligner", "verticalFilter"); exit(1); } } @@ -284,12 +268,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 @@ -303,11 +287,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; } } } @@ -349,17 +333,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 @@ -368,9 +352,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]; } } } @@ -382,7 +366,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); } @@ -397,7 +381,7 @@ vector Maligner::extractHighestPath(vector > } catch(exception& e) { - errorOut(e, "Maligner", "extractHighestPath"); + m->errorOut(e, "Maligner", "extractHighestPath"); exit(1); } } @@ -440,7 +424,7 @@ vector Maligner::mapTraceRegionsToAlignment(vector p } catch(exception& e) { - errorOut(e, "Maligner", "mapTraceRegionsToAlignment"); + m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment"); exit(1); } } @@ -458,7 +442,7 @@ string Maligner::constructChimericSeq(vector trace, vectorerrorOut(e, "Maligner", "constructChimericSeq"); exit(1); } } @@ -467,9 +451,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; } @@ -494,14 +478,15 @@ 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 { -cout << q->getName() << endl; + try { + indexes.clear(); + vector refResults; //generate blastdb Database* database = new BlastDB(-2.0, -1.0, matchScore, misMatchPenalty); for (int i = 0; i < db.size(); i++) { database->addSequence(*db[i]); } @@ -509,68 +494,131 @@ cout << q->getName() << endl; database->setNumSeqs(db.size()); //get parts of query - string queryAligned = q->getAligned(); - string leftQuery = queryAligned.substr(0, (queryAligned.length() / 3)); //first 1/3 of the sequence - string rightQuery = queryAligned.substr(((queryAligned.length() / 3)*2)); //last 1/3 of the sequence - + 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(), leftQuery); Sequence* queryRight = new Sequence(q->getName(), rightQuery); - map tempIndexesRight = database->findClosest(queryRight, num); - map tempIndexesLeft = database->findClosest(queryLeft, num); + vector tempIndexesRight = database->findClosestMegaBlast(queryRight, num+1); + vector tempIndexesLeft = database->findClosestMegaBlast(queryLeft, num+1); - //merge results - vector mergedResults; + //if ((tempIndexesRight.size() != (num+1)) || (tempIndexesLeft.size() != (num+1))) { 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 porcess sequence " + q->getName()); m->mothurOutEndLine(); return refResults; } - map::iterator it; - map::iterator it2; + vector smaller; + vector larger; - //add in right guys merging common finds - for (it = tempIndexesRight.begin(); it != tempIndexesRight.end(); it++) { - it2 = tempIndexesLeft.find(it->first); - - if (it2 == tempIndexesLeft.end()) { //result only present in right - rank temp(it->first, it->second); - mergedResults.push_back(temp); + 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]; + } - }else { //result present in both save best score - float bestscore; - if (it->second > it2->second) { bestscore = it->second; } - else { bestscore = it2->second; } - - rank temp(it->first, bestscore); - mergedResults.push_back(temp); - - tempIndexesLeft.erase(it2); + //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]; } } - //add in unique left guys - for (it = tempIndexesLeft.begin(); it != tempIndexesLeft.end(); it++) { - rank temp(it->first, it->second); - mergedResults.push_back(temp); + 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; + delete database; + + return refResults; + } + catch(exception& e) { + m->errorOut(e, "Maligner", "getBlastSeqs"); + exit(1); + } +} +//*************************************************************************************************************** +vector Maligner::getKmerSeqs(Sequence* q, int num) { + try { + indexes.clear(); - sort(mergedResults.begin(), mergedResults.end(), compareMembers); + //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(), leftQuery); + Sequence* queryRight = new Sequence(q->getName(), rightQuery); + + vector tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted); + vector tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted); + + //merge results + map seen; + map::iterator it; + + vector mergedResults; + for (int i = 0; i < tempIndexesLeft.size(); i++) { + //add left if you havent already + it = seen.find(tempIndexesLeft[i]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesLeft[i]); + seen[tempIndexesLeft[i]] = tempIndexesLeft[i]; + } + + //add right if you havent already + it = seen.find(tempIndexesRight[i]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesRight[i]); + seen[tempIndexesRight[i]] = tempIndexesRight[i]; + } + } +//cout << q->getName() << endl; vector refResults; for (int i = 0; i < numWanted; i++) { - Sequence* temp = new Sequence(db[mergedResults[i].num]->getName(), db[mergedResults[i].num]->getAligned()); +//cout << db[mergedResults[i]]->getName() << endl; + Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); refResults.push_back(temp); -cout << db[mergedResults[i].num]->getName() << endl; + 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); } } - //***************************************************************************************************************