some minor changes

[mothur.git] / maligner.cpp

diff --git a/maligner.cpp b/maligner.cpp
index 76689e9f11bf1ab3547ef9a33053830a7b436d14..30e993678a92895146ac289f132ab9dc2f7ac73b 100644 (file)
--- a/maligner.cpp
+++ b/maligner.cpp
@@ -34,6 +34,8 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
                        refSeqs = getKmerSeqs(query, numWanted); //fills indexes
                }else { m->mothurOut("not valid search."); exit(1);  } //should never get here
                
+               if (m->control_pressed) { return chimera;  }
+               
                refSeqs = minCoverageFilter(refSeqs);
 
                if (refSeqs.size() < 2)  { 
@@ -46,7 +48,8 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
 
                //fills outputResults
                chimera = chimeraMaligner(chimeraPenalty, decalc);
-       
+               
+               if (m->control_pressed) { return chimera;  }
                                
                //free memory
                delete query;
@@ -68,6 +71,9 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
                
                //trims seqs to first non gap char in all seqs and last non gap char in all seqs
                spotMap = decalc->trimSeqs(query, refSeqs);
+               
+               //you trimmed the whole sequence, skip
+               if (query->getAligned() == "") { return "no"; }
 
                vector<Sequence*> temp = refSeqs;
                temp.push_back(query);
@@ -77,10 +83,14 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
 
                vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
                
+               if (m->control_pressed) { return chimera;  }
+               
                fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
        
                vector<score_struct> path = extractHighestPath(matrix);
                
+               if (m->control_pressed) { return chimera;  }
+               
                vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
        
                if (trace.size() > 1) {         chimera = "yes";        }
@@ -95,6 +105,8 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
        
                percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
                
+               if (m->control_pressed) { return chimera;  }
+               
                //save output results
                for (int i = 0; i < trace.size(); i++) {
                        int regionStart = trace[i].col;
@@ -166,6 +178,8 @@ vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
                        //if coverage above minimum
                        if (coverage > minCoverage) {
                                newRefs.push_back(ref[i]);
+                       }else {
+                               delete ref[i];
                        }
                }
                
@@ -487,24 +501,19 @@ vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
        try {   
                indexes.clear();
                vector<Sequence*> refResults;
-               //generate blastdb
-               Database* database = new BlastDB(-2.0, -1.0, matchScore, misMatchPenalty);
-               for (int i = 0; i < db.size(); i++) {   database->addSequence(*db[i]);  }
-               database->generateDB();
-               database->setNumSeqs(db.size());
-               
+                               
                //get parts of query
                string queryUnAligned = q->getUnaligned();
                string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
                string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
 
-               Sequence* queryLeft = new Sequence(q->getName(), leftQuery);
-               Sequence* queryRight = new Sequence(q->getName(), rightQuery);
-               
-               vector<int> tempIndexesRight = database->findClosestMegaBlast(queryRight, num+1);
-               vector<int> tempIndexesLeft = database->findClosestMegaBlast(queryLeft, num+1);
+               Sequence* queryLeft = new Sequence(q->getName()+"left", leftQuery);
+               Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery);
                
-               //if ((tempIndexesRight.size() != (num+1)) || (tempIndexesLeft.size() != (num+1)))  {  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 porcess sequence " + q->getName()); m->mothurOutEndLine(); return refResults; }
+               vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
+               vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
+                       
+               //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; }
                
                vector<int> smaller;
                vector<int> larger;
@@ -555,8 +564,7 @@ vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
 //cout << endl;                
                delete queryRight;
                delete queryLeft;
-               delete database;
-               
+                       
                return refResults;
        }
        catch(exception& e) {
@@ -616,7 +624,7 @@ vector<Sequence*> Maligner::getKmerSeqs(Sequence* q, int num) {
                return refResults;
        }
        catch(exception& e) {
-               m->errorOut(e, "Maligner", "getBlastSeqs");
+               m->errorOut(e, "Maligner", "getKmerSeqs");
                exit(1);
        }
 }