- //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<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
- vector<float> leftScores = databaseLeft->getMegaBlastSearchScores();
- vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
- vector<float> rightScores = databaseLeft->getMegaBlastSearchScores();
-
- //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<int> smaller;
- vector<float> smallerScores;
- vector<int> larger;
- vector<float> largerScores;
-
- if (tempIndexesRight.size() < tempIndexesLeft.size()) { smaller = tempIndexesRight; smallerScores = rightScores; larger = tempIndexesLeft; largerScores = leftScores; }
- else { smaller = tempIndexesLeft; smallerScores = leftScores; larger = tempIndexesRight; largerScores = rightScores; }
-
- //for (int i = 0; i < smaller.size(); i++) { cout << "smaller = " << smaller[i] << '\t' << smallerScores[i] << endl; }
- //cout << endl;
- //for (int i = 0; i < larger.size(); i++) { cout << "larger = " << larger[i] << '\t' << largerScores[i] << endl; }
-
- //merge results
- map<int, int> seen;
- map<int, int>::iterator it;
- float lastSmaller = smallerScores[0];
- float lastLarger = largerScores[0];
- int lasti = 0;
- vector<int> 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];
- lastSmaller = smallerScores[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];
- lastLarger = largerScores[i];
- }
-
- lasti = i;
- if (mergedResults.size() > num) { break; }
- }
-
- //save lasti for smaller ties below
- lasti++;
- int iSmaller = lasti;
-
- if (!(mergedResults.size() > num)) { //if we still need more results.
- 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];
- lastLarger = largerScores[i];