5 * Created by westcott on 9/23/09.
6 * Copyright 2009 Schloss Lab. All rights reserved.
12 #include "blastdb.hpp"
14 /***********************************************************************/
15 Maligner::Maligner(vector<Sequence*> temp, int num, int match, int misMatch, float div, int ms, int minCov, string mode, Database* dataLeft, Database* dataRight) :
16 db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov), searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight) { m = MothurOut::getInstance(); }
17 /***********************************************************************/
18 string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
21 outputResults.clear();
23 //make copy so trimming doesn't destroy query from calling class - remember to deallocate
24 query = new Sequence(q->getName(), q->getAligned());
28 if (searchMethod == "distance") {
29 //find closest seqs to query in template - returns copies of seqs so trim does not destroy - remember to deallocate
30 refSeqs = decalc->findClosest(query, db, numWanted, indexes);
31 }else if (searchMethod == "blast") {
32 refSeqs = getBlastSeqs(query, numWanted); //fills indexes
33 }else if (searchMethod == "kmer") {
34 refSeqs = getKmerSeqs(query, numWanted); //fills indexes
35 }else { m->mothurOut("not valid search."); exit(1); } //should never get here
37 if (m->control_pressed) { return chimera; }
39 refSeqs = minCoverageFilter(refSeqs);
41 if (refSeqs.size() < 2) {
42 for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
43 percentIdenticalQueryChimera = 0.0;
47 int chimeraPenalty = computeChimeraPenalty();
50 chimera = chimeraMaligner(chimeraPenalty, decalc);
52 if (m->control_pressed) { return chimera; }
57 for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
62 m->errorOut(e, "Maligner", "getResults");
66 /***********************************************************************/
67 string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
72 //trims seqs to first non gap char in all seqs and last non gap char in all seqs
73 spotMap = decalc->trimSeqs(query, refSeqs);
75 vector<Sequence*> temp = refSeqs;
76 temp.push_back(query);
79 //for (int i = 0; i < temp.size(); i++) { cout << temp[i]->getName() << '\n' << temp[i]->getAligned() << endl; } return "no";
81 vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
83 if (m->control_pressed) { return chimera; }
85 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
87 vector<score_struct> path = extractHighestPath(matrix);
89 if (m->control_pressed) { return chimera; }
91 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
93 if (trace.size() > 1) { chimera = "yes"; }
94 else { chimera = "no"; }
96 int traceStart = path[0].col;
97 int traceEnd = path[path.size()-1].col;
98 string queryInRange = query->getAligned();
99 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
101 string chimeraSeq = constructChimericSeq(trace, refSeqs);
103 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
105 if (m->control_pressed) { return chimera; }
107 //save output results
108 for (int i = 0; i < trace.size(); i++) {
109 int regionStart = trace[i].col;
110 int regionEnd = trace[i].oldCol;
111 int seqIndex = trace[i].row;
115 temp.parent = refSeqs[seqIndex]->getName();
116 temp.parentAligned = db[indexes[seqIndex]]->getAligned();
117 temp.nastRegionStart = spotMap[regionStart];
118 temp.nastRegionEnd = spotMap[regionEnd];
119 temp.regionStart = regionStart;
120 temp.regionEnd = regionEnd;
122 string parentInRange = refSeqs[seqIndex]->getAligned();
123 parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
125 temp.queryToParent = computePercentID(queryInRange, parentInRange);
126 temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
128 string queryInRegion = query->getAligned();
129 queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
131 string parentInRegion = refSeqs[seqIndex]->getAligned();
132 parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
134 temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
136 outputResults.push_back(temp);
141 catch(exception& e) {
142 m->errorOut(e, "Maligner", "chimeraMaligner");
146 /***********************************************************************/
147 //removes top matches that do not have minimum coverage with query.
148 vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
150 vector<Sequence*> newRefs;
152 string queryAligned = query->getAligned();
154 for (int i = 0; i < ref.size(); i++) {
156 string refAligned = ref[i]->getAligned();
162 for (int j = 0; j < queryAligned.length(); j++) {
164 if (isalpha(queryAligned[j])) {
167 if (isalpha(refAligned[j])) {
173 int coverage = ((numCovered/(float)numBases)*100);
175 //if coverage above minimum
176 if (coverage > minCoverage) {
177 newRefs.push_back(ref[i]);
183 catch(exception& e) {
184 m->errorOut(e, "Maligner", "minCoverageFilter");
188 /***********************************************************************/
189 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
190 int Maligner::computeChimeraPenalty() {
193 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
195 int penalty = int(numAllowable + 1) * misMatchPenalty;
200 catch(exception& e) {
201 m->errorOut(e, "Maligner", "computeChimeraPenalty");
205 /***********************************************************************/
206 //this is a vertical filter
207 void Maligner::verticalFilter(vector<Sequence*> seqs) {
209 vector<int> gaps; gaps.resize(query->getAligned().length(), 0);
211 string filterString = (string(query->getAligned().length(), '1'));
214 for (int i = 0; i < seqs.size(); i++) {
216 string seqAligned = seqs[i]->getAligned();
218 for (int j = 0; j < seqAligned.length(); j++) {
219 //if this spot is a gap
220 if ((seqAligned[j] == '-') || (seqAligned[j] == '.')) { gaps[j]++; }
224 //zero out spot where all sequences have blanks
225 int numColRemoved = 0;
226 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
227 if(gaps[i] == seqs.size()) { filterString[i] = '0'; numColRemoved++; }
230 map<int, int> newMap;
232 for (int i = 0; i < seqs.size(); i++) {
234 string seqAligned = seqs[i]->getAligned();
235 string newAligned = "";
238 for (int j = 0; j < seqAligned.length(); j++) {
239 //if this spot is not a gap
240 if (filterString[j] == '1') {
241 newAligned += seqAligned[j];
242 newMap[count] = spotMap[j];
247 seqs[i]->setAligned(newAligned);
252 catch(exception& e) {
253 m->errorOut(e, "Maligner", "verticalFilter");
257 //***************************************************************************************************************
258 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
261 vector< vector<score_struct> > m; m.resize(rows);
263 for (int i = 0; i < m.size(); i++) {
264 for (int j = 0; j < cols; j++) {
266 //initialize each cell
269 temp.score = -9999999;
273 m[i].push_back(temp);
279 catch(exception& e) {
280 m->errorOut(e, "Maligner", "buildScoreMatrix");
284 //***************************************************************************************************************
285 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
288 //get matrix dimensions
289 int numCols = query->getAligned().length();
290 int numRows = seqs.size();
292 //initialize first col
293 string queryAligned = query->getAligned();
294 for (int i = 0; i < numRows; i++) {
295 string subjectAligned = seqs[i]->getAligned();
298 if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
300 }else if (queryAligned[0] == subjectAligned[0]) {
301 ms[i][0].score = matchScore;
307 //fill rest of matrix
308 for (int j = 1; j < numCols; j++) { //iterate through matrix columns
310 for (int i = 0; i < numRows; i++) { //iterate through matrix rows
312 string subjectAligned = seqs[i]->getAligned();
314 int matchMisMatchScore = 0;
316 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
318 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
320 }else if (queryAligned[j] == subjectAligned[j]) {
321 matchMisMatchScore = matchScore;
322 }else if (queryAligned[j] != subjectAligned[j]) {
323 matchMisMatchScore = misMatchPenalty;
326 //compute score based on previous columns scores
327 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
329 int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
331 //you are not at yourself
332 if (prevIndex != i) { sumScore += penalty; }
333 if (sumScore < 0) { sumScore = 0; }
335 if (sumScore > ms[i][j].score) {
336 ms[i][j].score = sumScore;
337 ms[i][j].prev = prevIndex;
344 catch(exception& e) {
345 m->errorOut(e, "Maligner", "fillScoreMatrix");
349 //***************************************************************************************************************
350 vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
353 //get matrix dimensions
354 int numCols = query->getAligned().length();
355 int numRows = ms.size();
358 //find highest score scoring matrix
359 score_struct highestStruct;
360 int highestScore = 0;
362 for (int i = 0; i < numRows; i++) {
363 for (int j = 0; j < numCols; j++) {
364 if (ms[i][j].score > highestScore) {
365 highestScore = ms[i][j].score;
366 highestStruct = ms[i][j];
371 vector<score_struct> path;
373 int rowIndex = highestStruct.row;
374 int pos = highestStruct.col;
375 int score = highestStruct.score;
377 while (pos >= 0 && score > 0) {
378 score_struct temp = ms[rowIndex][pos];
381 if (score > 0) { path.push_back(temp); }
383 rowIndex = temp.prev;
387 reverse(path.begin(), path.end());
392 catch(exception& e) {
393 m->errorOut(e, "Maligner", "extractHighestPath");
397 //***************************************************************************************************************
398 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
400 vector<trace_struct> trace;
402 int region_index = path[0].row;
403 int region_start = path[0].col;
405 for (int i = 1; i < path.size(); i++) {
407 int next_region_index = path[i].row;
409 if (next_region_index != region_index) {
412 int col_index = path[i].col;
414 temp.col = region_start;
415 temp.oldCol = col_index-1;
416 temp.row = region_index;
418 trace.push_back(temp);
420 region_index = path[i].row;
421 region_start = col_index;
427 temp.col = region_start;
428 temp.oldCol = path[path.size()-1].col;
429 temp.row = region_index;
430 trace.push_back(temp);
435 catch(exception& e) {
436 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
440 //***************************************************************************************************************
441 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
445 for (int i = 0; i < trace.size(); i++) {
446 string seqAlign = seqs[trace[i].row]->getAligned();
447 seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
453 catch(exception& e) {
454 m->errorOut(e, "Maligner", "constructChimericSeq");
458 //***************************************************************************************************************
459 float Maligner::computePercentID(string queryAlign, string chimera) {
462 if (queryAlign.length() != chimera.length()) {
463 m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
464 m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine();
465 m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
471 int numIdentical = 0;
473 for (int i = 0; i < queryAlign.length(); i++) {
474 if ((isalpha(queryAlign[i])) || (isalpha(chimera[i]))) {
476 if (queryAlign[i] == chimera[i]) {
482 if (numBases == 0) { return 0; }
484 float percentIdentical = (numIdentical/(float)numBases) * 100;
486 return percentIdentical;
489 catch(exception& e) {
490 m->errorOut(e, "Maligner", "computePercentID");
494 //***************************************************************************************************************
495 vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
498 vector<Sequence*> refResults;
501 string queryUnAligned = q->getUnaligned();
502 string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
503 string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
505 Sequence* queryLeft = new Sequence(q->getName()+"left", leftQuery);
506 Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery);
508 vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
509 vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
511 //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; }
516 if (tempIndexesRight.size() < tempIndexesLeft.size()) { smaller = tempIndexesRight; larger = tempIndexesLeft; }
517 else { smaller = tempIndexesLeft; larger = tempIndexesRight; }
521 map<int, int>::iterator it;
523 vector<int> mergedResults;
524 for (int i = 0; i < smaller.size(); i++) {
525 //add left if you havent already
526 it = seen.find(smaller[i]);
527 if (it == seen.end()) {
528 mergedResults.push_back(smaller[i]);
529 seen[smaller[i]] = smaller[i];
532 //add right if you havent already
533 it = seen.find(larger[i]);
534 if (it == seen.end()) {
535 mergedResults.push_back(larger[i]);
536 seen[larger[i]] = larger[i];
540 for (int i = smaller.size(); i < larger.size(); i++) {
541 it = seen.find(larger[i]);
542 if (it == seen.end()) {
543 mergedResults.push_back(larger[i]);
544 seen[larger[i]] = larger[i];
548 if (mergedResults.size() < numWanted) { numWanted = mergedResults.size(); }
549 //cout << q->getName() << endl;
550 for (int i = 0; i < numWanted; i++) {
551 //cout << db[mergedResults[i]]->getName() << endl;
552 if (db[mergedResults[i]]->getName() != q->getName()) {
553 Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
554 refResults.push_back(temp);
555 indexes.push_back(mergedResults[i]);
557 //cout << mergedResults[i] << endl;
565 catch(exception& e) {
566 m->errorOut(e, "Maligner", "getBlastSeqs");
570 //***************************************************************************************************************
571 vector<Sequence*> Maligner::getKmerSeqs(Sequence* q, int num) {
576 string queryUnAligned = q->getUnaligned();
577 string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
578 string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
580 Sequence* queryLeft = new Sequence(q->getName(), leftQuery);
581 Sequence* queryRight = new Sequence(q->getName(), rightQuery);
583 vector<int> tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted);
584 vector<int> tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted);
588 map<int, int>::iterator it;
590 vector<int> mergedResults;
591 for (int i = 0; i < tempIndexesLeft.size(); i++) {
592 //add left if you havent already
593 it = seen.find(tempIndexesLeft[i]);
594 if (it == seen.end()) {
595 mergedResults.push_back(tempIndexesLeft[i]);
596 seen[tempIndexesLeft[i]] = tempIndexesLeft[i];
599 //add right if you havent already
600 it = seen.find(tempIndexesRight[i]);
601 if (it == seen.end()) {
602 mergedResults.push_back(tempIndexesRight[i]);
603 seen[tempIndexesRight[i]] = tempIndexesRight[i];
607 //cout << q->getName() << endl;
608 vector<Sequence*> refResults;
609 for (int i = 0; i < numWanted; i++) {
610 //cout << db[mergedResults[i]]->getName() << endl;
611 Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
612 refResults.push_back(temp);
613 indexes.push_back(mergedResults[i]);
621 catch(exception& e) {
622 m->errorOut(e, "Maligner", "getBlastSeqs");
626 //***************************************************************************************************************