added warning about merging with something above cutoff to cluster. working on chimeras

[mothur.git] / bellerophon.cpp

diff --git a/bellerophon.cpp b/bellerophon.cpp
index dc021c064c27bf2550da911589cc51fd5f381458..17859dc710309b7535e234aec56791bc7a8ee1b2 100644 (file)
--- a/bellerophon.cpp
+++ b/bellerophon.cpp
@@ -15,9 +15,10 @@
 
 //***************************************************************************************************************
 
-Bellerophon::Bellerophon(string name) {
+Bellerophon::Bellerophon(string name, string o)  {
        try {
                fastafile = name;
+               outputDir = o;
        }
        catch(exception& e) {
                errorOut(e, "Bellerophon", "Bellerophon");
@@ -26,7 +27,7 @@ Bellerophon::Bellerophon(string name) {
 }
 
 //***************************************************************************************************************
-void Bellerophon::print(ostream& out) {
+void Bellerophon::print(ostream& out, ostream& outAcc) {
        try {
                int above1 = 0;
                out << "Name\tScore\tLeft\tRight\t" << endl;
@@ -37,6 +38,7 @@ void Bellerophon::print(ostream& out) {
                        //calc # of seqs with preference above 1.0
                        if (pref[i].score[0] > 1.0) { 
                                above1++; 
+                               outAcc << pref[i].name << endl;
                                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();
                        }
@@ -75,18 +77,22 @@ 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;
                        
                        //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,6 +100,8 @@ void Bellerophon::getChimeras() {
                //read in sequences
                seqs = readSeqs(fastafile);
                
+               if (unaligned) { mothurOut("Your sequences need to be aligned when you use the bellerophon method."); mothurOutEndLine(); return 1;  }
+               
                int numSeqs = seqs.size();
                
                if (numSeqs == 0) { mothurOut("Error in reading you sequences."); mothurOutEndLine(); exit(1); }
@@ -115,7 +123,7 @@ void Bellerophon::getChimeras() {
                                
                        }else{ increment = 0; }
                }
-cout << "increment = " << increment << endl;           
+               
                if (increment == 0) { iters = 1; }
                else { iters = ((seqs[0]->getAlignLength() - (2*window)) / increment); }
                
@@ -138,21 +146,22 @@ cout << "increment = " << increment << endl;
                                vector<Sequence> left;  vector<Sequence> right;
                                
                                for (int i = 0; i < seqs.size(); i++) {
-//cout << "whole = " << seqs[i].getAligned() << endl;
+//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
@@ -193,7 +202,6 @@ cout << "increment = " << increment << endl;
                                delete SparseLeft;
                                delete SparseRight;
                                
-                               
                                //fill preference structure
                                generatePreferences(distMapLeft, distMapRight, midpoint);
                                
@@ -216,14 +224,14 @@ cout << "increment = " << increment << endl;
                        
                        //how much higher or lower is this than expected
                        pref[i].score[0] = pref[i].score[0] / expectedPercent;
-                       
-                       
-                       
+               
                }
                
                //sort Preferences highest to lowest
                sort(pref.begin(), pref.end(), comparePref);
-
+               
+               return 0;
+               
        }
        catch(exception& e) {
                errorOut(e, "Bellerophon", "getChimeras");
@@ -337,7 +345,6 @@ void Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right,
                
                  
                //calculate the dme
-               
                int count0 = 0;
                for (int i = 0; i < pref.size(); i++) {  dme += pref[i].score[1];  if (pref[i].score[1] == 0.0) { count0++; }  }