*/
#include "maligner.h"
-#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, Database* dataLeft, Database* dataRight) :
- db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov), searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight) {
-
+/***********************************************************************/ //int num, int match, int misMatch, , string mode, Database* dataLeft, Database* dataRight
+Maligner::Maligner(vector<Sequence*> temp, int match, int misMatch, float div, int ms, int minCov) : db(temp), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov) {
+ //numWanted(num), , 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 == "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
-
- if (m->control_pressed) { return chimera; }
+ //copy refSeqs so that filter does not effect original
+ for(int i = 0; i < db.size(); i++) {
+ Sequence* newSeq = new Sequence(db[i]->getName(), db[i]->getAligned());
+ refSeqs.push_back(newSeq);
+ }
refSeqs = minCoverageFilter(refSeqs);
-
-
if (refSeqs.size() < 2) {
for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
percentIdenticalQueryChimera = 0.0;
}
int chimeraPenalty = computeChimeraPenalty();
- //cout << "chimeraPenalty = " << chimeraPenalty << endl;
+
//fills outputResults
chimera = chimeraMaligner(chimeraPenalty, decalc);
if (m->control_pressed) { return chimera; }
fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
-
+
+ vector<score_struct> path = extractHighestPath(matrix);
+
+ if (m->control_pressed) { return chimera; }
+
+ vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
+
+ if (trace.size() > 1) { chimera = "yes"; }
+ else { chimera = "no"; return chimera; }
+
+ int traceStart = path[0].col;
+ int traceEnd = path[path.size()-1].col;
+ string queryInRange = query->getAligned();
+ queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart));
+
+ string chimeraSeq = constructChimericSeq(trace, refSeqs);
+
+ percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
+
+ /*
vector<trace_struct> trace = extractHighestPath(matrix);
//cout << "traces\n";
int traceStart = trace[0].col;
int traceEnd = trace[trace.size()-1].oldCol;
string queryInRange = query->getAligned();
- queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
+ queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart));*/
if (m->control_pressed) { return chimera; }
results temp;
temp.parent = refSeqs[seqIndex]->getName();
- temp.parentAligned = db[indexes[seqIndex]]->getAligned();
+ temp.parentAligned = db[seqIndex]->getAligned();
temp.nastRegionStart = spotMap[regionStart];
temp.nastRegionEnd = spotMap[regionEnd];
temp.regionStart = regionStart;
temp.regionEnd = regionEnd;
string parentInRange = refSeqs[seqIndex]->getAligned();
- parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
+ parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart));
temp.queryToParent = computePercentID(queryInRange, parentInRange);
temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
string queryInRegion = query->getAligned();
- queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
+ queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart));
string parentInRegion = refSeqs[seqIndex]->getAligned();
- parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
+ parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart));
temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
-
+
+ cout << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << endl;
outputResults.push_back(temp);
}
exit(1);
}
}
-
//***************************************************************************************************************
+vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
+ try {
+
+ //get matrix dimensions
+ int numCols = query->getAligned().length();
+ int numRows = ms.size();
+
+
+ //find highest score scoring matrix
+ score_struct highestStruct;
+ int highestScore = 0;
+
+ for (int i = 0; i < numRows; i++) {
+ for (int j = 0; j < numCols; j++) {
+ if (ms[i][j].score > highestScore) {
+ highestScore = ms[i][j].score;
+ highestStruct = ms[i][j];
+ }
+ }
+ }
+
+ vector<score_struct> path;
+
+ int rowIndex = highestStruct.row;
+ int pos = highestStruct.col;
+ int score = highestStruct.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());
+
+ return path;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "Maligner", "extractHighestPath");
+ exit(1);
+ }
+}
+//***************************************************************************************************************
+vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
+ try {
+ vector<trace_struct> trace;
+
+ int region_index = path[0].row;
+ int region_start = path[0].col;
+
+ for (int i = 1; i < path.size(); i++) {
+
+ int next_region_index = path[i].row;
+
+ if (next_region_index != region_index) {
+
+ // add trace region
+ int col_index = path[i].col;
+ trace_struct temp;
+ temp.col = region_start;
+ temp.oldCol = col_index-1;
+ temp.row = region_index;
+
+ trace.push_back(temp);
+
+ region_index = path[i].row;
+ region_start = col_index;
+ }
+ }
+
+ // get last one
+ trace_struct temp;
+ temp.col = region_start;
+ temp.oldCol = path[path.size()-1].col;
+ temp.row = region_index;
+ trace.push_back(temp);
+
+ return trace;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
+ exit(1);
+ }
+}
+
+/***************************************************************************************************************
vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
try {
exit(1);
}
}
-
+*/
//***************************************************************************************************************
string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
chimera += seqAlign;
}
-
+
+ if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); }
return chimera;
}
catch(exception& e) {
if (queryAlign.length() != chimera.length()) {
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(queryAlign.length())); m->mothurOutEndLine();
m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
return -1.0;
}
-
- int numBases = 0;
int numIdentical = 0;
-
+ int countA = 0;
+ int countB = 0;
for (int i = 0; i < queryAlign.length(); i++) {
- if ((isalpha(queryAlign[i])) || (isalpha(chimera[i]))) {
- numBases++;
- if (queryAlign[i] == chimera[i]) {
- numIdentical++;
+ if (((queryAlign[i] != 'G') && (queryAlign[i] != 'T') && (queryAlign[i] != 'A') && (queryAlign[i] != 'C')&& (queryAlign[i] != '.') && (queryAlign[i] != '-')) ||
+ ((chimera[i] != 'G') && (chimera[i] != 'T') && (chimera[i] != 'A') && (chimera[i] != 'C')&& (chimera[i] != '.') && (chimera[i] != '-'))) {}
+ else {
+
+ bool charA = false; bool charB = false;
+ if ((queryAlign[i] == 'G') || (queryAlign[i] == 'T') || (queryAlign[i] == 'A') || (queryAlign[i] == 'C')) { charA = true; }
+ if ((chimera[i] == 'G') || (chimera[i] == 'T') || (chimera[i] == 'A') || (chimera[i] == 'C')) { charB = true; }
+
+ if (charA || charB) {
+
+ if (charA) { countA++; }
+ if (charB) { countB++; }
+
+ if (queryAlign[i] == chimera[i]) {
+ numIdentical++;
+ }
}
}
}
-
+
+ float numBases = (countA + countB) /(float) 2;
+
if (numBases == 0) { return 0; }
float percentIdentical = (numIdentical/(float)numBases) * 100;
}
}
//***************************************************************************************************************
-vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
- try {
- indexes.clear();
- vector<Sequence*> 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()+"left", leftQuery);
- Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery);
-
- vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
- vector<float> leftScores = databaseLeft->getSearchScores();
- vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
- vector<float> rightScores = databaseLeft->getSearchScores();
-
- //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];
- }
-
- lasti = i;
- if (mergedResults.size() > num) { break; }
- }
- }
-
-
- //add in any ties from smaller
- while (iSmaller < smaller.size()) {
- if (smallerScores[iSmaller] == lastSmaller) {
- it = seen.find(smaller[iSmaller]);
-
- 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]);
-
- if (it == seen.end()) { //we haven't already seen it
- mergedResults.push_back(larger[lasti]);
- seen[larger[lasti]] = larger[lasti];
- }
- }
- else { break; }
- lasti++;
- }
- */
-
- for (int i = smaller.size(); i < larger.size(); 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];
- }
- }
- 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 << db[mergedResults[i]]->getName() << endl;
- }
-
-//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);
-
- vector<int> tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted);
- vector<int> tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted);
- vector<float> scoresLeft = databaseLeft->getSearchScores();
- vector<float> scoresRight = databaseRight->getSearchScores();
-
- //merge results
- map<int, int> seen;
- map<int, int>::iterator it;
- float lastRight = scoresRight[0];
- float lastLeft = scoresLeft[0];
- //int lasti = 0;
- 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];
- lastLeft = scoresLeft[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];
- lastRight = scoresRight[i];
- }
-
- //if (mergedResults.size() > numWanted) { lasti = i; break; } //you have enough results
- }
-
- //add in sequences with same distance as last sequence added
- /*lasti++;
- int i = lasti;
- while (i < tempIndexesLeft.size()) {
- if (scoresLeft[i] == lastLeft) {
- it = seen.find(tempIndexesLeft[i]);
-
- if (it == seen.end()) {
- mergedResults.push_back(tempIndexesLeft[i]);
- seen[tempIndexesLeft[i]] = tempIndexesLeft[i];
- }
- }
- else { break; }
- i++;
- }
-
- // cout << "lastRight\t" << lastRight << endl;
- //add in sequences with same distance as last sequence added
- i = lasti;
- while (i < tempIndexesRight.size()) {
- if (scoresRight[i] == lastRight) {
- it = seen.find(tempIndexesRight[i]);
-
- if (it == seen.end()) {
- mergedResults.push_back(tempIndexesRight[i]);
- seen[tempIndexesRight[i]] = tempIndexesRight[i];
- }
- }
- else { break; }
- i++;
- }*/
-
- numWanted = mergedResults.size();
-
- if (numWanted > mergedResults.size()) {
- //m->mothurOut("numwanted is larger than the number of template sequences, adjusting numwanted."); m->mothurOutEndLine();
- numWanted = mergedResults.size();
- }
-
-
-//cout << q->getName() << endl;
- vector<Sequence*> refResults;
- 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 << endl;
- delete queryRight;
- delete queryLeft;
-
- return refResults;
- }
- catch(exception& e) {
- m->errorOut(e, "Maligner", "getKmerSeqs");
- exit(1);
- }
-}
-//***************************************************************************************************************
-
projects / mothur.git / blobdiff