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;
77 // for(int i=0;i<refSeqs.size();i++){
78 // cout << refSeqs[i]->getName() << endl;
81 temp.push_back(query);
85 //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl ; }//<< refSeqs[i]->getAligned() << endl
87 vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
89 if (m->control_pressed) { return chimera; }
91 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
93 vector<score_struct> path = extractHighestPath(matrix);
95 if (m->control_pressed) { return chimera; }
97 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
99 if (trace.size() > 1) { chimera = "yes"; }
100 else { chimera = "no"; return chimera; }
102 int traceStart = path[0].col;
103 int traceEnd = path[path.size()-1].col;
104 string queryInRange = query->getAligned();
105 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
106 // cout << queryInRange << endl;
107 string chimeraSeq = constructChimericSeq(trace, refSeqs);
108 // cout << chimeraSeq << endl;
110 // cout << queryInRange.length() << endl;
111 // cout << chimeraSeq.length() << endl;
113 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
115 // cout << percentIdenticalQueryChimera << endl;
117 vector<trace_struct> trace = extractHighestPath(matrix);
119 //cout << "traces\n";
120 //for(int i=0;i<trace.size();i++){
121 // cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
124 if (trace.size() > 1) { chimera = "yes"; }
125 else { chimera = "no"; return chimera; }
127 int traceStart = trace[0].col;
128 int traceEnd = trace[trace.size()-1].oldCol;
129 string queryInRange = query->getAligned();
130 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart));*/
132 if (m->control_pressed) { return chimera; }
134 //save output results
135 for (int i = 0; i < trace.size(); i++) {
136 int regionStart = trace[i].col;
137 int regionEnd = trace[i].oldCol;
138 int seqIndex = trace[i].row;
140 // cout << regionStart << '\t' << regionEnd << '\t' << seqIndex << endl;
143 temp.parent = refSeqs[seqIndex]->getName();
144 temp.parentAligned = db[seqIndex]->getAligned();
145 temp.nastRegionStart = spotMap[regionStart];
146 temp.nastRegionEnd = spotMap[regionEnd];
147 temp.regionStart = regionStart;
148 temp.regionEnd = regionEnd;
150 string parentInRange = refSeqs[seqIndex]->getAligned();
151 parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
153 temp.queryToParent = computePercentID(queryInRange, parentInRange);
154 temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
156 string queryInRegion = query->getAligned();
157 queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
159 string parentInRegion = refSeqs[seqIndex]->getAligned();
160 parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
162 temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
164 // cout << query->getName() << '\t' << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << ", " << temp.divR << endl;
166 outputResults.push_back(temp);
171 catch(exception& e) {
172 m->errorOut(e, "Maligner", "chimeraMaligner");
176 /***********************************************************************/
177 //removes top matches that do not have minimum coverage with query.
178 vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
180 vector<Sequence*> newRefs;
182 string queryAligned = query->getAligned();
184 for (int i = 0; i < ref.size(); i++) {
186 string refAligned = ref[i]->getAligned();
192 for (int j = 0; j < queryAligned.length(); j++) {
194 if (isalpha(queryAligned[j])) {
197 if (isalpha(refAligned[j])) {
203 int coverage = ((numCovered/(float)numBases)*100);
205 //if coverage above minimum
206 if (coverage > minCoverage) {
207 newRefs.push_back(ref[i]);
215 catch(exception& e) {
216 m->errorOut(e, "Maligner", "minCoverageFilter");
220 /***********************************************************************/
221 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
222 int Maligner::computeChimeraPenalty() {
225 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
227 // if(numAllowable < 1){ numAllowable = 1; }
229 int penalty = int(numAllowable + 1) * misMatchPenalty;
234 catch(exception& e) {
235 m->errorOut(e, "Maligner", "computeChimeraPenalty");
239 /***********************************************************************/
240 //this is a vertical filter
241 void Maligner::verticalFilter(vector<Sequence*> seqs) {
243 vector<int> gaps; gaps.resize(query->getAligned().length(), 0);
245 string filterString = (string(query->getAligned().length(), '1'));
248 for (int i = 0; i < seqs.size(); i++) {
250 string seqAligned = seqs[i]->getAligned();
252 for (int j = 0; j < seqAligned.length(); j++) {
253 //if this spot is a gap
254 if ((seqAligned[j] == '-') || (seqAligned[j] == '.')) { gaps[j]++; }
258 //zero out spot where all sequences have blanks
259 int numColRemoved = 0;
260 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
261 if(gaps[i] == seqs.size()) { filterString[i] = '0'; numColRemoved++; }
264 map<int, int> newMap;
266 for (int i = 0; i < seqs.size(); i++) {
268 string seqAligned = seqs[i]->getAligned();
269 string newAligned = "";
272 for (int j = 0; j < seqAligned.length(); j++) {
273 //if this spot is not a gap
274 if (filterString[j] == '1') {
275 newAligned += seqAligned[j];
276 newMap[count] = spotMap[j];
281 seqs[i]->setAligned(newAligned);
286 catch(exception& e) {
287 m->errorOut(e, "Maligner", "verticalFilter");
291 //***************************************************************************************************************
292 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
295 vector< vector<score_struct> > m(rows);
297 for (int i = 0; i < rows; i++) {
298 for (int j = 0; j < cols; j++) {
300 //initialize each cell
303 temp.score = -9999999;
307 m[i].push_back(temp);
313 catch(exception& e) {
314 m->errorOut(e, "Maligner", "buildScoreMatrix");
319 //***************************************************************************************************************
321 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
324 //get matrix dimensions
325 int numCols = query->getAligned().length();
326 int numRows = seqs.size();
328 // cout << numRows << endl;
330 //initialize first col
331 string queryAligned = query->getAligned();
332 for (int i = 0; i < numRows; i++) {
333 string subjectAligned = seqs[i]->getAligned();
336 if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
338 // ms[i][0].mismatches = 0;
339 }else if (queryAligned[0] == subjectAligned[0]) { //|| subjectAligned[0] == 'N')
340 ms[i][0].score = matchScore;
341 // ms[i][0].mismatches = 0;
344 // ms[i][0].mismatches = 1;
348 //fill rest of matrix
349 for (int j = 1; j < numCols; j++) { //iterate through matrix columns
351 // for (int i = 0; i < 1; i++) { //iterate through matrix rows
352 for (int i = 0; i < numRows; i++) { //iterate through matrix rows
354 string subjectAligned = seqs[i]->getAligned();
356 int matchMisMatchScore = 0;
358 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
360 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
362 }else if (queryAligned[j] == subjectAligned[j]) {
363 matchMisMatchScore = matchScore;
364 }else if (queryAligned[j] != subjectAligned[j]) {
365 matchMisMatchScore = misMatchPenalty;
368 //compute score based on previous columns scores
369 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
371 int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
373 //you are not at yourself
374 if (prevIndex != i) { sumScore += penalty; }
375 if (sumScore < 0) { sumScore = 0; }
377 if (sumScore > ms[i][j].score) {
378 ms[i][j].score = sumScore;
379 ms[i][j].prev = prevIndex;
382 // cout << i << '\t' << j << '\t' << queryAligned[j] << '\t' << subjectAligned[j] << '\t' << matchMisMatchScore << '\t' << ms[i][j].score << endl;
392 // cout << numRows << '\t' << numCols << endl;
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].score;
402 // for(int i=0;i<numRows;i++){
403 // cout << seqs[i]->getName();
404 // for(int j=0;j<numCols;j++){
405 // cout << '\t' << ms[i][j].prev;
412 catch(exception& e) {
413 m->errorOut(e, "Maligner", "fillScoreMatrix");
417 //***************************************************************************************************************
418 vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
421 //get matrix dimensions
422 int numCols = query->getAligned().length();
423 int numRows = ms.size();
426 //find highest score scoring matrix
427 score_struct highestStruct;
428 int highestScore = 0;
430 for (int i = 0; i < numRows; i++) {
431 for (int j = 0; j < numCols; j++) {
432 if (ms[i][j].score > highestScore) {
433 highestScore = ms[i][j].score;
434 highestStruct = ms[i][j];
439 // cout << highestScore << endl;
440 vector<score_struct> path;
442 int rowIndex = highestStruct.row;
443 int pos = highestStruct.col;
444 int score = highestStruct.score;
446 while (pos >= 0 && score > 0) {
447 score_struct temp = ms[rowIndex][pos];
450 if (score > 0) { path.push_back(temp); }
452 rowIndex = temp.prev;
456 reverse(path.begin(), path.end());
461 catch(exception& e) {
462 m->errorOut(e, "Maligner", "extractHighestPath");
466 //***************************************************************************************************************
467 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
469 vector<trace_struct> trace;
471 int region_index = path[0].row;
472 int region_start = path[0].col;
474 for (int i = 1; i < path.size(); i++) {
476 int next_region_index = path[i].row;
477 // cout << i << '\t' << next_region_index << endl;
479 if (next_region_index != region_index) {
482 int col_index = path[i].col;
484 temp.col = region_start;
485 temp.oldCol = col_index-1;
486 temp.row = region_index;
488 trace.push_back(temp);
490 region_index = path[i].row;
491 region_start = col_index;
497 temp.col = region_start;
498 temp.oldCol = path[path.size()-1].col;
499 temp.row = region_index;
500 trace.push_back(temp);
505 catch(exception& e) {
506 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
511 /***************************************************************************************************************
513 vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
517 //get matrix dimensions
518 int numCols = query->getAligned().length();
519 int numRows = ms.size();
522 //find highest score scoring matrix
523 vector<score_struct> highestStruct;
524 int highestScore = 0;
526 for (int i = 0; i < numRows; i++) {
527 for (int j = 0; j < numCols; j++) {
528 if (ms[i][j].score > highestScore) {
529 highestScore = ms[i][j].score;
530 highestStruct.resize(0);
531 highestStruct.push_back(ms[i][j]);
533 else if(ms[i][j].score == highestScore){
534 highestStruct.push_back(ms[i][j]);
539 //cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl;
541 vector<trace_struct> maxTrace;
542 double maxPercentIdenticalQueryAntiChimera = 0;
544 for(int i=0;i<highestStruct.size();i++){
546 vector<score_struct> path;
548 int rowIndex = highestStruct[i].row;
549 int pos = highestStruct[i].col;
550 int score = highestStruct[i].score;
552 while (pos >= 0 && score > 0) {
553 score_struct temp = ms[rowIndex][pos];
556 if (score > 0) { path.push_back(temp); }
558 rowIndex = temp.prev;
562 reverse(path.begin(), path.end());
564 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
566 //cout << "traces\n";
567 //for(int j=0;j<trace.size();j++){
568 // cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
571 int traceStart = path[0].col;
572 int traceEnd = path[path.size()-1].col;
573 // cout << "traceStart/End\t" << traceStart << '\t' << traceEnd << endl;
575 string queryInRange = query->getAligned();
576 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
577 // cout << "here" << endl;
578 string chimeraSeq = constructChimericSeq(trace, refSeqs);
579 string antiChimeraSeq = constructAntiChimericSeq(trace, refSeqs);
581 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
582 double percentIdenticalQueryAntiChimera = computePercentID(queryInRange, antiChimeraSeq);
583 // cout << i << '\t' << percentIdenticalQueryChimera << '\t' << percentIdenticalQueryAntiChimera << endl;
585 if(percentIdenticalQueryAntiChimera > maxPercentIdenticalQueryAntiChimera){
586 maxPercentIdenticalQueryAntiChimera = percentIdenticalQueryAntiChimera;
590 // cout << maxPercentIdenticalQueryAntiChimera << endl;
595 catch(exception& e) {
596 m->errorOut(e, "Maligner", "extractHighestPath");
601 //***************************************************************************************************************
603 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
605 vector<trace_struct> trace;
607 int region_index = path[0].row;
608 int region_start = path[0].col;
610 for (int i = 1; i < path.size(); i++) {
612 int next_region_index = path[i].row;
614 if (next_region_index != region_index) {
617 int col_index = path[i].col;
619 temp.col = region_start;
620 temp.oldCol = col_index-1;
621 temp.row = region_index;
623 trace.push_back(temp);
625 region_index = path[i].row;
626 region_start = col_index;
632 temp.col = region_start;
633 temp.oldCol = path[path.size()-1].col;
634 temp.row = region_index;
635 trace.push_back(temp);
638 // cout << trace.size() << endl;
639 // for(int i=0;i<trace.size();i++){
640 // cout << seqs[trace[i].row]->getName() << endl;
647 catch(exception& e) {
648 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
653 //***************************************************************************************************************
655 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
659 for (int i = 0; i < trace.size(); i++) {
660 // cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
662 string seqAlign = seqs[trace[i].row]->getAligned();
663 seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
666 // cout << chimera << endl;
667 // if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); } //this was introducing a fence post error
668 // cout << chimera << endl;
671 catch(exception& e) {
672 m->errorOut(e, "Maligner", "constructChimericSeq");
677 //***************************************************************************************************************
679 string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
681 string antiChimera = "";
683 for (int i = 0; i < trace.size(); i++) {
684 // cout << i << '\t' << (trace.size() - i - 1) << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
686 int oppositeIndex = trace.size() - i - 1;
688 string seqAlign = seqs[trace[oppositeIndex].row]->getAligned();
689 seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
690 antiChimera += seqAlign;
695 catch(exception& e) {
696 m->errorOut(e, "Maligner", "constructChimericSeq");
701 //***************************************************************************************************************
702 float Maligner::computePercentID(string queryAlign, string chimera) {
705 if (queryAlign.length() != chimera.length()) {
706 m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
707 m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine();
708 m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
712 // cout << queryAlign.length() << endl;
713 int numIdentical = 0;
716 for (int i = 0; i < queryAlign.length(); i++) {
717 if (((queryAlign[i] != 'G') && (queryAlign[i] != 'T') && (queryAlign[i] != 'A') && (queryAlign[i] != 'C')&& (queryAlign[i] != '.') && (queryAlign[i] != '-')) ||
718 ((chimera[i] != 'G') && (chimera[i] != 'T') && (chimera[i] != 'A') && (chimera[i] != 'C')&& (chimera[i] != '.') && (chimera[i] != '-'))) {}
721 bool charA = false; bool charB = false;
722 if ((queryAlign[i] == 'G') || (queryAlign[i] == 'T') || (queryAlign[i] == 'A') || (queryAlign[i] == 'C')) { charA = true; }
723 if ((chimera[i] == 'G') || (chimera[i] == 'T') || (chimera[i] == 'A') || (chimera[i] == 'C')) { charB = true; }
726 if (charA || charB) {
728 if (charA) { countA++; }
729 if (charB) { countB++; }
731 if (queryAlign[i] == chimera[i]) {
735 // cout << queryAlign[i] << '\t' << chimera[i] << '\t' << countA << '\t' << countB << endl;
740 // cout << "pat\t" << countA << '\t' << countB << '\t' << numIdentical << endl;
743 float numBases = (countA + countB) /(float) 2;
745 if (numBases == 0) { return 0; }
747 // cout << numIdentical << '\t' << numBases << endl;
749 float percentIdentical = (numIdentical/(float)numBases) * 100;
751 // cout << percentIdentical << endl;
753 return percentIdentical;
756 catch(exception& e) {
757 m->errorOut(e, "Maligner", "computePercentID");
761 //***************************************************************************************************************