changes for chimera slayer

[mothur.git] / slayer.cpp

diff --git a/slayer.cpp b/slayer.cpp
index 0ef795309a3f758c4f20e47e2e7504386aeaeb80..a244ea42634fc9b31218aceaafd3d76f249d92e8 100644 (file)
--- a/slayer.cpp
+++ b/slayer.cpp
@@ -28,24 +28,18 @@ string Slayer::getResults(Sequence* query, vector<Sequence*> refSeqs) {
                                Sequence* q = new Sequence(query->getName(), query->getAligned());
                                Sequence* leftParent = new Sequence(refSeqs[i]->getName(), refSeqs[i]->getAligned());
                                Sequence* rightParent = new Sequence(refSeqs[j]->getName(), refSeqs[j]->getAligned());
-               //cout << "parents: (" << refSeqs[i]->getName() << ", " << refSeqs[j]->getName() << ")\n";
+       
                                map<int, int> spots;  //map from spot in original sequence to spot in filtered sequence for query and both parents
                                vector<data_struct> divs = runBellerophon(q, leftParent, rightParent, spots);
-                               //cout << divs.size() << endl;
-                               if (m->control_pressed) { 
-                                       delete q;
-                                       delete leftParent;
-                                       delete rightParent;
-                                       return "no"; 
-                               }
+       
+                               if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; }
                                        
-//                             cout << divs.size() << endl;
                                vector<data_struct> selectedDivs;
                                for (int k = 0; k < divs.size(); k++) {
                                        
                                        vector<snps> snpsLeft = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winLStart, divs[k].winLEnd);
                                        vector<snps> snpsRight = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winRStart, divs[k].winREnd);
-                                       //cout << refSeqs[i]->getName() << '\t' << refSeqs[j]->getName() << '\t' << k << divs[k].parentA.getAligned() << endl << divs[k].parentB.getAligned() << endl;  
+       
                                        if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; }
                                        
                                        int numSNPSLeft = snpsLeft.size();
@@ -432,21 +426,41 @@ float Slayer::snpAB(vector<snps> data) {
        }
 }
 /***********************************************************************/
-float Slayer::computePercentID(string queryFrag, string parent, int left, int right) {
+float Slayer::computePercentID(string queryAlign, string chimera, int left, int right) {
        try {
-               int total = 0;
-               int matches = 0;
-
+                               
+               int numIdentical = 0;
+               int countA = 0;
+               int countB = 0;
                for (int i = left; i <= right; i++) {
-                       total++;
-                       if (queryFrag[i] == parent[i]) {
-                               matches++;
+                       if (((queryAlign[i] != 'G') && (queryAlign[i] != 'T') && (queryAlign[i] != 'A') && (queryAlign[i] != 'C')&& (queryAlign[i] != '.') && (queryAlign[i] != '-')) ||
+                               ((chimera[i] != 'G') && (chimera[i] != 'T') && (chimera[i] != 'A') && (chimera[i] != 'C')&& (chimera[i] != '.') && (chimera[i] != '-'))) {}
+                       else {
+                               
+                               bool charA = false; bool charB = false;
+                               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++; }
+                                       if (charB) { countB++; }
+                                       
+                                       if (queryAlign[i] == chimera[i]) {
+                                               numIdentical++;
+                                       }
+                               }
                        }
                }
-
-               float percentID =( matches/(float)total) * 100;
                
-               return percentID;
+               float numBases = (countA + countB) /(float) 2;
+               
+               if (numBases == 0) { return 0; }
+               
+               float percentIdentical = (numIdentical/(float)numBases) * 100;
+               
+               return percentIdentical;
+               
        }
        catch(exception& e) {
                m->errorOut(e, "Slayer", "computePercentID");