]> git.donarmstrong.com Git - mothur.git/blobdiff - bellerophon.cpp
testing mpi
[mothur.git] / bellerophon.cpp
index 9aa4f7ac7aa566f39a9b68c7b9e0103f516e0ba4..54dfb9b67ecf5e764e3c159dffcc47d6fd7004dc 100644 (file)
 
 //***************************************************************************************************************
 
-Bellerophon::Bellerophon(string name) {
+Bellerophon::Bellerophon(string name, string o)  {
        try {
                fastafile = name;
+               outputDir = o;
        }
        catch(exception& e) {
-               errorOut(e, "Bellerophon", "Bellerophon");
+               m->errorOut(e, "Bellerophon", "Bellerophon");
                exit(1);
        }
 }
 
 //***************************************************************************************************************
-void Bellerophon::print(ostream& out) {
+int Bellerophon::print(ostream& out, ostream& outAcc) {
        try {
                int above1 = 0;
                out << "Name\tScore\tLeft\tRight\t" << endl;
                //output prefenence structure to .chimeras file
                for (int i = 0; i < pref.size(); i++) {
+                       
+                       if (m->control_pressed) {  return 0; }
+                       
                        out << pref[i].name << '\t' << setprecision(3) << pref[i].score[0] << '\t' << pref[i].leftParent[0] << '\t' << pref[i].rightParent[0] << endl;
                        
                        //calc # of seqs with preference above 1.0
                        if (pref[i].score[0] > 1.0) { 
                                above1++; 
-                               mothurOut(pref[i].name + " is a suspected chimera at breakpoint " + toString(pref[i].midpoint)); mothurOutEndLine();
-                               mothurOut("It's score is " + toString(pref[i].score[0]) + " with suspected left parent " + pref[i].leftParent[0] + " and right parent " + pref[i].rightParent[0]); mothurOutEndLine();
+                               outAcc << pref[i].name << endl;
+                               m->mothurOut(pref[i].name + " is a suspected chimera at breakpoint " + toString(pref[i].midpoint)); m->mothurOutEndLine();
+                               m->mothurOut("It's score is " + toString(pref[i].score[0]) + " with suspected left parent " + pref[i].leftParent[0] + " and right parent " + pref[i].rightParent[0]); m->mothurOutEndLine();
                        }
                }
                
                //output results to screen
-               mothurOutEndLine();
-               mothurOut("Sequence with preference score above 1.0: " + toString(above1)); mothurOutEndLine();
+               m->mothurOutEndLine();
+               m->mothurOut("Sequence with preference score above 1.0: " + toString(above1)); m->mothurOutEndLine();
                int spot;
                spot = pref.size()-1;
-               mothurOut("Minimum:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+               m->mothurOut("Minimum:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
                spot = pref.size() * 0.975;
-               mothurOut("2.5%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+               m->mothurOut("2.5%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
                spot = pref.size() * 0.75;
-               mothurOut("25%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+               m->mothurOut("25%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
                spot = pref.size() * 0.50;
-               mothurOut("Median: \t" + toString(pref[spot].score[0])); mothurOutEndLine();
+               m->mothurOut("Median: \t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
                spot = pref.size() * 0.25;
-               mothurOut("75%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+               m->mothurOut("75%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
                spot = pref.size() * 0.025;
-               mothurOut("97.5%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+               m->mothurOut("97.5%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
                spot = 0;
-               mothurOut("Maximum:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+               m->mothurOut("Maximum:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
+               
+               return 1;
 
        }
        catch(exception& e) {
-               errorOut(e, "Bellerophon", "print");
+               m->errorOut(e, "Bellerophon", "print");
                exit(1);
        }
 }
@@ -75,18 +82,24 @@ inline bool comparePref(Preference left, Preference right){
 }
 
 //***************************************************************************************************************
-void Bellerophon::getChimeras() {
+int Bellerophon::getChimeras() {
        try {
                
                //do soft filter
                if (filter)  {
                        string optionString = "fasta=" + fastafile + ", soft=50";
+                       if (outputDir != "") { optionString += ", outputdir=" + outputDir; }
+                       
                        filterSeqs = new FilterSeqsCommand(optionString);
                        filterSeqs->execute();
                        delete filterSeqs;
                        
+                       if (m->control_pressed) { return 0; }
+                       
                        //reset fastafile to filtered file
-                       fastafile = getRootName(fastafile) + "filter.fasta";
+                       if (outputDir == "") { fastafile = getRootName(fastafile) + "filter.fasta"; }
+                       else                             { fastafile = outputDir + getRootName(getSimpleName(fastafile)) + "filter.fasta"; }
+                       
                }
                
                distCalculator = new eachGapDist();
@@ -94,22 +107,26 @@ void Bellerophon::getChimeras() {
                //read in sequences
                seqs = readSeqs(fastafile);
                
+               if (m->control_pressed) { return 0; }
+       
+               if (unaligned) { m->mothurOut("Your sequences need to be aligned when you use the bellerophon method."); m->mothurOutEndLine(); return 1;  }
+               
                int numSeqs = seqs.size();
                
-               if (numSeqs == 0) { mothurOut("Error in reading you sequences."); mothurOutEndLine(); exit(1); }
+               if (numSeqs == 0) { m->mothurOut("Error in reading you sequences."); m->mothurOutEndLine(); exit(1); }
                
                //set default window to 25% of sequence length
                string seq0 = seqs[0]->getAligned();
                if (window == 0) { window = seq0.length() / 4;  }
                else if (window > (seq0.length() / 2)) {  
-                       mothurOut("Your sequence length is = " + toString(seq0.length()) + ". You have selected a window size greater than the length of half your aligned sequence. I will run it with a window size of " + toString((seq0.length() / 2))); mothurOutEndLine();
+                       m->mothurOut("Your sequence length is = " + toString(seq0.length()) + ". You have selected a window size greater than the length of half your aligned sequence. I will run it with a window size of " + toString((seq0.length() / 2))); m->mothurOutEndLine();
                        window = (seq0.length() / 2);
                }
                
                if (increment > (seqs[0]->getAlignLength() - (2*window))) { 
                        if (increment != 10) {
                        
-                               mothurOut("You have selected a increment that is too large. I will use the default."); mothurOutEndLine();
+                               m->mothurOut("You have selected a increment that is too large. I will use the default."); m->mothurOutEndLine();
                                increment = 10;
                                if (increment > (seqs[0]->getAlignLength() - (2*window))) {  increment = 0;  }
                                
@@ -134,25 +151,31 @@ void Bellerophon::getChimeras() {
                int count = 0;
                while (count < iters) {
                                
+                               if (m->control_pressed) { return 0; }
+                               
                                //create 2 vectors of sequences, 1 for left side and one for right side
                                vector<Sequence> left;  vector<Sequence> right;
                                
                                for (int i = 0; i < seqs.size(); i++) {
-//cout << "whole = " << seqs[i].getAligned() << endl;
+                               
+                                       if (m->control_pressed) { return 0; }
+                                       
+//cout << "midpoint = " << midpoint << "\twindow = " << window << endl;
+//cout << "whole = " << seqs[i]->getAligned().length() << endl;
                                        //save left side
                                        string seqLeft = seqs[i]->getAligned().substr(midpoint-window, window);
                                        Sequence tempLeft;
                                        tempLeft.setName(seqs[i]->getName());
                                        tempLeft.setAligned(seqLeft);
                                        left.push_back(tempLeft);
-//cout << "left = " << tempLeft.getAligned() << endl;                  
+//cout << "left = " << tempLeft.getAligned().length() << endl;                 
                                        //save right side
                                        string seqRight = seqs[i]->getAligned().substr(midpoint, window);
                                        Sequence tempRight;
                                        tempRight.setName(seqs[i]->getName());
                                        tempRight.setAligned(seqRight);
                                        right.push_back(tempRight);
-//cout << "right = " << seqRight << endl;      
+//cout << "right = " << seqRight.length() << endl;     
                                }
                                
                                //adjust midpoint by increment
@@ -167,8 +190,14 @@ void Bellerophon::getChimeras() {
                                SparseMatrix* SparseRight = new SparseMatrix();
                                
                                createSparseMatrix(0, left.size(), SparseLeft, left);
+                               
+                               if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
+                               
                                createSparseMatrix(0, right.size(), SparseRight, right);
                                
+                               if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
+                               
+                               left.clear(); right.clear();
                                vector<SeqMap> distMapRight;
                                vector<SeqMap> distMapLeft;
                                
@@ -182,11 +211,13 @@ void Bellerophon::getChimeras() {
                                //cout << "left" << endl << endl;
                                for (MatData currentCell = SparseLeft->begin(); currentCell != SparseLeft->end(); currentCell++) {
                                        distMapLeft[currentCell->row][currentCell->column] = currentCell->dist;
+                                       if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
                                        //cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
                                }
                                //cout << "right" << endl << endl;
                                for (MatData currentCell = SparseRight->begin(); currentCell != SparseRight->end(); currentCell++) {
                                        distMapRight[currentCell->row][currentCell->column] = currentCell->dist;
+                                       if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
                                        //cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
                                }
                                
@@ -220,10 +251,14 @@ void Bellerophon::getChimeras() {
                
                //sort Preferences highest to lowest
                sort(pref.begin(), pref.end(), comparePref);
-
+               
+               for (int i = 0; i < seqs.size(); i++) { delete seqs[i];  }  seqs.clear();
+               
+               return 0;
+               
        }
        catch(exception& e) {
-               errorOut(e, "Bellerophon", "getChimeras");
+               m->errorOut(e, "Bellerophon", "getChimeras");
                exit(1);
        }
 }
@@ -235,6 +270,8 @@ int Bellerophon::createSparseMatrix(int startSeq, int endSeq, SparseMatrix* spar
                for(int i=startSeq; i<endSeq; i++){
                        
                        for(int j=0;j<i;j++){
+                               
+                               if (m->control_pressed) { return 0; }
                        
                                distCalculator->calcDist(s[i], s[j]);
                                float dist = distCalculator->getDist();
@@ -248,12 +285,12 @@ int Bellerophon::createSparseMatrix(int startSeq, int endSeq, SparseMatrix* spar
                return 1;
        }
        catch(exception& e) {
-               errorOut(e, "Bellerophon", "createSparseMatrix");
+               m->errorOut(e, "Bellerophon", "createSparseMatrix");
                exit(1);
        }
 }
 /***************************************************************************************************************/
-void Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right, int mid){
+int Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right, int mid){
        try {
                
                float dme = 0.0;
@@ -277,6 +314,8 @@ void Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right,
                        SeqMap currentRight = right[i];         // same as left but with distances on the right side.
                        
                        for (int j = 0; j < i; j++) {
+                       
+                               if (m->control_pressed) {  return 0; }
                                
                                itL = currentLeft.find(j);
                                itR = currentRight.find(j);
@@ -341,6 +380,8 @@ void Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right,
 //cout << endl << "dme = " << dme << endl;
                //recalculate prefernences based on dme
                for (int i = 0; i < pref.size(); i++) {
+               
+                       if (m->control_pressed) {  return 0; }
 //cout << "unadjusted pref " << i << " = " << pref[i].score[1] << endl;        
                        // gives the actual percentage of the dme this seq adds
                        pref[i].score[1] = pref[i].score[1] / dme;
@@ -357,6 +398,8 @@ void Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right,
                //is this score bigger then the last score
                for (int i = 0; i < pref.size(); i++) {  
                        
+                       if (m->control_pressed) {  return 0; }
+                       
                        //update biggest score
                        if (pref[i].score[1] > pref[i].score[0]) {
                                pref[i].score[0] = pref[i].score[1];
@@ -368,10 +411,12 @@ void Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right,
                        }
                        
                }
+               
+               return 1;
 
        }
        catch(exception& e) {
-               errorOut(e, "Bellerophon", "generatePreferences");
+               m->errorOut(e, "Bellerophon", "generatePreferences");
                exit(1);
        }
 }