-vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
- try {
- indexes.clear();
- vector<Sequence*> 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]); }
- database->generateDB();
- database->setNumSeqs(db.size());
-
- //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<int> tempIndexesRight = database->findClosestMegaBlast(queryRight, num+1);
- vector<int> tempIndexesLeft = database->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<int> smaller;
- vector<int> larger;
-
- if (tempIndexesRight.size() < tempIndexesLeft.size()) { smaller = tempIndexesRight; larger = tempIndexesLeft; }
- else { smaller = tempIndexesLeft; larger = tempIndexesRight; }
-
- //merge results
- map<int, int> seen;
- map<int, int>::iterator it;
-
- 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];
- }
-
- //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;
- delete database;
-
- return refResults;
- }
- catch(exception& e) {
- m->errorOut(e, "Maligner", "getBlastSeqs");
- exit(1);
- }
-}
-//***************************************************************************************************************
-vector<Sequence*> 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<int> tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted);
- vector<int> tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted);
-
- //merge results
- map<int, int> seen;
- map<int, int>::iterator it;
-
- vector<int> 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<Sequence*> 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);
- }
-}
-//***************************************************************************************************************
-