X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=maligner.cpp;h=5da4b9f9cece7a9e64dd4a5dd621477fc365e94e;hb=8699fd2c88f10abda5fe32c89be061a89d673fd6;hp=c7f9be02054510ef2d7a995139fab5d273f91809;hpb=58cc09a375d1e1afceef3b036574ff21394ccc4d;p=mothur.git diff --git a/maligner.cpp b/maligner.cpp index c7f9be0..5da4b9f 100644 --- a/maligner.cpp +++ b/maligner.cpp @@ -8,26 +8,36 @@ */ #include "maligner.h" +#include "kmerdb.hpp" +#include "blastdb.hpp" /***********************************************************************/ -Maligner::Maligner(vector temp, int num, int match, int misMatch, float div, int minCov) : - db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minCoverage(minCov) {} +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) { +string Maligner::getResults(Sequence* q, DeCalculator* decalc) { try { + outputResults.clear(); + //make copy so trimming doesn't destroy query from calling class - remember to deallocate query = new Sequence(q->getName(), q->getAligned()); string chimera; - - decalc = new DeCalculator(); - //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); + 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, 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; @@ -35,37 +45,64 @@ string Maligner::getResults(Sequence* q) { } int chimeraPenalty = computeChimeraPenalty(); - - //trims seqs to first non gap char in all seqs and last non gap char in all seqs - decalc->trimSeqs(query, refSeqs); + + //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) { + m->errorOut(e, "Maligner", "getResults"); + exit(1); + } +} +/***********************************************************************/ +string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { + try { + + string chimera; + //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); - delete decalc; + if (m->control_pressed) { return chimera; } //save output results for (int i = 0; i < trace.size(); i++) { @@ -76,6 +113,9 @@ string Maligner::getResults(Sequence* q) { results temp; temp.parent = refSeqs[seqIndex]->getName(); + temp.parentAligned = db[indexes[seqIndex]]->getAligned(); + temp.nastRegionStart = spotMap[regionStart]; + temp.nastRegionEnd = spotMap[regionEnd]; temp.regionStart = regionStart; temp.regionEnd = regionEnd; @@ -95,15 +135,11 @@ string Maligner::getResults(Sequence* q) { outputResults.push_back(temp); } - - //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", "chimeraMaligner"); exit(1); } } @@ -145,7 +181,7 @@ vector Maligner::minCoverageFilter(vector ref){ return newRefs; } catch(exception& e) { - errorOut(e, "Maligner", "minCoverageFilter"); + m->errorOut(e, "Maligner", "minCoverageFilter"); exit(1); } } @@ -162,7 +198,7 @@ int Maligner::computeChimeraPenalty() { } catch(exception& e) { - errorOut(e, "Maligner", "computeChimeraPenalty"); + m->errorOut(e, "Maligner", "computeChimeraPenalty"); exit(1); } } @@ -191,24 +227,30 @@ void Maligner::verticalFilter(vector seqs) { if(gaps[i] == seqs.size()) { filterString[i] = '0'; numColRemoved++; } } + map newMap; //for each sequence for (int i = 0; i < seqs.size(); i++) { string seqAligned = seqs[i]->getAligned(); string newAligned = ""; + int count = 0; for (int j = 0; j < seqAligned.length(); j++) { //if this spot is not a gap - if (filterString[j] == '1') { newAligned += seqAligned[j]; } + if (filterString[j] == '1') { + newAligned += seqAligned[j]; + newMap[count] = spotMap[j]; + count++; + } } seqs[i]->setAligned(newAligned); } - + spotMap = newMap; } catch(exception& e) { - errorOut(e, "Maligner", "verticalFilter"); + m->errorOut(e, "Maligner", "verticalFilter"); exit(1); } } @@ -235,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 @@ -254,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; } } } @@ -300,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 @@ -319,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]; } } } @@ -333,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); } @@ -348,7 +390,7 @@ vector Maligner::extractHighestPath(vector > } catch(exception& e) { - errorOut(e, "Maligner", "extractHighestPath"); + m->errorOut(e, "Maligner", "extractHighestPath"); exit(1); } } @@ -391,7 +433,7 @@ vector Maligner::mapTraceRegionsToAlignment(vector p } catch(exception& e) { - errorOut(e, "Maligner", "mapTraceRegionsToAlignment"); + m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment"); exit(1); } } @@ -409,7 +451,7 @@ string Maligner::constructChimericSeq(vector trace, vectorerrorOut(e, "Maligner", "constructChimericSeq"); exit(1); } } @@ -418,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; } @@ -445,7 +487,139 @@ 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 { + 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(), leftQuery); + Sequence* queryRight = new Sequence(q->getName(), rightQuery); + + vector tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1); + vector tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1); + + //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; } + + 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(); + 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++) { +//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; + + return refResults; + } + catch(exception& e) { + m->errorOut(e, "Maligner", "getBlastSeqs"); exit(1); } }