*/
#include "maligner.h"
-#include "database.hpp"
+#include "kmerdb.hpp"
#include "blastdb.hpp"
/***********************************************************************/
-Maligner::Maligner(vector<Sequence*> 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<Sequence*> 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();
+
+// cout << matchScore << '\t' << misMatchPenalty << endl;
+//
+// matchScore = 1;
+// misMatchPenalty = -1;
+
+ }
+
/***********************************************************************/
string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
try {
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
+
+ if (m->control_pressed) { return chimera; }
-//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;
}
int chimeraPenalty = computeChimeraPenalty();
-//cout << chimeraPenalty << endl;
+
//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);
}
}
//trims seqs to first non gap char in all seqs and last non gap char in all seqs
spotMap = decalc->trimSeqs(query, refSeqs);
+ //you trimmed the whole sequence, skip
+ if (query->getAligned() == "") { return "no"; }
+
vector<Sequence*> temp = refSeqs;
temp.push_back(query);
-
+
+
verticalFilter(temp);
vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
- fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
-
- vector<score_struct> path = extractHighestPath(matrix);
+ if (m->control_pressed) { return chimera; }
- vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
+ fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
+
+ vector<trace_struct> trace = extractHighestPath(matrix);
+
+// cout << "traces\n";
+// for(int i=0;i<trace.size();i++){
+// cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
+// }
if (trace.size() > 1) { chimera = "yes"; }
- else { chimera = "no"; }
-
- int traceStart = path[0].col;
- int traceEnd = path[path.size()-1].col;
+ else { chimera = "no"; return chimera; }
+ int traceStart = trace[0].col;
+ int traceEnd = trace[trace.size()-1].oldCol;
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++) {
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.queryToParent = computePercentID(queryInRange, parentInRange);
temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
-
+
string queryInRegion = query->getAligned();
queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
outputResults.push_back(temp);
}
-
+
return chimera;
}
catch(exception& e) {
- errorOut(e, "Maligner", "chimeraMaligner");
+ m->errorOut(e, "Maligner", "chimeraMaligner");
exit(1);
}
}
//if coverage above minimum
if (coverage > minCoverage) {
newRefs.push_back(ref[i]);
+ }else {
+ delete ref[i];
}
}
return newRefs;
}
catch(exception& e) {
- errorOut(e, "Maligner", "minCoverageFilter");
+ m->errorOut(e, "Maligner", "minCoverageFilter");
exit(1);
}
}
try {
int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
-
+
+// if(numAllowable < 1){ numAllowable = 1; }
+
int penalty = int(numAllowable + 1) * misMatchPenalty;
-
+
return penalty;
}
catch(exception& e) {
- errorOut(e, "Maligner", "computeChimeraPenalty");
+ m->errorOut(e, "Maligner", "computeChimeraPenalty");
exit(1);
}
}
spotMap = newMap;
}
catch(exception& e) {
- errorOut(e, "Maligner", "verticalFilter");
+ m->errorOut(e, "Maligner", "verticalFilter");
exit(1);
}
}
vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
try{
- vector< vector<score_struct> > m; m.resize(rows);
+ vector< vector<score_struct> > m(rows);
- for (int i = 0; i < m.size(); i++) {
+ for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
//initialize each cell
return m;
}
catch(exception& e) {
- errorOut(e, "Maligner", "buildScoreMatrix");
+ m->errorOut(e, "Maligner", "buildScoreMatrix");
exit(1);
}
}
+
//***************************************************************************************************************
-void Maligner::fillScoreMatrix(vector<vector<score_struct> >& m, vector<Sequence*> seqs, int penalty) {
+
+void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
try{
//get matrix dimensions
//are you both gaps?
if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
- m[i][0].score = 0;
- }else if (queryAligned[0] == subjectAligned[0]) {
- m[i][0].score = matchScore;
+ ms[i][0].score = 0;
+// ms[i][0].mismatches = 0;
+ }else if (queryAligned[0] == subjectAligned[0] || subjectAligned[0] == 'N') {
+ ms[i][0].score = matchScore;
+// ms[i][0].mismatches = 0;
}else{
- m[i][0].score = 0;
+ ms[i][0].score = 0;
+// ms[i][0].mismatches = 1;
}
}
if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
//leave the same
}else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
+ matchMisMatchScore = matchScore;
//leave the same
}else if (queryAligned[j] == subjectAligned[j]) {
matchMisMatchScore = matchScore;
+// ms[i][j].mismatches = ms[i][j-1].mismatches;
}else if (queryAligned[j] != subjectAligned[j]) {
matchMisMatchScore = misMatchPenalty;
+// ms[i][j].mismatches = ms[i][j-1].mismatches + 1;
}
//compute score based on previous columns scores
for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
- int sumScore = matchMisMatchScore + m[prevIndex][j-1].score;
+ int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
//you are not at yourself
if (prevIndex != i) { sumScore += penalty; }
if (sumScore < 0) { sumScore = 0; }
- if (sumScore > 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;
}
}
}
}
+// for(int i=0;i<numRows;i++){
+// cout << seqs[i]->getName();
+// for(int j=0;j<numCols;j++){
+// cout << '\t' << ms[i][j].mismatches;
+// }
+// cout << endl;
+// }
+// cout << endl;
+// for(int i=0;i<numRows;i++){
+// cout << seqs[i]->getName();
+// for(int j=0;j<numCols;j++){
+// cout << '\t' << ms[i][j].score;
+// }
+// cout << endl;
+// }
+// cout << endl;
+// for(int i=0;i<numRows;i++){
+// cout << seqs[i]->getName();
+// for(int j=0;j<numCols;j++){
+// cout << '\t' << ms[i][j].prev;
+// }
+// cout << endl;
+// }
+
+
}
catch(exception& e) {
- errorOut(e, "Maligner", "fillScoreMatrix");
+ m->errorOut(e, "Maligner", "fillScoreMatrix");
exit(1);
}
}
+
//***************************************************************************************************************
-vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > m) {
+
+vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
try {
+
//get matrix dimensions
int numCols = query->getAligned().length();
- int numRows = m.size();
+ int numRows = ms.size();
//find highest score scoring matrix
- score_struct highestStruct;
+ vector<score_struct> highestStruct;
int highestScore = 0;
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.resize(0);
+ highestStruct.push_back(ms[i][j]);
+ }
+ else if(ms[i][j].score == highestScore){
+ highestStruct.push_back(ms[i][j]);
}
}
}
-
- vector<score_struct> path;
+
+// cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl;
- int rowIndex = highestStruct.row;
- int pos = highestStruct.col;
- int score = highestStruct.score;
+ vector<trace_struct> maxTrace;
+ double maxPercentIdenticalQueryAntiChimera = 0;
- while (pos >= 0 && score > 0) {
- score_struct temp = m[rowIndex][pos];
- score = temp.score;
+ for(int i=0;i<highestStruct.size();i++){
+
+ vector<score_struct> path;
+
+ int rowIndex = highestStruct[i].row;
+ int pos = highestStruct[i].col;
+ int score = highestStruct[i].score;
+
+ while (pos >= 0 && score > 0) {
+ score_struct temp = ms[rowIndex][pos];
+ score = temp.score;
+
+ if (score > 0) { path.push_back(temp); }
+
+ rowIndex = temp.prev;
+ pos--;
+ }
+
+ reverse(path.begin(), path.end());
+
+ vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
+
+// cout << "traces\n";
+// for(int j=0;j<trace.size();j++){
+// cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
+// }
+
+ int traceStart = path[0].col;
+ int traceEnd = path[path.size()-1].col;
+// cout << "traceStart/End\t" << traceStart << '\t' << traceEnd << endl;
- if (score > 0) { path.push_back(temp); }
+ string queryInRange = query->getAligned();
+ queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
+// cout << "here" << endl;
+ string chimeraSeq = constructChimericSeq(trace, refSeqs);
+ string antiChimeraSeq = constructAntiChimericSeq(trace, refSeqs);
+
+ percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
+ double percentIdenticalQueryAntiChimera = computePercentID(queryInRange, antiChimeraSeq);
+// cout << i << '\t' << percentIdenticalQueryChimera << '\t' << percentIdenticalQueryAntiChimera << endl;
- rowIndex = temp.prev;
- pos--;
+ if(percentIdenticalQueryAntiChimera > maxPercentIdenticalQueryAntiChimera){
+ maxPercentIdenticalQueryAntiChimera = percentIdenticalQueryAntiChimera;
+ maxTrace = trace;
+ }
}
-
- reverse(path.begin(), path.end());
-
- return path;
+// cout << maxPercentIdenticalQueryAntiChimera << endl;
+ return maxTrace;
}
catch(exception& e) {
- errorOut(e, "Maligner", "extractHighestPath");
+ m->errorOut(e, "Maligner", "extractHighestPath");
exit(1);
}
}
+
//***************************************************************************************************************
+
vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
try {
vector<trace_struct> trace;
temp.row = region_index;
trace.push_back(temp);
+// cout << endl;
+// cout << trace.size() << endl;
+// for(int i=0;i<trace.size();i++){
+// cout << seqs[trace[i].row]->getName() << endl;
+// }
+// cout << endl;
+
return trace;
}
catch(exception& e) {
- errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
+ m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
exit(1);
}
}
+
//***************************************************************************************************************
+
string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
try {
string chimera = "";
for (int i = 0; i < trace.size(); i++) {
+// cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
+
string seqAlign = seqs[trace[i].row]->getAligned();
seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
chimera += seqAlign;
return chimera;
}
catch(exception& e) {
- errorOut(e, "Maligner", "constructChimericSeq");
+ m->errorOut(e, "Maligner", "constructChimericSeq");
+ exit(1);
+ }
+}
+
+//***************************************************************************************************************
+
+string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
+ try {
+ string antiChimera = "";
+
+ for (int i = 0; i < trace.size(); i++) {
+// cout << i << '\t' << (trace.size() - i - 1) << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
+
+ int oppositeIndex = trace.size() - i - 1;
+
+ string seqAlign = seqs[trace[oppositeIndex].row]->getAligned();
+ seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
+ antiChimera += seqAlign;
+ }
+
+ return antiChimera;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "Maligner", "constructChimericSeq");
exit(1);
}
}
+
//***************************************************************************************************************
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;
}
}
catch(exception& e) {
- errorOut(e, "Maligner", "computePercentID");
+ m->errorOut(e, "Maligner", "computePercentID");
exit(1);
}
}
//***************************************************************************************************************
vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
- try {
-cout << q->getName() << endl;
- //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());
-
+ try {
+ indexes.clear();
+ vector<Sequence*> refResults;
+
//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
-
- Sequence* queryLeft = new Sequence(q->getName(), leftQuery);
- Sequence* queryRight = new Sequence(q->getName(), rightQuery);
+ 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
- map<int, float> tempIndexesRight = database->findClosest(queryRight, num);
- map<int, float> tempIndexesLeft = database->findClosest(queryLeft, num);
+ 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; }
+ }
- //merge results
- vector<rank> mergedResults;
+ //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];
+ }
+
+ lasti = i;
+ if (mergedResults.size() > num) { break; }
+ }
+ }
- map<int, float>::iterator it;
- map<int, float>::iterator it2;
- //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);
-
- }else { //result present in both save best score
- float bestscore;
- if (it->second > it2->second) { bestscore = it->second; }
- else { bestscore = it2->second; }
+ //add in any ties from smaller
+ while (iSmaller < smaller.size()) {
+ if (smallerScores[iSmaller] == lastSmaller) {
+ it = seen.find(smaller[iSmaller]);
- rank temp(it->first, bestscore);
- mergedResults.push_back(temp);
+ if (it == seen.end()) {
+ mergedResults.push_back(smaller[iSmaller]);
+ seen[smaller[iSmaller]] = smaller[iSmaller];
+ }
+ }
+ else { break; }
+ iSmaller++;
+ }
+
+ lasti++;
+ //add in any ties from larger
+ while (lasti < larger.size()) {
+ if (largerScores[lasti] == lastLarger) { //is it a tie
+ it = seen.find(larger[lasti]);
- tempIndexesLeft.erase(it2);
+ if (it == seen.end()) { //we haven't already seen it
+ mergedResults.push_back(larger[lasti]);
+ seen[larger[lasti]] = larger[lasti];
+ }
}
+ else { break; }
+ lasti++;
}
- //add in unique left guys
- for (it = tempIndexesLeft.begin(); it != tempIndexesLeft.end(); it++) {
- rank temp(it->first, it->second);
- mergedResults.push_back(temp);
+ numWanted = mergedResults.size();
+
+ if (mergedResults.size() < numWanted) { numWanted = mergedResults.size(); }
+//cout << q->getName() << " merged results size = " << mergedResults.size() << '\t' << "numwanted = " << numWanted << endl;
+ for (int i = 0; i < numWanted; i++) {
+//cout << db[mergedResults[i]]->getName() << '\t' << mergedResults[i] << 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 << "done " << q->getName() << endl;
+ delete queryRight;
+ delete queryLeft;
+
+ 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);
- sort(mergedResults.begin(), mergedResults.end(), compareMembers);
+ 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++) {
- 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", "getKmerSeqs");
exit(1);
}
}
-
//***************************************************************************************************************
projects / mothur.git / blobdiff