test

[mothur.git] / pintail.cpp

diff --git a/pintail.cpp b/pintail.cpp
index 5f575f2bc41721e7fc3b52ec7d6fb68163b43131..6adfdf91ebeaf30847bee05d90f2c9b22d3c5946 100644 (file)
--- a/pintail.cpp
+++ b/pintail.cpp
@@ -18,7 +18,7 @@ inline bool compareQuanMembers(quanMember left, quanMember right){
 } 
 //***************************************************************************************************************
 
-Pintail::Pintail(string filename, string temp) {  fastafile = filename;  templateFile = temp;  }
+Pintail::Pintail(string filename, string temp, string o) {  fastafile = filename;  templateFile = temp; outputDir = o; }
 //***************************************************************************************************************
 
 Pintail::~Pintail() {
@@ -74,7 +74,7 @@ void Pintail::print(ostream& out) {
 }
 
 //***************************************************************************************************************
-void Pintail::getChimeras() {
+int Pintail::getChimeras() {
        try {
                
                //read in query sequences and subject sequences
@@ -85,6 +85,8 @@ void Pintail::getChimeras() {
                
                int numSeqs = querySeqs.size();
                
+               if (unaligned) { mothurOut("Your sequences need to be aligned when you use the pintail method."); mothurOutEndLine(); return 1;  }
+               
                obsDistance.resize(numSeqs);
                expectedDistance.resize(numSeqs);
                seqCoef.resize(numSeqs);
@@ -146,20 +148,20 @@ void Pintail::getChimeras() {
                        mothurOut("Done."); mothurOutEndLine();
                }else {         createProcessesPairs();         }
                
-string o = "closestmatch.eachgap.fasta";
-ofstream out7;
-openOutputFile(o, out7);
+//string o = "closestmatch.eachgap.fasta";
+//ofstream out7;
+//openOutputFile(o, out7);
 
-for (int i = 0; i < bestfit.size(); i++) {
-       out7 << ">" << querySeqs[i]->getName() << "-"<< bestfit[i]->getName() << endl;
-       out7 << bestfit[i]->getAligned() << endl;
-}              
-out7.close();  
+//for (int i = 0; i < bestfit.size(); i++) {
+       //out7 << ">" << querySeqs[i]->getName() << "-"<< bestfit[i]->getName() << endl;
+       //out7 << bestfit[i]->getAligned() << endl;
+//}            
+//out7.close();        
                //find P
                mothurOut("Getting conservation... "); cout.flush();
                if (consfile == "") { 
                        mothurOut("Calculating probability of conservation for your template sequences.  This can take a while...  I will output the frequency of the highest base in each position to a .freq file so that you can input them using the conservation parameter next time you run this command.  Providing the .freq file will improve speed.    "); cout.flush();
-                       probabilityProfile = decalc->calcFreq(templateSeqs, templateFile); 
+                       probabilityProfile = decalc->calcFreq(templateSeqs, outputDir + getSimpleName(templateFile)); 
                        mothurOut("Done."); mothurOutEndLine();
                }else                           {   probabilityProfile = readFreq();                      }
 
@@ -279,13 +281,13 @@ out7.close();
                        string noOutliers, outliers;
                        
                        if ((!filter) && (seqMask == "")) {
-                               noOutliers = getRootName(templateFile) + "pintail.quan";
+                               noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.quan";
                        }else if ((filter) && (seqMask == "")) { 
-                               noOutliers = getRootName(templateFile) + "pintail.filtered.quan";
+                               noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.filtered.quan";
                        }else if ((!filter) && (seqMask != "")) { 
-                               noOutliers = getRootName(templateFile) + "pintail.masked.quan";
+                               noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.masked.quan";
                        }else if ((filter) && (seqMask != "")) { 
-                               noOutliers = getRootName(templateFile) + "pintail.filtered.masked.quan";
+                               noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.filtered.masked.quan";
                        }
 
                        //outliers = getRootName(templateFile) + "pintail.quanYESOUTLIERS";
@@ -382,6 +384,8 @@ out7.close();
                        
                delete distcalculator;
                delete decalc;
+               
+               return 0;
        }
        catch(exception& e) {
                errorOut(e, "Pintail", "getChimeras");
@@ -439,28 +443,9 @@ vector<Sequence*> Pintail::findPairs(int start, int end) {
                vector<Sequence*> seqsMatches;  
                
                for(int i = start; i < end; i++){
-               
-                       float smallest = 10000.0;
-                       Sequence query = *(querySeqs[i]);
-                       Sequence* match;
-                       
-                       for(int j = 0; j < templateSeqs.size(); j++){
-                               
-                               Sequence temp = *(templateSeqs[j]);
-                               
-                               distcalculator->calcDist(query, temp);
-                               float dist = distcalculator->getDist();
-                               
-                               if (dist < smallest) { 
-                                       match = templateSeqs[j];
-                                       smallest = dist;
-                               }
-                       }
-                       
-                       //make copy so trimSeqs doesn't overtrim
-                       Sequence* copy = new Sequence(match->getName(), match->getAligned());
                        
-                       seqsMatches.push_back(copy);
+                       vector<Sequence*> copy = decalc->findClosest(querySeqs[i], templateSeqs, 1);
+                       seqsMatches.push_back(copy[0]);
                }
                
                return seqsMatches;