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.setName(q.getName()); query.setAligned(q.getAligned());
30 //copy refSeqs so that filter does not effect original
31 for(int i = 0; i < db.size(); i++) {
32 Sequence newSeq(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) {
68 //trims seqs to first non gap char in all seqs and last non gap char in all seqs
69 spotMap = decalc.trimSeqs(query, refSeqs);
71 //you trimmed the whole sequence, skip
72 if (query.getAligned() == "") { return "no"; }
74 vector<Sequence> temp = refSeqs;
75 temp.push_back(query);
77 temp = verticalFilter(temp);
78 query = temp[temp.size()-1];
79 for (int i = 0; i < temp.size()-1; i++) { refSeqs[i] = temp[i]; }
81 //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl ; }//<< refSeqs[i]->getAligned() << endl
83 vector< vector<score_struct> > matrix = buildScoreMatrix(query.getAligned().length(), refSeqs.size()); //builds and initializes
85 if (m->control_pressed) { return chimera; }
87 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
89 vector<score_struct> path = extractHighestPath(matrix);
91 if (m->control_pressed) { return chimera; }
93 vector<trace_struct> trace = mapTraceRegionsToAlignment(path);
95 if (trace.size() > 1) { chimera = "yes"; }
96 else { chimera = "no"; return chimera; }
98 int traceStart = path[0].col;
99 int traceEnd = path[path.size()-1].col;
100 string queryInRange = query.getAligned();
101 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
102 // cout << queryInRange << endl;
103 string chimeraSeq = constructChimericSeq(trace, refSeqs);
104 // cout << chimeraSeq << endl;
106 // cout << queryInRange.length() << endl;
107 // cout << chimeraSeq.length() << endl;
109 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
111 // cout << percentIdenticalQueryChimera << endl;
113 vector<trace_struct> trace = extractHighestPath(matrix);
115 //cout << "traces\n";
116 //for(int i=0;i<trace.size();i++){
117 // cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
120 if (trace.size() > 1) { chimera = "yes"; }
121 else { chimera = "no"; return chimera; }
123 int traceStart = trace[0].col;
124 int traceEnd = trace[trace.size()-1].oldCol;
125 string queryInRange = query->getAligned();
126 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart));*/
128 if (m->control_pressed) { return chimera; }
130 //save output results
131 for (int i = 0; i < trace.size(); i++) {
132 int regionStart = trace[i].col;
133 int regionEnd = trace[i].oldCol;
134 int seqIndex = trace[i].row;
136 // cout << regionStart << '\t' << regionEnd << '\t' << seqIndex << endl;
139 temp.parent = refSeqs[seqIndex].getName();
140 temp.parentAligned = db[seqIndex].getAligned();
141 temp.nastRegionStart = spotMap[regionStart];
142 temp.nastRegionEnd = spotMap[regionEnd];
143 temp.regionStart = unalignedMap[regionStart];
144 temp.regionEnd = unalignedMap[regionEnd];
146 string parentInRange = refSeqs[seqIndex].getAligned();
147 parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
149 temp.queryToParent = computePercentID(queryInRange, parentInRange);
150 temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
152 string queryInRegion = query.getAligned();
153 queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
155 string parentInRegion = refSeqs[seqIndex].getAligned();
156 parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
158 temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
160 //cout << query->getName() << '\t' << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << ", " << temp.divR << endl;
162 outputResults.push_back(temp);
167 catch(exception& e) {
168 m->errorOut(e, "Maligner", "chimeraMaligner");
172 /***********************************************************************/
173 //removes top matches that do not have minimum coverage with query.
174 vector<Sequence> Maligner::minCoverageFilter(vector<Sequence> ref){
176 vector<Sequence> newRefs;
178 string queryAligned = query.getAligned();
180 for (int i = 0; i < ref.size(); i++) {
182 string refAligned = ref[i].getAligned();
188 for (int j = 0; j < queryAligned.length(); j++) {
190 if (isalpha(queryAligned[j])) {
193 if (isalpha(refAligned[j])) {
199 int coverage = ((numCovered/(float)numBases)*100);
201 //if coverage above minimum
202 if (coverage > minCoverage) {
203 newRefs.push_back(ref[i]);
211 catch(exception& e) {
212 m->errorOut(e, "Maligner", "minCoverageFilter");
216 /***********************************************************************/
217 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
218 int Maligner::computeChimeraPenalty() {
221 int numAllowable = ((1.0 - (1.0/minDivR)) * query.getNumBases());
223 // if(numAllowable < 1){ numAllowable = 1; }
225 int penalty = int(numAllowable + 1) * misMatchPenalty;
230 catch(exception& e) {
231 m->errorOut(e, "Maligner", "computeChimeraPenalty");
235 /***********************************************************************/
236 //this is a vertical filter
237 vector<Sequence> Maligner::verticalFilter(vector<Sequence> seqs) {
239 vector<int> gaps; gaps.resize(query.getAligned().length(), 0);
241 string filterString = (string(query.getAligned().length(), '1'));
244 for (int i = 0; i < seqs.size(); i++) {
246 string seqAligned = seqs[i].getAligned();
248 for (int j = 0; j < seqAligned.length(); j++) {
249 //if this spot is a gap
250 if ((seqAligned[j] == '-') || (seqAligned[j] == '.')) { gaps[j]++; }
254 //zero out spot where all sequences have blanks
255 int numColRemoved = 0;
256 for(int i = 0; i < seqs[0].getAligned().length(); i++){
257 if(gaps[i] == seqs.size()) { filterString[i] = '0'; numColRemoved++; }
260 map<int, int> newMap;
262 for (int i = 0; i < seqs.size(); i++) {
264 string seqAligned = seqs[i].getAligned();
265 string newAligned = "";
268 for (int j = 0; j < seqAligned.length(); j++) {
269 //if this spot is not a gap
270 if (filterString[j] == '1') {
271 newAligned += seqAligned[j];
272 newMap[count] = spotMap[j];
277 seqs[i].setAligned(newAligned);
280 string query = seqs[seqs.size()-1].getAligned();
281 int queryLength = query.length();
283 unalignedMap.resize(queryLength, 0);
286 for(int i=1;i<queryLength;i++){
287 if(query[i] != '.' && query[i] != '-'){
288 unalignedMap[i] = unalignedMap[i-1] + 1;
291 unalignedMap[i] = unalignedMap[i-1];
299 catch(exception& e) {
300 m->errorOut(e, "Maligner", "verticalFilter");
304 //***************************************************************************************************************
305 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
308 vector< vector<score_struct> > m(rows);
310 for (int i = 0; i < rows; i++) {
311 for (int j = 0; j < cols; j++) {
313 //initialize each cell
316 temp.score = -9999999;
320 m[i].push_back(temp);
326 catch(exception& e) {
327 m->errorOut(e, "Maligner", "buildScoreMatrix");
332 //***************************************************************************************************************
334 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence> seqs, int penalty) {
337 //get matrix dimensions
338 int numCols = query.getAligned().length();
339 int numRows = seqs.size();
341 // cout << numRows << endl;
343 //initialize first col
344 string queryAligned = query.getAligned();
345 for (int i = 0; i < numRows; i++) {
346 string subjectAligned = seqs[i].getAligned();
349 if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
351 // ms[i][0].mismatches = 0;
352 }else if (queryAligned[0] == subjectAligned[0]) { //|| subjectAligned[0] == 'N')
353 ms[i][0].score = matchScore;
354 // ms[i][0].mismatches = 0;
357 // ms[i][0].mismatches = 1;
361 //fill rest of matrix
362 for (int j = 1; j < numCols; j++) { //iterate through matrix columns
364 // for (int i = 0; i < 1; i++) { //iterate through matrix rows
365 for (int i = 0; i < numRows; i++) { //iterate through matrix rows
367 string subjectAligned = seqs[i].getAligned();
369 int matchMisMatchScore = 0;
371 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
373 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
375 }else if (queryAligned[j] == subjectAligned[j]) {
376 matchMisMatchScore = matchScore;
377 }else if (queryAligned[j] != subjectAligned[j]) {
378 matchMisMatchScore = misMatchPenalty;
381 //compute score based on previous columns scores
382 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
384 int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
386 //you are not at yourself
387 if (prevIndex != i) { sumScore += penalty; }
388 if (sumScore < 0) { sumScore = 0; }
390 if (sumScore > ms[i][j].score) {
391 ms[i][j].score = sumScore;
392 ms[i][j].prev = prevIndex;
395 // cout << i << '\t' << j << '\t' << queryAligned[j] << '\t' << subjectAligned[j] << '\t' << matchMisMatchScore << '\t' << ms[i][j].score << endl;
405 // cout << numRows << '\t' << numCols << endl;
406 // for(int i=0;i<numRows;i++){
407 // cout << seqs[i]->getName();
408 // for(int j=0;j<numCols;j++){
409 // cout << '\t' << ms[i][j].score;
415 // for(int i=0;i<numRows;i++){
416 // cout << seqs[i]->getName();
417 // for(int j=0;j<numCols;j++){
418 // cout << '\t' << ms[i][j].prev;
425 catch(exception& e) {
426 m->errorOut(e, "Maligner", "fillScoreMatrix");
430 //***************************************************************************************************************
431 vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
434 //get matrix dimensions
435 int numCols = query.getAligned().length();
436 int numRows = ms.size();
439 //find highest score scoring matrix
440 score_struct highestStruct;
441 int highestScore = 0;
443 for (int i = 0; i < numRows; i++) {
444 for (int j = 0; j < numCols; j++) {
445 if (ms[i][j].score > highestScore) {
446 highestScore = ms[i][j].score;
447 highestStruct = ms[i][j];
452 vector<score_struct> path;
454 int rowIndex = highestStruct.row;
455 int pos = highestStruct.col;
456 int score = highestStruct.score;
458 // cout << rowIndex << '\t' << pos << '\t' << score << endl;
460 while (pos >= 0 && score > 0) {
461 score_struct temp = ms[rowIndex][pos];
464 if (score > 0) { path.push_back(temp); }
466 rowIndex = temp.prev;
470 reverse(path.begin(), path.end());
475 catch(exception& e) {
476 m->errorOut(e, "Maligner", "extractHighestPath");
480 //***************************************************************************************************************
481 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path) {
483 vector<trace_struct> trace;
485 int region_index = path[0].row;
486 int region_start = path[0].col;
488 for (int i = 1; i < path.size(); i++) {
490 int next_region_index = path[i].row;
491 //cout << i << '\t' << next_region_index << endl;
493 if (next_region_index != region_index) {
496 int col_index = path[i].col;
498 temp.col = region_start;
499 temp.oldCol = col_index-1;
500 temp.row = region_index;
502 trace.push_back(temp);
504 region_index = path[i].row;
505 region_start = col_index;
511 temp.col = region_start;
512 temp.oldCol = path[path.size()-1].col;
513 temp.row = region_index;
514 trace.push_back(temp);
519 catch(exception& e) {
520 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
525 /***************************************************************************************************************
527 vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
531 //get matrix dimensions
532 int numCols = query->getAligned().length();
533 int numRows = ms.size();
536 //find highest score scoring matrix
537 vector<score_struct> highestStruct;
538 int highestScore = 0;
540 for (int i = 0; i < numRows; i++) {
541 for (int j = 0; j < numCols; j++) {
542 if (ms[i][j].score > highestScore) {
543 highestScore = ms[i][j].score;
544 highestStruct.resize(0);
545 highestStruct.push_back(ms[i][j]);
547 else if(ms[i][j].score == highestScore){
548 highestStruct.push_back(ms[i][j]);
553 //cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl;
555 vector<trace_struct> maxTrace;
556 double maxPercentIdenticalQueryAntiChimera = 0;
558 for(int i=0;i<highestStruct.size();i++){
560 vector<score_struct> path;
562 int rowIndex = highestStruct[i].row;
563 int pos = highestStruct[i].col;
564 int score = highestStruct[i].score;
566 while (pos >= 0 && score > 0) {
567 score_struct temp = ms[rowIndex][pos];
570 if (score > 0) { path.push_back(temp); }
572 rowIndex = temp.prev;
576 reverse(path.begin(), path.end());
578 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
580 //cout << "traces\n";
581 //for(int j=0;j<trace.size();j++){
582 // cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
585 int traceStart = path[0].col;
586 int traceEnd = path[path.size()-1].col;
587 // cout << "traceStart/End\t" << traceStart << '\t' << traceEnd << endl;
589 string queryInRange = query->getAligned();
590 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
591 // cout << "here" << endl;
592 string chimeraSeq = constructChimericSeq(trace, refSeqs);
593 string antiChimeraSeq = constructAntiChimericSeq(trace, refSeqs);
595 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
596 double percentIdenticalQueryAntiChimera = computePercentID(queryInRange, antiChimeraSeq);
597 // cout << i << '\t' << percentIdenticalQueryChimera << '\t' << percentIdenticalQueryAntiChimera << endl;
599 if(percentIdenticalQueryAntiChimera > maxPercentIdenticalQueryAntiChimera){
600 maxPercentIdenticalQueryAntiChimera = percentIdenticalQueryAntiChimera;
604 // cout << maxPercentIdenticalQueryAntiChimera << endl;
609 catch(exception& e) {
610 m->errorOut(e, "Maligner", "extractHighestPath");
615 //***************************************************************************************************************
617 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
619 vector<trace_struct> trace;
621 int region_index = path[0].row;
622 int region_start = path[0].col;
624 for (int i = 1; i < path.size(); i++) {
626 int next_region_index = path[i].row;
628 if (next_region_index != region_index) {
631 int col_index = path[i].col;
633 temp.col = region_start;
634 temp.oldCol = col_index-1;
635 temp.row = region_index;
637 trace.push_back(temp);
639 region_index = path[i].row;
640 region_start = col_index;
646 temp.col = region_start;
647 temp.oldCol = path[path.size()-1].col;
648 temp.row = region_index;
649 trace.push_back(temp);
652 // cout << trace.size() << endl;
653 // for(int i=0;i<trace.size();i++){
654 // cout << seqs[trace[i].row]->getName() << endl;
661 catch(exception& e) {
662 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
667 //***************************************************************************************************************
669 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence> seqs) {
673 for (int i = 0; i < trace.size(); i++) {
674 // cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
676 string seqAlign = seqs[trace[i].row].getAligned();
677 seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
680 // cout << chimera << endl;
681 // if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); } //this was introducing a fence post error
682 // cout << chimera << endl;
685 catch(exception& e) {
686 m->errorOut(e, "Maligner", "constructChimericSeq");
691 //***************************************************************************************************************
693 string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence> seqs) {
695 string antiChimera = "";
697 for (int i = 0; i < trace.size(); i++) {
698 // cout << i << '\t' << (trace.size() - i - 1) << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
700 int oppositeIndex = trace.size() - i - 1;
702 string seqAlign = seqs[trace[oppositeIndex].row].getAligned();
703 seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
704 antiChimera += seqAlign;
709 catch(exception& e) {
710 m->errorOut(e, "Maligner", "constructChimericSeq");
715 //***************************************************************************************************************
716 float Maligner::computePercentID(string queryAlign, string chimera) {
719 if (queryAlign.length() != chimera.length()) {
720 m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
721 m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine();
722 m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
726 // cout << queryAlign.length() << endl;
727 int numIdentical = 0;
730 for (int i = 0; i < queryAlign.length(); i++) {
731 if (((queryAlign[i] != 'G') && (queryAlign[i] != 'T') && (queryAlign[i] != 'A') && (queryAlign[i] != 'C')&& (queryAlign[i] != '.') && (queryAlign[i] != '-')) ||
732 ((chimera[i] != 'G') && (chimera[i] != 'T') && (chimera[i] != 'A') && (chimera[i] != 'C')&& (chimera[i] != '.') && (chimera[i] != '-'))) {}
735 bool charA = false; bool charB = false;
736 if ((queryAlign[i] == 'G') || (queryAlign[i] == 'T') || (queryAlign[i] == 'A') || (queryAlign[i] == 'C')) { charA = true; }
737 if ((chimera[i] == 'G') || (chimera[i] == 'T') || (chimera[i] == 'A') || (chimera[i] == 'C')) { charB = true; }
740 if (charA || charB) {
742 if (charA) { countA++; }
743 if (charB) { countB++; }
745 if (queryAlign[i] == chimera[i]) {
749 // cout << queryAlign[i] << '\t' << chimera[i] << '\t' << countA << '\t' << countB << endl;
754 // cout << "pat\t" << countA << '\t' << countB << '\t' << numIdentical << endl;
757 float numBases = (countA + countB) /(float) 2;
759 if (numBases == 0) { return 0; }
761 // cout << numIdentical << '\t' << numBases << endl;
763 float percentIdentical = (numIdentical/(float)numBases) * 100;
765 // cout << percentIdentical << endl;
767 return percentIdentical;
770 catch(exception& e) {
771 m->errorOut(e, "Maligner", "computePercentID");
775 //***************************************************************************************************************