chimera.slayer debugging

[mothur.git] / maligner.cpp

diff --git a/maligner.cpp b/maligner.cpp
index 599c78be2c5c69570cb0d96febfa861f90bb8eb5..1205cef5e36b8bb4381f5ae8bb516c30e934318d 100644 (file)
--- a/maligner.cpp
+++ b/maligner.cpp
@@ -73,11 +73,16 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
                if (query->getAligned() == "") { return "no"; }
 
                vector<Sequence*> temp = refSeqs;
+               
+//             for(int i=0;i<refSeqs.size();i++){
+//                     cout << refSeqs[i]->getName() << endl;
+//             }
+               
                temp.push_back(query);
                        
                verticalFilter(temp);
                
-               //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl << refSeqs[i]->getAligned() << endl; }
+               //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl ; }//<< refSeqs[i]->getAligned() << endl
 
                vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
                
@@ -90,19 +95,24 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
                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));
-               
+               queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
+//             cout << queryInRange << endl;
                string chimeraSeq = constructChimericSeq(trace, refSeqs);
+//             cout << chimeraSeq << endl;
+               
+//             cout << queryInRange.length() << endl;
+//             cout << chimeraSeq.length() << endl;
                
                percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
                
+//             cout << percentIdenticalQueryChimera << endl;
        /*      
                vector<trace_struct> trace = extractHighestPath(matrix);
                                
@@ -127,6 +137,7 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
                        int regionEnd = trace[i].oldCol;
                        int seqIndex = trace[i].row;
                        
+//                     cout << regionStart << '\t' << regionEnd << '\t' << seqIndex << endl;
                        results temp;
                        
                        temp.parent = refSeqs[seqIndex]->getName();
@@ -137,20 +148,21 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
                        temp.regionEnd = regionEnd;
                        
                        string parentInRange = refSeqs[seqIndex]->getAligned();
-                       parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart));
+                       parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
                        
                        temp.queryToParent = computePercentID(queryInRange, parentInRange);
                        temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
 
                        string queryInRegion = query->getAligned();
-                       queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart));
+                       queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
                        
                        string parentInRegion = refSeqs[seqIndex]->getAligned();
-                       parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart));
+                       parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
                        
                        temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
                        
-                       cout << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << endl;
+//                     cout << query->getName() << '\t' << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << ", " <<  temp.divR << endl;
+
                        outputResults.push_back(temp);
                }
                
@@ -313,6 +325,8 @@ void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequenc
                int numCols = query->getAligned().length();
                int numRows = seqs.size();
                
+//             cout << numRows << endl;
+               
                //initialize first col
                string queryAligned = query->getAligned();
                for (int i = 0; i < numRows; i++) {
@@ -334,6 +348,7 @@ void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequenc
                //fill rest of matrix
                for (int j = 1; j < numCols; j++) {  //iterate through matrix columns
                
+//                     for (int i = 0; i < 1; i++) {  //iterate through matrix rows
                        for (int i = 0; i < numRows; i++) {  //iterate through matrix rows
                                
                                string subjectAligned = seqs[i]->getAligned();
@@ -343,16 +358,13 @@ void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequenc
                                if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
                                        //leave the same
                                }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
-                                       //matchMisMatchScore = matchScore;
                                        //leave the same
                                }else if (queryAligned[j] == subjectAligned[j]) {
                                        matchMisMatchScore = matchScore;
-//                                     ms[i][j].mismatches = ms[i][j-1].mismatches;
                                }else if (queryAligned[j] != subjectAligned[j]) {
                                        matchMisMatchScore = misMatchPenalty;
-//                                     ms[i][j].mismatches = ms[i][j-1].mismatches + 1;
                                }
-                               
+
                                //compute score based on previous columns scores
                                for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
                                        
@@ -367,35 +379,33 @@ void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequenc
                                                ms[i][j].prev = prevIndex;
                                        }
                                }
+//                             cout << i << '\t' << j << '\t' << queryAligned[j] << '\t' << subjectAligned[j] << '\t' << matchMisMatchScore << '\t' << ms[i][j].score << endl;
+
                        }
+                       
                }
                
-       /*      for(int i=0;i<numRows;i++){
-                       cout << seqs[i]->getName();
-                       for(int j=0;j<numCols;j++){
-                               cout << '\t' << ms[i][j].mismatches;
-                       }
-                       cout << endl;
-               }
-               cout << endl;*/
-               /*cout << numRows << '\t' << numCols << endl;
-               for(int i=0;i<numRows;i++){
-                       cout << seqs[i]->getName() << endl << seqs[i]->getAligned() << endl << endl;
-                       if ((seqs[i]->getName() == "S000003470") || (seqs[i]->getName() == "S000383265") || (seqs[i]->getName() == "7000004128191054")) {
-                       for(int j=0;j<numCols;j++){
-                               cout << '\t' << ms[i][j].score;
-                       }
-                       cout << endl;
-                       }
-               }
-               cout << endl;*/
-               /*for(int i=0;i<numRows;i++){
-                       cout << seqs[i]->getName();
-                       for(int j=0;j<numCols;j++){
-                               cout << '\t' << ms[i][j].prev;
-                       }
-                       cout << endl;
-               }*/
+               
+               
+               
+               
+//             cout << numRows << '\t' << numCols << endl;
+//             for(int i=0;i<numRows;i++){
+//                     cout << seqs[i]->getName();
+//                     for(int j=0;j<numCols;j++){
+//                             cout << '\t' << ms[i][j].score;
+//                     }
+//                     cout << endl;
+//             }
+//             cout << endl;
+//             
+//             for(int i=0;i<numRows;i++){
+//                     cout << seqs[i]->getName();
+//                     for(int j=0;j<numCols;j++){
+//                             cout << '\t' << ms[i][j].prev;
+//                     }
+//                     cout << endl;
+//             }
                
                
        }
@@ -426,6 +436,7 @@ vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> >
                        }
                }
                
+//             cout << highestScore << endl;
                vector<score_struct> path;
                
                int rowIndex = highestStruct.row;
@@ -463,7 +474,8 @@ vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> p
                for (int i = 1; i < path.size(); i++) {
                        
                        int next_region_index = path[i].row;
-                       
+//                     cout << i << '\t' << next_region_index << endl;
+
                        if (next_region_index != region_index) {
                                
                                // add trace region
@@ -651,8 +663,9 @@ string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequenc
                        seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
                        chimera += seqAlign;
                }
-               
-               if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); }
+//             cout << chimera << endl;
+//             if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); }       //this was introducing a fence post error
+//             cout << chimera << endl;
                return chimera;
        }
        catch(exception& e) {
@@ -696,6 +709,7 @@ float Maligner::computePercentID(string queryAlign, string chimera) {
                        return -1.0;
                }
 
+//             cout << queryAlign.length() << endl;
                int numIdentical = 0;
                int countA = 0;
                int countB = 0;
@@ -708,6 +722,7 @@ float Maligner::computePercentID(string queryAlign, string chimera) {
                                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++; }
@@ -717,15 +732,24 @@ float Maligner::computePercentID(string queryAlign, string chimera) {
                                                numIdentical++;
                                        }
                                }
+//                             cout << queryAlign[i] << '\t' << chimera[i] << '\t' << countA << '\t' << countB << endl;
+
                        }
                }
                
+//             cout << "pat\t" << countA << '\t' << countB << '\t' << numIdentical << endl;
+
+               
                float numBases = (countA + countB) /(float) 2;
                
                if (numBases == 0) { return 0; }
        
+//             cout << numIdentical << '\t' << numBases << endl;
+               
                float percentIdentical = (numIdentical/(float)numBases) * 100;
-
+               
+//             cout << percentIdentical << endl;
+               
                return percentIdentical;
                
        }