5 * Created by westcott on 9/23/09.
6 * Copyright 2009 Schloss Lab. All rights reserved.
12 /***********************************************************************/ //int num, int match, int misMatch, , string mode, Database* dataLeft, Database* dataRight
13 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) {
14 //numWanted(num), , searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight)
16 m = MothurOut::getInstance();
19 /***********************************************************************/
20 string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
23 outputResults.clear();
25 //make copy so trimming doesn't destroy query from calling class - remember to deallocate
26 query = new Sequence(q->getName(), q->getAligned());
30 //copy refSeqs so that filter does not effect original
31 for(int i = 0; i < db.size(); i++) {
32 Sequence* newSeq = new Sequence(db[i]->getName(), db[i]->getAligned());
33 refSeqs.push_back(newSeq);
36 refSeqs = minCoverageFilter(refSeqs);
38 if (refSeqs.size() < 2) {
39 for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
40 percentIdenticalQueryChimera = 0.0;
44 int chimeraPenalty = computeChimeraPenalty();
47 chimera = chimeraMaligner(chimeraPenalty, decalc);
49 if (m->control_pressed) { return chimera; }
54 for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
59 m->errorOut(e, "Maligner", "getResults");
63 /***********************************************************************/
64 string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
69 //trims seqs to first non gap char in all seqs and last non gap char in all seqs
70 spotMap = decalc->trimSeqs(query, refSeqs);
72 //you trimmed the whole sequence, skip
73 if (query->getAligned() == "") { return "no"; }
75 vector<Sequence*> temp = refSeqs;
76 temp.push_back(query);
80 //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl ; }//<< refSeqs[i]->getAligned() << endl
82 vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
84 if (m->control_pressed) { return chimera; }
86 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
88 vector<score_struct> path = extractHighestPath(matrix);
90 if (m->control_pressed) { return chimera; }
92 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
94 if (trace.size() > 1) { chimera = "yes"; }
95 else { chimera = "no"; return chimera; }
97 int traceStart = path[0].col;
98 int traceEnd = path[path.size()-1].col;
99 string queryInRange = query->getAligned();
100 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart));
102 string chimeraSeq = constructChimericSeq(trace, refSeqs);
104 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
107 vector<trace_struct> trace = extractHighestPath(matrix);
109 //cout << "traces\n";
110 //for(int i=0;i<trace.size();i++){
111 // cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
114 if (trace.size() > 1) { chimera = "yes"; }
115 else { chimera = "no"; return chimera; }
117 int traceStart = trace[0].col;
118 int traceEnd = trace[trace.size()-1].oldCol;
119 string queryInRange = query->getAligned();
120 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart));*/
122 if (m->control_pressed) { return chimera; }
124 //save output results
125 for (int i = 0; i < trace.size(); i++) {
126 int regionStart = trace[i].col;
127 int regionEnd = trace[i].oldCol;
128 int seqIndex = trace[i].row;
132 temp.parent = refSeqs[seqIndex]->getName();
133 temp.parentAligned = db[seqIndex]->getAligned();
134 temp.nastRegionStart = spotMap[regionStart];
135 temp.nastRegionEnd = spotMap[regionEnd];
136 temp.regionStart = regionStart;
137 temp.regionEnd = regionEnd;
139 string parentInRange = refSeqs[seqIndex]->getAligned();
140 parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart));
142 temp.queryToParent = computePercentID(queryInRange, parentInRange);
143 temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
145 string queryInRegion = query->getAligned();
146 queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart));
148 string parentInRegion = refSeqs[seqIndex]->getAligned();
149 parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart));
151 temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
153 //cout << query->getName() << '\t' << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << ", " << temp.divR << endl;
155 outputResults.push_back(temp);
160 catch(exception& e) {
161 m->errorOut(e, "Maligner", "chimeraMaligner");
165 /***********************************************************************/
166 //removes top matches that do not have minimum coverage with query.
167 vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
169 vector<Sequence*> newRefs;
171 string queryAligned = query->getAligned();
173 for (int i = 0; i < ref.size(); i++) {
175 string refAligned = ref[i]->getAligned();
181 for (int j = 0; j < queryAligned.length(); j++) {
183 if (isalpha(queryAligned[j])) {
186 if (isalpha(refAligned[j])) {
192 int coverage = ((numCovered/(float)numBases)*100);
194 //if coverage above minimum
195 if (coverage > minCoverage) {
196 newRefs.push_back(ref[i]);
204 catch(exception& e) {
205 m->errorOut(e, "Maligner", "minCoverageFilter");
209 /***********************************************************************/
210 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
211 int Maligner::computeChimeraPenalty() {
214 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
216 // if(numAllowable < 1){ numAllowable = 1; }
218 int penalty = int(numAllowable + 1) * misMatchPenalty;
223 catch(exception& e) {
224 m->errorOut(e, "Maligner", "computeChimeraPenalty");
228 /***********************************************************************/
229 //this is a vertical filter
230 void Maligner::verticalFilter(vector<Sequence*> seqs) {
232 vector<int> gaps; gaps.resize(query->getAligned().length(), 0);
234 string filterString = (string(query->getAligned().length(), '1'));
237 for (int i = 0; i < seqs.size(); i++) {
239 string seqAligned = seqs[i]->getAligned();
241 for (int j = 0; j < seqAligned.length(); j++) {
242 //if this spot is a gap
243 if ((seqAligned[j] == '-') || (seqAligned[j] == '.')) { gaps[j]++; }
247 //zero out spot where all sequences have blanks
248 int numColRemoved = 0;
249 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
250 if(gaps[i] == seqs.size()) { filterString[i] = '0'; numColRemoved++; }
253 map<int, int> newMap;
255 for (int i = 0; i < seqs.size(); i++) {
257 string seqAligned = seqs[i]->getAligned();
258 string newAligned = "";
261 for (int j = 0; j < seqAligned.length(); j++) {
262 //if this spot is not a gap
263 if (filterString[j] == '1') {
264 newAligned += seqAligned[j];
265 newMap[count] = spotMap[j];
270 seqs[i]->setAligned(newAligned);
275 catch(exception& e) {
276 m->errorOut(e, "Maligner", "verticalFilter");
280 //***************************************************************************************************************
281 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
284 vector< vector<score_struct> > m(rows);
286 for (int i = 0; i < rows; i++) {
287 for (int j = 0; j < cols; j++) {
289 //initialize each cell
292 temp.score = -9999999;
296 m[i].push_back(temp);
302 catch(exception& e) {
303 m->errorOut(e, "Maligner", "buildScoreMatrix");
308 //***************************************************************************************************************
310 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
313 //get matrix dimensions
314 int numCols = query->getAligned().length();
315 int numRows = seqs.size();
317 //initialize first col
318 string queryAligned = query->getAligned();
319 for (int i = 0; i < numRows; i++) {
320 string subjectAligned = seqs[i]->getAligned();
323 if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
325 // ms[i][0].mismatches = 0;
326 }else if (queryAligned[0] == subjectAligned[0]) { //|| subjectAligned[0] == 'N')
327 ms[i][0].score = matchScore;
328 // ms[i][0].mismatches = 0;
331 // ms[i][0].mismatches = 1;
335 //fill rest of matrix
336 for (int j = 1; j < numCols; j++) { //iterate through matrix columns
338 for (int i = 0; i < numRows; i++) { //iterate through matrix rows
340 string subjectAligned = seqs[i]->getAligned();
342 int matchMisMatchScore = 0;
344 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
346 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
347 //matchMisMatchScore = matchScore;
349 }else if (queryAligned[j] == subjectAligned[j]) {
350 matchMisMatchScore = matchScore;
351 // ms[i][j].mismatches = ms[i][j-1].mismatches;
352 }else if (queryAligned[j] != subjectAligned[j]) {
353 matchMisMatchScore = misMatchPenalty;
354 // ms[i][j].mismatches = ms[i][j-1].mismatches + 1;
357 //compute score based on previous columns scores
358 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
360 int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
362 //you are not at yourself
363 if (prevIndex != i) { sumScore += penalty; }
364 if (sumScore < 0) { sumScore = 0; }
366 if (sumScore > ms[i][j].score) {
367 ms[i][j].score = sumScore;
368 ms[i][j].prev = prevIndex;
374 /* for(int i=0;i<numRows;i++){
375 cout << seqs[i]->getName();
376 for(int j=0;j<numCols;j++){
377 cout << '\t' << ms[i][j].mismatches;
382 /*cout << numRows << '\t' << numCols << endl;
383 for(int i=0;i<numRows;i++){
384 cout << seqs[i]->getName() << endl << seqs[i]->getAligned() << endl << endl;
385 if ((seqs[i]->getName() == "S000003470") || (seqs[i]->getName() == "S000383265") || (seqs[i]->getName() == "7000004128191054")) {
386 for(int j=0;j<numCols;j++){
387 cout << '\t' << ms[i][j].score;
393 /*for(int i=0;i<numRows;i++){
394 cout << seqs[i]->getName();
395 for(int j=0;j<numCols;j++){
396 cout << '\t' << ms[i][j].prev;
403 catch(exception& e) {
404 m->errorOut(e, "Maligner", "fillScoreMatrix");
408 //***************************************************************************************************************
409 vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
412 //get matrix dimensions
413 int numCols = query->getAligned().length();
414 int numRows = ms.size();
417 //find highest score scoring matrix
418 score_struct highestStruct;
419 int highestScore = 0;
421 for (int i = 0; i < numRows; i++) {
422 for (int j = 0; j < numCols; j++) {
423 if (ms[i][j].score > highestScore) {
424 highestScore = ms[i][j].score;
425 highestStruct = ms[i][j];
430 vector<score_struct> path;
432 int rowIndex = highestStruct.row;
433 int pos = highestStruct.col;
434 int score = highestStruct.score;
436 while (pos >= 0 && score > 0) {
437 score_struct temp = ms[rowIndex][pos];
440 if (score > 0) { path.push_back(temp); }
442 rowIndex = temp.prev;
446 reverse(path.begin(), path.end());
451 catch(exception& e) {
452 m->errorOut(e, "Maligner", "extractHighestPath");
456 //***************************************************************************************************************
457 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
459 vector<trace_struct> trace;
461 int region_index = path[0].row;
462 int region_start = path[0].col;
464 for (int i = 1; i < path.size(); i++) {
466 int next_region_index = path[i].row;
468 if (next_region_index != region_index) {
471 int col_index = path[i].col;
473 temp.col = region_start;
474 temp.oldCol = col_index-1;
475 temp.row = region_index;
477 trace.push_back(temp);
479 region_index = path[i].row;
480 region_start = col_index;
486 temp.col = region_start;
487 temp.oldCol = path[path.size()-1].col;
488 temp.row = region_index;
489 trace.push_back(temp);
494 catch(exception& e) {
495 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
500 /***************************************************************************************************************
502 vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
506 //get matrix dimensions
507 int numCols = query->getAligned().length();
508 int numRows = ms.size();
511 //find highest score scoring matrix
512 vector<score_struct> highestStruct;
513 int highestScore = 0;
515 for (int i = 0; i < numRows; i++) {
516 for (int j = 0; j < numCols; j++) {
517 if (ms[i][j].score > highestScore) {
518 highestScore = ms[i][j].score;
519 highestStruct.resize(0);
520 highestStruct.push_back(ms[i][j]);
522 else if(ms[i][j].score == highestScore){
523 highestStruct.push_back(ms[i][j]);
528 //cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl;
530 vector<trace_struct> maxTrace;
531 double maxPercentIdenticalQueryAntiChimera = 0;
533 for(int i=0;i<highestStruct.size();i++){
535 vector<score_struct> path;
537 int rowIndex = highestStruct[i].row;
538 int pos = highestStruct[i].col;
539 int score = highestStruct[i].score;
541 while (pos >= 0 && score > 0) {
542 score_struct temp = ms[rowIndex][pos];
545 if (score > 0) { path.push_back(temp); }
547 rowIndex = temp.prev;
551 reverse(path.begin(), path.end());
553 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
555 //cout << "traces\n";
556 //for(int j=0;j<trace.size();j++){
557 // cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
560 int traceStart = path[0].col;
561 int traceEnd = path[path.size()-1].col;
562 // cout << "traceStart/End\t" << traceStart << '\t' << traceEnd << endl;
564 string queryInRange = query->getAligned();
565 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
566 // cout << "here" << endl;
567 string chimeraSeq = constructChimericSeq(trace, refSeqs);
568 string antiChimeraSeq = constructAntiChimericSeq(trace, refSeqs);
570 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
571 double percentIdenticalQueryAntiChimera = computePercentID(queryInRange, antiChimeraSeq);
572 // cout << i << '\t' << percentIdenticalQueryChimera << '\t' << percentIdenticalQueryAntiChimera << endl;
574 if(percentIdenticalQueryAntiChimera > maxPercentIdenticalQueryAntiChimera){
575 maxPercentIdenticalQueryAntiChimera = percentIdenticalQueryAntiChimera;
579 // cout << maxPercentIdenticalQueryAntiChimera << endl;
584 catch(exception& e) {
585 m->errorOut(e, "Maligner", "extractHighestPath");
590 //***************************************************************************************************************
592 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
594 vector<trace_struct> trace;
596 int region_index = path[0].row;
597 int region_start = path[0].col;
599 for (int i = 1; i < path.size(); i++) {
601 int next_region_index = path[i].row;
603 if (next_region_index != region_index) {
606 int col_index = path[i].col;
608 temp.col = region_start;
609 temp.oldCol = col_index-1;
610 temp.row = region_index;
612 trace.push_back(temp);
614 region_index = path[i].row;
615 region_start = col_index;
621 temp.col = region_start;
622 temp.oldCol = path[path.size()-1].col;
623 temp.row = region_index;
624 trace.push_back(temp);
627 // cout << trace.size() << endl;
628 // for(int i=0;i<trace.size();i++){
629 // cout << seqs[trace[i].row]->getName() << endl;
636 catch(exception& e) {
637 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
642 //***************************************************************************************************************
644 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
648 for (int i = 0; i < trace.size(); i++) {
649 // cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
651 string seqAlign = seqs[trace[i].row]->getAligned();
652 seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
656 if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); }
659 catch(exception& e) {
660 m->errorOut(e, "Maligner", "constructChimericSeq");
665 //***************************************************************************************************************
667 string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
669 string antiChimera = "";
671 for (int i = 0; i < trace.size(); i++) {
672 // cout << i << '\t' << (trace.size() - i - 1) << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
674 int oppositeIndex = trace.size() - i - 1;
676 string seqAlign = seqs[trace[oppositeIndex].row]->getAligned();
677 seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
678 antiChimera += seqAlign;
683 catch(exception& e) {
684 m->errorOut(e, "Maligner", "constructChimericSeq");
689 //***************************************************************************************************************
690 float Maligner::computePercentID(string queryAlign, string chimera) {
693 if (queryAlign.length() != chimera.length()) {
694 m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
695 m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine();
696 m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
700 int numIdentical = 0;
703 for (int i = 0; i < queryAlign.length(); i++) {
704 if (((queryAlign[i] != 'G') && (queryAlign[i] != 'T') && (queryAlign[i] != 'A') && (queryAlign[i] != 'C')&& (queryAlign[i] != '.') && (queryAlign[i] != '-')) ||
705 ((chimera[i] != 'G') && (chimera[i] != 'T') && (chimera[i] != 'A') && (chimera[i] != 'C')&& (chimera[i] != '.') && (chimera[i] != '-'))) {}
708 bool charA = false; bool charB = false;
709 if ((queryAlign[i] == 'G') || (queryAlign[i] == 'T') || (queryAlign[i] == 'A') || (queryAlign[i] == 'C')) { charA = true; }
710 if ((chimera[i] == 'G') || (chimera[i] == 'T') || (chimera[i] == 'A') || (chimera[i] == 'C')) { charB = true; }
712 if (charA || charB) {
714 if (charA) { countA++; }
715 if (charB) { countB++; }
717 if (queryAlign[i] == chimera[i]) {
724 float numBases = (countA + countB) /(float) 2;
726 if (numBases == 0) { return 0; }
728 float percentIdentical = (numIdentical/(float)numBases) * 100;
730 return percentIdentical;
733 catch(exception& e) {
734 m->errorOut(e, "Maligner", "computePercentID");
738 //***************************************************************************************************************