working on slayer bug

diff --git a/blastdb.cpp b/blastdb.cpp
index 026c4ca3054f0d9db16eb53ec1327115f67239d4..f162cdfd254d084e13ae14f3ff866db4d0805c1e 100644 (file)
--- a/blastdb.cpp
+++ b/blastdb.cpp
@@ -114,6 +114,7 @@ vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
 vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n) {
        try{
                vector<int> topMatches;
+               Scores.clear();
                
                ofstream queryFile;
 
@@ -146,7 +147,7 @@ vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n) {
                        
                        m->gobble(m8FileHandle);
                        topMatches.push_back(templateAccession);
-                       megaScores.push_back(searchScore);
+                       Scores.push_back(searchScore);
 //cout << templateAccession << endl;
                }
                m8FileHandle.close();

diff --git a/chimera.cpp b/chimera.cpp
index ae9ad9a46ef1ce5634c867219904857731bf3d6e..54124872b7a5e87ce4aa545585fcb4b94c759044 100644 (file)
--- a/chimera.cpp
+++ b/chimera.cpp
@@ -186,6 +186,8 @@ vector<Sequence*> Chimera::readSeqs(string file) {
        
                m->mothurOut("Done."); m->mothurOutEndLine();
                
+               filterString = (string(container[0]->getAligned().length(), '1'));
+               
                return container;
        }
        catch(exception& e) {

diff --git a/chimera.h b/chimera.h
index f82734675dcbe23f7d38e7465b42d3b6fe4e0239..58d637ec9698fa33d7c522c27cbc539595d42e5f 100644 (file)
--- a/chimera.h
+++ b/chimera.h
@@ -136,7 +136,7 @@ class Chimera {
 
        public:
        
-               Chimera(){ m = MothurOut::getInstance(); length = 0; unaligned = false; }
+               Chimera(){ m = MothurOut::getInstance(); length = 0; unaligned = false;  }
                virtual ~Chimera(){     for (int i = 0; i < templateSeqs.size(); i++) { delete templateSeqs[i];  } };
                virtual bool getUnaligned()                             {       return unaligned;                       }
                virtual int getLength()                                 {   return length;      }

diff --git a/chimeraslayer.cpp b/chimeraslayer.cpp
index 73037c783be84270a28b0dacf1795c51a3668d0e..7aa7cd4271f0ccbdcb443ae580b7ddd046b255cf 100644 (file)
--- a/chimeraslayer.cpp
+++ b/chimeraslayer.cpp
@@ -85,22 +85,22 @@ ChimeraSlayer::ChimeraSlayer(string file, string temp, bool trim, map<string, in
 //***************************************************************************************************************
 int ChimeraSlayer::doPrep() {
        try {
+               if (searchMethod == "distance") { 
+                       //read in all query seqs
+                       vector<Sequence*> tempQuerySeqs = readSeqs(fastafile);
                
-               //read in all query seqs
-               vector<Sequence*> tempQuerySeqs = readSeqs(fastafile);
+                       vector<Sequence*> temp = templateSeqs;
+                       for (int i = 0; i < tempQuerySeqs.size(); i++) {  temp.push_back(tempQuerySeqs[i]);  }
                
-               vector<Sequence*> temp = templateSeqs;
-               for (int i = 0; i < tempQuerySeqs.size(); i++) {  temp.push_back(tempQuerySeqs[i]);  }
+                       createFilter(temp, 0.0); //just removed columns where all seqs have a gap
                
-               createFilter(temp, 0.0); //just removed columns where all seqs have a gap
+                       for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i];  }
                
-               for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i];  }
-               
-               if (m->control_pressed) {  return 0; } 
-               
-               //run filter on template
-               for (int i = 0; i < templateSeqs.size(); i++) {  if (m->control_pressed) {  return 0; }  runFilter(templateSeqs[i]);  }
+                       if (m->control_pressed) {  return 0; } 
                
+                       //run filter on template
+                       for (int i = 0; i < templateSeqs.size(); i++) {  if (m->control_pressed) {  return 0; }  runFilter(templateSeqs[i]);  }
+               }
                string  kmerDBNameLeft;
                string  kmerDBNameRight;
        
@@ -751,7 +751,12 @@ int ChimeraSlayer::getChimeras(Sequence* query) {
                if (m->control_pressed) {  return 0;  }
                
                vector<results> Results = maligner.getOutput();
-
+               
+               //cout << query->getName() << endl;
+               //for (int i = 0; i < Results.size(); i++) {
+                       //cout << Results[i].parent << '\t' << Results[i].regionStart << '\t' << Results[i].regionEnd << '\t' << Results[i].nastRegionStart << '\t' << Results[i].nastRegionEnd << '\t' << Results[i].queryToParent << '\t' << Results[i].queryToParentLocal << endl;
+               //}
+               //cout << "done\n" << endl;
                if (chimeraFlag == "yes") {
                        
                        if (realign) {

diff --git a/cluster.cpp b/cluster.cpp
index 440562c07bcfc19d426779e4f627b84f7b0df3a1..ac9f4482da12110796c54f15c4137681a3d29038 100644 (file)
--- a/cluster.cpp
+++ b/cluster.cpp
@@ -214,6 +214,7 @@ void Cluster::update(double& cutOFF){
                // The vector has to be traversed in reverse order to preserve the index
                // for faster removal in removeCell()
                for (int i=nRowCells-1;i>=0;i--) {
+                       //if you are not the smallCell
                        if (!((rowCells[i]->row == smallRow) && (rowCells[i]->column == smallCol))) {
                                if (rowCells[i]->row == smallRow) {
                                        search = rowCells[i]->column;

diff --git a/database.hpp b/database.hpp
index a31075b4fddd46ec6d02d32186723a0608b7784d..f9c0c4878c31a67e2df415387f5ee5cf3bbda939 100644 (file)
--- a/database.hpp
+++ b/database.hpp
@@ -52,7 +52,7 @@ public:
        virtual vector<int> findClosestSequences(Sequence*, int) = 0;  // returns indexes of n closest sequences to query
        virtual vector<int> findClosestMegaBlast(Sequence*, int){return results;}
        virtual float getSearchScore();
-       virtual vector<float> getMegaBlastSearchScores() { return megaScores; } //assumes you already called findClosestMegaBlast
+       virtual vector<float> getSearchScores() { return Scores; } //assumes you already called findClosestMegaBlast
        virtual int getLongestBase(); 
        virtual void readKmerDB(ifstream&){};
        virtual void setNumSeqs(int i) {        numSeqs = i;    }
@@ -64,7 +64,7 @@ protected:
        int numSeqs, longest;
        float searchScore;
        vector<int> results;
-       vector<float> megaScores;
+       vector<float> Scores;
 };
 /**************************************************************************************************/
 #endif

diff --git a/decalc.cpp b/decalc.cpp
index 2a1736a7344e312e4d2cdb4371fc943335d6970b..52607fcb668dd9d90b903007f6ca978755798ec5 100644 (file)
--- a/decalc.cpp
+++ b/decalc.cpp
@@ -796,8 +796,8 @@ vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*
                vector<SeqDist> dists;
                float lastRight = distsRight[0].dist;
                float lastLeft = distsLeft[0].dist;
-               int lasti = 0;
-               for (int i = 0; i < distsLeft.size(); i++) {
+               //int lasti = 0;
+               for (int i = 0; i < numWanted+1; i++) {
                        //add left if you havent already
                        it = seen.find(db[distsLeft[i].index]->getName());
                        if (it == seen.end()) {  
@@ -816,13 +816,13 @@ vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*
 //                             cout << "loop-right\t" << db[distsRight[i].index]->getName() << '\t' << distsRight[i].dist << endl;
                        }
                        
-                       if (dists.size() > numWanted) { lasti = i; break; } //you have enough results
+                       //if (dists.size() > numWanted) { lasti = i; break; } //you have enough results
                }
                
 //             cout << "lastLeft\t" << lastLeft << endl;
                
                //add in sequences with same distance as last sequence added
-               lasti++;
+       /*      lasti++;
                int i = lasti;
                while (i < distsLeft.size()) {  
                        if (distsLeft[i].dist == lastLeft) {
@@ -856,8 +856,8 @@ vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*
                        else { break; }
                        i++;
                }
-
-               numWanted = seen.size();
+*/
+               numWanted = dists.size();
                
                if (numWanted > dists.size()) { 
                        //m->mothurOut("numwanted is larger than the number of template sequences, adjusting numwanted."); m->mothurOutEndLine(); 

diff --git a/distancedb.cpp b/distancedb.cpp
index 8d0c6298578c4311a59407610ae28a87f48ee493..c1bf7e7fb08dc3fae5b3cf17f575580269677148 100644 (file)
--- a/distancedb.cpp
+++ b/distancedb.cpp
@@ -49,6 +49,7 @@ void DistanceDB::addSequence(Sequence seq) {
 vector<int> DistanceDB::findClosestSequences(Sequence* query, int numWanted){
        try {
                vector<int> topMatches;
+               Scores.clear();
                bool templateSameLength = true;
                string sequence = query->getAligned();
                vector<seqDist> dists; 
@@ -87,6 +88,7 @@ vector<int> DistanceDB::findClosestSequences(Sequence* query, int numWanted){
                                //fill topmatches with numwanted closest sequences indexes
                                for (int i = 0; i < numWanted; i++) {
                                        topMatches.push_back(dists[i].seq2);
+                                       Scores.push_back(dists[i].dist);
                                }
                        }else {
                                int bestIndex = 0;
@@ -103,6 +105,7 @@ vector<int> DistanceDB::findClosestSequences(Sequence* query, int numWanted){
                                }
                                searchScore = smallDist;
                                topMatches.push_back(bestIndex);
+                               Scores.push_back(smallDist);
                        }
                
                }else{

diff --git a/kmerdb.cpp b/kmerdb.cpp
index e7e8ab2acbbf04c8bc30c6adb84fbd3880f7b2d8..2703e16d481ceccd6019ff642b97432ba5cf819f 100644 (file)
--- a/kmerdb.cpp
+++ b/kmerdb.cpp
@@ -60,6 +60,7 @@ vector<int> KmerDB::findClosestSequences(Sequence* candidateSeq, int num){
                vector<int> topMatches;
                Kmer kmer(kmerSize);
                searchScore = 0;
+               Scores.clear();
                
                vector<int> matches(numSeqs, 0);                                                //      a record of the sequences with shared kmers
                vector<int> timesKmerFound(kmerLocations.size()+1, 0);  //      a record of the kmers that we have already found
@@ -92,6 +93,8 @@ vector<int> KmerDB::findClosestSequences(Sequence* candidateSeq, int num){
                        //save top matches
                        for (int i = 0; i < num; i++) {
                                topMatches.push_back(seqMatches[i].seq);
+                               float thisScore = 100 * seqMatches[i].match / (float) numKmers;
+                               Scores.push_back(thisScore);
                        }
                }else{
                        int bestIndex = 0;
@@ -107,6 +110,7 @@ vector<int> KmerDB::findClosestSequences(Sequence* candidateSeq, int num){
                        searchScore = bestMatch;
                        searchScore = 100 * searchScore / (float) numKmers;             //      return the Sequence object corresponding to the db
                        topMatches.push_back(bestIndex);
+                       Scores.push_back(searchScore);
                }
                return topMatches;              
        }

diff --git a/maligner.cpp b/maligner.cpp
index 7d7145f4ca6bd9fac7a316745ace4432e4c26fe0..c5cc83a4d216b9c5a89eb6f60de8a15e5906d06c 100644 (file)
--- a/maligner.cpp
+++ b/maligner.cpp
@@ -48,7 +48,9 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
                if (m->control_pressed) { return chimera;  }
                
                refSeqs = minCoverageFilter(refSeqs);
-
+               
+               
+               
                if (refSeqs.size() < 2)  { 
                        for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
                        percentIdenticalQueryChimera = 0.0;
@@ -56,7 +58,7 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
                }
                
                int chimeraPenalty = computeChimeraPenalty();
-
+               //cout << "chimeraPenalty = " << chimeraPenalty << endl;
                //fills outputResults
                chimera = chimeraMaligner(chimeraPenalty, decalc);
                
@@ -89,8 +91,9 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
                vector<Sequence*> temp = refSeqs;
                temp.push_back(query);
                        
-                       
                verticalFilter(temp);
+               
+               //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
                
@@ -100,10 +103,10 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
        
                vector<trace_struct> trace = extractHighestPath(matrix);
                                
-//             cout << "traces\n";
-//             for(int i=0;i<trace.size();i++){
-//                     cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
-//             }
+               //cout << "traces\n";
+               //for(int i=0;i<trace.size();i++){
+               //      cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
+               //}
                
                if (trace.size() > 1) {         chimera = "yes";        }
                else { chimera = "no";  return chimera; }
@@ -315,7 +318,7 @@ void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequenc
                        if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
                                ms[i][0].score = 0;
 //                             ms[i][0].mismatches = 0;
-                       }else if (queryAligned[0] == subjectAligned[0] || subjectAligned[0] == 'N') {
+                       }else if (queryAligned[0] == subjectAligned[0])  { //|| subjectAligned[0] == 'N')
                                ms[i][0].score = matchScore;
 //                             ms[i][0].mismatches = 0;
                        }else{
@@ -336,7 +339,7 @@ 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;
+                                       //matchMisMatchScore = matchScore;
                                        //leave the same
                                }else if (queryAligned[j] == subjectAligned[j]) {
                                        matchMisMatchScore = matchScore;
@@ -363,29 +366,32 @@ void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequenc
                        }
                }
                
-//             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;
-//             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;
-//             }
+       /*      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;
+               }*/
                
                
        }
@@ -423,7 +429,7 @@ vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> >
                        }
                }
                        
-//             cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl;   
+               //cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl; 
                
                vector<trace_struct> maxTrace;
                double maxPercentIdenticalQueryAntiChimera = 0;
@@ -450,10 +456,10 @@ vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> >
 
                        vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
                
-//                     cout << "traces\n";
-//                     for(int j=0;j<trace.size();j++){
-//                             cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
-//                     }
+                       //cout << "traces\n";
+                       //for(int j=0;j<trace.size();j++){
+                       //      cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
+                       //}
                                                
                        int traceStart = path[0].col;
                        int traceEnd = path[path.size()-1].col; 
@@ -476,6 +482,7 @@ vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> >
                }
 //             cout << maxPercentIdenticalQueryAntiChimera << endl;
                return maxTrace;
+       
                
        }
        catch(exception& e) {
@@ -633,9 +640,9 @@ vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
                Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery);
 
                vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
-               vector<float> leftScores = databaseLeft->getMegaBlastSearchScores();
+               vector<float> leftScores = databaseLeft->getSearchScores();
                vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
-               vector<float> rightScores = databaseLeft->getMegaBlastSearchScores();
+               vector<float> rightScores = databaseLeft->getSearchScores();
 
                //if ((tempIndexesRight.size() == 0) && (tempIndexesLeft.size() == 0))  {  m->mothurOut("megablast returned " + toString(tempIndexesRight.size()) + " results for the right end, and " + toString(tempIndexesLeft.size()) + " for the left end. Needed " + toString(num+1) + ". Unable to process sequence " + q->getName()); m->mothurOutEndLine(); return refResults; }
                
@@ -676,11 +683,11 @@ vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
                        }
                        
                        lasti = i;
-                       if (mergedResults.size() > num) { break; }
+                       //if (mergedResults.size() > num) { break; }
                }
                
                //save lasti for smaller ties below
-               lasti++;
+               /*lasti++;
                int iSmaller = lasti;
                
                if (!(mergedResults.size() > num)) { //if we still need more results.  
@@ -726,7 +733,17 @@ vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
                        else { break; }
                        lasti++;                        
                }
+               */
                
+               for (int i = smaller.size(); i < larger.size(); i++) {
+                       //add right if you havent already
+                       it = seen.find(larger[i]);
+                       if (it == seen.end()) {  
+                               mergedResults.push_back(larger[i]);
+                               seen[larger[i]] = larger[i];
+                               lastLarger = largerScores[i];
+                       }
+               }
                numWanted = mergedResults.size();
                
                if (mergedResults.size() < numWanted) { numWanted = mergedResults.size(); }
@@ -738,6 +755,7 @@ vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
                                Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
                                refResults.push_back(temp);
                                indexes.push_back(mergedResults[i]);
+                               //cout << db[mergedResults[i]]->getName() << endl;
                        }
                        
 //cout << mergedResults[i] << endl;
@@ -768,11 +786,15 @@ vector<Sequence*> Maligner::getKmerSeqs(Sequence* q, int num) {
                
                vector<int> tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted);
                vector<int> tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted);
+               vector<float> scoresLeft = databaseLeft->getSearchScores();
+               vector<float> scoresRight = databaseRight->getSearchScores();
                
                //merge results         
                map<int, int> seen;
                map<int, int>::iterator it;
-               
+               float lastRight = scoresRight[0];
+               float lastLeft = scoresLeft[0];
+               //int lasti = 0;
                vector<int> mergedResults;
                for (int i = 0; i < tempIndexesLeft.size(); i++) {
                        //add left if you havent already
@@ -780,6 +802,7 @@ vector<Sequence*> Maligner::getKmerSeqs(Sequence* q, int num) {
                        if (it == seen.end()) {  
                                mergedResults.push_back(tempIndexesLeft[i]);
                                seen[tempIndexesLeft[i]] = tempIndexesLeft[i];
+                               lastLeft = scoresLeft[i];
                        }
 
                        //add right if you havent already
@@ -787,16 +810,61 @@ vector<Sequence*> Maligner::getKmerSeqs(Sequence* q, int num) {
                        if (it == seen.end()) {  
                                mergedResults.push_back(tempIndexesRight[i]);
                                seen[tempIndexesRight[i]] = tempIndexesRight[i];
+                               lastRight = scoresRight[i];
+                       }
+                       
+                       //if (mergedResults.size() > numWanted) { lasti = i; break; } //you have enough results
+               }
+               
+               //add in sequences with same distance as last sequence added
+               /*lasti++;
+               int i = lasti;
+               while (i < tempIndexesLeft.size()) {  
+                       if (scoresLeft[i] == lastLeft) {
+                               it = seen.find(tempIndexesLeft[i]);
+                               
+                               if (it == seen.end()) {  
+                                       mergedResults.push_back(tempIndexesLeft[i]);
+                                       seen[tempIndexesLeft[i]] = tempIndexesLeft[i];
+                               }
+                       }
+                       else { break; }
+                       i++;
+               }
+               
+               //              cout << "lastRight\t" << lastRight << endl;
+               //add in sequences with same distance as last sequence added
+               i = lasti;
+               while (i < tempIndexesRight.size()) {  
+                       if (scoresRight[i] == lastRight) {
+                               it = seen.find(tempIndexesRight[i]);
+                               
+                               if (it == seen.end()) {  
+                                       mergedResults.push_back(tempIndexesRight[i]);
+                                       seen[tempIndexesRight[i]] = tempIndexesRight[i];
+                               }
                        }
+                       else { break; }
+                       i++;
+               }*/
+               
+               numWanted = mergedResults.size();
+               
+               if (numWanted > mergedResults.size()) { 
+                       //m->mothurOut("numwanted is larger than the number of template sequences, adjusting numwanted."); m->mothurOutEndLine(); 
+                       numWanted = mergedResults.size();
                }
                
+               
 //cout << q->getName() << endl;                
                vector<Sequence*> refResults;
                for (int i = 0; i < numWanted; i++) {
 //cout << db[mergedResults[i]]->getName() << endl;     
-                       Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
-                       refResults.push_back(temp);
-                       indexes.push_back(mergedResults[i]);
+                       if (db[mergedResults[i]]->getName() != q->getName()) { 
+                               Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
+                               refResults.push_back(temp);
+                               indexes.push_back(mergedResults[i]);
+                       }
                }
 //cout << endl;                
                delete queryRight;

diff --git a/slayer.cpp b/slayer.cpp
index ba1bf337d4084e58416b06c48cb6a359bb5b05b1..8575faf4312450c78cc8db249259b036edb91e0d 100644 (file)
--- a/slayer.cpp
+++ b/slayer.cpp
@@ -297,8 +297,8 @@ int Slayer::bootstrapSNPS(vector<snps> left, vector<snps> right, float& BSA, flo
                int count_A = 0; // sceneario QLA,QRB supported
                int count_B = 0; // sceneario QLB,QRA supported
        
-               int numLeft = max(1, int(left.size() * (percentSNPSample/(float)100) + 0.5));
-               int numRight = max(1, int(right.size() * (percentSNPSample/(float)100) + 0.5));
+               int numLeft = max(1, int(left.size() * percentSNPSample/(float)100 + 0.5));
+               int numRight = max(1, int(right.size() * percentSNPSample/(float)100 + 0.5));
 
                for (int i = 0; i < iters; i++) {
                        //random sampling with replacement.
@@ -365,8 +365,8 @@ int Slayer::bootstrapSNPS(vector<snps> left, vector<snps> right, float& BSA, flo
 
 
 
-               BSA = ((float) count_A / (float) iters) * 100;
-               BSB = ((float) count_B / (float) iters) * 100;
+               BSA = (float) count_A / (float) iters * 100;
+               BSB = (float) count_B / (float) iters * 100;
 //cout << "bsa = " << BSA << " bsb = " << BSB << endl;
 
                return 0;