]> git.donarmstrong.com Git - mothur.git/blobdiff - shhhercommand.cpp
added count file to cluster.classic and cluster.split. modified splitting classes...
[mothur.git] / shhhercommand.cpp
index a09825e4b9a3cd11dbe77b3dd82e703481c650bc..c34f25de509c78b9e1b97b7c303142864e3c7e01 100644 (file)
@@ -47,6 +47,30 @@ string ShhherCommand::getHelpString(){
        }
 }
 //**********************************************************************************************************************
+string ShhherCommand::getOutputFileNameTag(string type, string inputName=""){  
+       try {
+        string outputFileName = "";
+               map<string, vector<string> >::iterator it;
+        
+        //is this a type this command creates
+        it = outputTypes.find(type);
+        if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
+        else {
+            if (type == "fasta")            {   outputFileName =  "shhh.fasta";   }
+            else if (type == "name")    {   outputFileName =  "shhh.names";   }
+            else if (type == "group")        {   outputFileName =  "shhh.groups";   }
+            else if (type == "counts")        {   outputFileName =  "shhh.counts";   }
+            else if (type == "qfile")        {   outputFileName =  "shhh.qual";   }
+            else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
+        }
+        return outputFileName;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ShhherCommand", "getOutputFileNameTag");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
 
 ShhherCommand::ShhherCommand(){        
        try {
@@ -54,8 +78,12 @@ ShhherCommand::ShhherCommand(){
                setParameters();
                
                //initialize outputTypes
-//             vector<string> tempOutNames;
-//             outputTypes["pn.dist"] = tempOutNames;
+               vector<string> tempOutNames;
+               outputTypes["fasta"] = tempOutNames;
+        outputTypes["name"] = tempOutNames;
+        outputTypes["group"] = tempOutNames;
+        outputTypes["counts"] = tempOutNames;
+        outputTypes["qfile"] = tempOutNames;
 
        }
        catch(exception& e) {
@@ -96,9 +124,13 @@ ShhherCommand::ShhherCommand(string option) {
                        }
                        
                        //initialize outputTypes
-                       vector<string> tempOutNames;
-//                     outputTypes["pn.dist"] = tempOutNames;
-                       //                      outputTypes["fasta"] = tempOutNames;
+            vector<string> tempOutNames;
+            outputTypes["fasta"] = tempOutNames;
+            outputTypes["name"] = tempOutNames;
+            outputTypes["group"] = tempOutNames;
+            outputTypes["counts"] = tempOutNames;
+            outputTypes["qfile"] = tempOutNames;
+
                        
                        //if the user changes the input directory command factory will send this info to us in the output parameter 
                        string inputDir = validParameter.validFile(parameters, "inputdir", false);              
@@ -150,14 +182,21 @@ ShhherCommand::ShhherCommand(string option) {
                        else{
                                ofstream temp;
                 
-                string thisoutputDir = m->hasPath(flowFilesFileName); //if user entered a file with a path then preserve it
+                string thisoutputDir = outputDir;
+                if (outputDir == "") {  thisoutputDir =  m->hasPath(flowFilesFileName); } //if user entered a file with a path then preserve it
+                
+                               //we want to rip off .files, and also .flow if its there
+                string fileroot = m->getRootName(m->getSimpleName(flowFilesFileName));
+                if (fileroot[fileroot.length()-1] == '.') {  fileroot = fileroot.substr(0, fileroot.length()-1); } //rip off dot
+                string extension = m->getExtension(fileroot);
+                if (extension == ".flow") { fileroot = m->getRootName(fileroot);  }
+                else { fileroot += "."; } //add back if needed
                 
-                               //flow.files = 9 character offset
-                               compositeFASTAFileName = thisoutputDir + m->getRootName(m->getSimpleName(flowFilesFileName)) + "shhh.fasta";
+                               compositeFASTAFileName = thisoutputDir + fileroot + "shhh.fasta";
                                m->openOutputFile(compositeFASTAFileName, temp);
                                temp.close();
                                
-                               compositeNamesFileName = thisoutputDir + m->getRootName(m->getSimpleName(flowFilesFileName)) + "shhh.names";
+                               compositeNamesFileName = thisoutputDir + fileroot + "shhh.names";
                                m->openOutputFile(compositeNamesFileName, temp);
                                temp.close();
                        }
@@ -763,7 +802,7 @@ string ShhherCommand::createNamesFile(){
                        duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
                }
                
-               string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
+               string nameFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("name");
                
                ofstream nameFile;
                m->openOutputFile(nameFileName, nameFile);
@@ -931,6 +970,8 @@ void ShhherCommand::initPyroCluster(){
     try{
         if (numOTUs < processors) { processors = 1; }
         
+        if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
+        
         dist.assign(numSeqs * numOTUs, 0);
         change.assign(numOTUs, 1);
         centroids.assign(numOTUs, -1);
@@ -940,6 +981,8 @@ void ShhherCommand::initPyroCluster(){
         nSeqsBreaks.assign(processors+1, 0);
         nOTUsBreaks.assign(processors+1, 0);
         
+        if (m->debug) { m->mothurOut("[DEBUG]: made it through the memory allocation.\n"); }
+        
         nSeqsBreaks[0] = 0;
         for(int i=0;i<processors;i++){
             nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
@@ -1603,7 +1646,7 @@ void ShhherCommand::writeQualities(vector<int> otuCounts){
     try {
         string thisOutputDir = outputDir;
         if (outputDir == "") {  thisOutputDir += m->hasPath(flowFileName);  }
-        string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.qual";
+        string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("qfile");
         
         ofstream qualityFile;
         m->openOutputFile(qualityFileName, qualityFile);
@@ -1710,7 +1753,7 @@ void ShhherCommand::writeSequences(vector<int> otuCounts){
     try {
         string thisOutputDir = outputDir;
         if (outputDir == "") {  thisOutputDir += m->hasPath(flowFileName);  }
-        string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.fasta";
+        string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("fasta");
         ofstream fastaFile;
         m->openOutputFile(fastaFileName, fastaFile);
         
@@ -1758,7 +1801,7 @@ void ShhherCommand::writeNames(vector<int> otuCounts){
     try {
         string thisOutputDir = outputDir;
         if (outputDir == "") {  thisOutputDir += m->hasPath(flowFileName);  }
-        string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.names";
+        string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("name");
         ofstream nameFile;
         m->openOutputFile(nameFileName, nameFile);
         
@@ -1796,14 +1839,17 @@ void ShhherCommand::writeGroups(){
     try {
         string thisOutputDir = outputDir;
         if (outputDir == "") {  thisOutputDir += m->hasPath(flowFileName);  }
-        string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
-        string groupFileName = fileRoot + "shhh.groups";
+        string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
+        int pos = fileRoot.find_first_of('.');
+        string fileGroup = fileRoot;
+        if (pos != string::npos) {  fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1)));  }
+        string groupFileName = thisOutputDir + fileRoot + getOutputFileNameTag("group");
         ofstream groupFile;
         m->openOutputFile(groupFileName, groupFile);
         
         for(int i=0;i<numSeqs;i++){
             if (m->control_pressed) { break; }
-            groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
+            groupFile << seqNameVector[i] << '\t' << fileGroup << endl;
         }
         groupFile.close();
         outputNames.push_back(groupFileName);
@@ -1821,7 +1867,7 @@ void ShhherCommand::writeClusters(vector<int> otuCounts){
     try {
         string thisOutputDir = outputDir;
         if (outputDir == "") {  thisOutputDir += m->hasPath(flowFileName);  }
-        string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.counts";
+        string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) +getOutputFileNameTag("counts");
         ofstream otuCountsFile;
         m->openOutputFile(otuCountsFileName, otuCountsFile);
         
@@ -1927,12 +1973,14 @@ int ShhherCommand::createProcesses(vector<string> filenames){
                
                //divide the groups between the processors
                vector<linePair> lines;
+        vector<int> numFilesToComplete;
                int numFilesPerProcessor = filenames.size() / processors;
                for (int i = 0; i < processors; i++) {
                        int startIndex =  i * numFilesPerProcessor;
                        int endIndex = (i+1) * numFilesPerProcessor;
                        if(i == (processors - 1)){      endIndex = filenames.size();    }
                        lines.push_back(linePair(startIndex, endIndex));
+            numFilesToComplete.push_back((endIndex-startIndex));
                }
                
         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)         
@@ -1946,6 +1994,14 @@ int ShhherCommand::createProcesses(vector<string> filenames){
                                process++;
                        }else if (pid == 0){
                                num = driver(filenames, compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName  + toString(getpid()) + ".temp", lines[process].start, lines[process].end);
+                
+                //pass numSeqs to parent
+                               ofstream out;
+                               string tempFile = compositeFASTAFileName + toString(getpid()) + ".num.temp";
+                               m->openOutputFile(tempFile, out);
+                               out << num << endl;
+                               out.close();
+                
                                exit(0);
                        }else { 
                                m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
@@ -2008,6 +2064,18 @@ int ShhherCommand::createProcesses(vector<string> filenames){
         #endif
         
         for (int i=0;i<processIDS.size();i++) { 
+            ifstream in;
+                       string tempFile =  compositeFASTAFileName + toString(processIDS[i]) + ".num.temp";
+                       m->openInputFile(tempFile, in);
+                       if (!in.eof()) { 
+                int tempNum = 0; 
+                in >> tempNum; 
+                if (tempNum != numFilesToComplete[i+1]) {
+                    m->mothurOut("[ERROR]: main process expected " + toString(processIDS[i]) + " to complete " + toString(numFilesToComplete[i+1]) + " files, and it only reported completing " + toString(tempNum) + ". This will cause file mismatches.  The flow files may be too large to process with multiple processors. \n");
+                }
+            }
+                       in.close(); m->mothurRemove(tempFile);
+            
             if (compositeFASTAFileName != "") {
                 m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName);
                 m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName);
@@ -2046,18 +2114,24 @@ vector<string> ShhherCommand::parseFlowFiles(string filename){
             out << thisNumFLows << endl;
             files.push_back(outputFileName);
             
+            int numLinesWrote = 0;
             for (int i = 0; i < largeSize; i++) {
                 if (in.eof()) { break; }
-                string line = m->getline(in);
+                string line = m->getline(in); m->gobble(in);
                 out << line << endl;
+                numLinesWrote++;
             }
             out.close();
+            
+            if (numLinesWrote == 0) {  m->mothurRemove(outputFileName); files.pop_back();  }
             count++;
         }
         in.close();
         
         if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); }  files.clear(); }
         
+        m->mothurOut("\nDivided " + filename + " into " + toString(files.size()) + " files.\n\n"); 
+        
         return files;
     }
        catch(exception& e) {
@@ -2070,6 +2144,8 @@ vector<string> ShhherCommand::parseFlowFiles(string filename){
 int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){
     try {
         
+        int numCompleted = 0;
+        
         for(int i=start;i<end;i++){
                        
                        if (m->control_pressed) { break; }
@@ -2100,6 +2176,7 @@ int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFil
                 vector<int> uniqueLengths;
                 int numFlowCells;
                 
+                if (m->debug) { m->mothurOut("[DEBUG]: About to read flowgrams.\n"); }
                 int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells);
                 
                 if (m->control_pressed) { break; }
@@ -2156,6 +2233,8 @@ int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFil
                 vector<int> nSeqsBreaks;
                 vector<int> nOTUsBreaks;
                 
+                if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
+                
                 dist.assign(numSeqs * numOTUs, 0);
                 change.assign(numOTUs, 1);
                 centroids.assign(numOTUs, -1);
@@ -2169,6 +2248,8 @@ int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFil
                 nSeqsBreaks[1] = numSeqs;
                 nOTUsBreaks[1] = numOTUs;
                 
+                if (m->debug) { m->mothurOut("[DEBUG]: done allocating memory, about to denoise.\n"); }
+                
                 if (m->control_pressed) { break; }
                 
                 double maxDelta = 0;
@@ -2228,7 +2309,7 @@ int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFil
                 if (m->control_pressed) { break; }
                 
                 vector<int> otuCounts(numOTUs, 0);
-                for(int i=0;i<numSeqs;i++)     {       otuCounts[otuData[i]]++;        }
+                for(int j=0;j<numSeqs;j++)     {       otuCounts[otuData[j]]++;        }
                 
                 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);        
                 
@@ -2237,43 +2318,47 @@ int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFil
                 if ((large) && (g == 0)) {  flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; }
                 string thisOutputDir = outputDir;
                 if (outputDir == "") {  thisOutputDir = m->hasPath(flowFileName);  }
-                string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.qual";
-                string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.fasta";
-                string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.names";
-                string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.counts";
-                string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
-                string groupFileName = fileRoot + "shhh.groups";
+                string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("qfile");
+                string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("fasta");
+                string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("name");
+                string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("counts");
+                string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
+                int pos = fileRoot.find_first_of('.');
+                string fileGroup = fileRoot;
+                if (pos != string::npos) {  fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1)));  }
+                string groupFileName = thisOutputDir + fileRoot + getOutputFileNameTag("group");
 
                 
                 writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; }
                 writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; }
                 writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU);                             if (m->control_pressed) { break; }
                 writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI);                 if (m->control_pressed) { break; }
-                writeGroups(groupFileName, fileRoot, numSeqs, seqNameVector);                                          if (m->control_pressed) { break; }
+                writeGroups(groupFileName, fileGroup, numSeqs, seqNameVector);                                         if (m->control_pressed) { break; }
                 
                 if (large) {
                     if (g > 0) {
-                        m->appendFiles(qualityFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.qual"));
+                        m->appendFiles(qualityFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("qfile")));
                         m->mothurRemove(qualityFileName);
-                        m->appendFiles(fastaFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.fasta"));
+                        m->appendFiles(fastaFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("fasta")));
                         m->mothurRemove(fastaFileName);
-                        m->appendFiles(nameFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.names"));
+                        m->appendFiles(nameFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("name")));
                         m->mothurRemove(nameFileName);
-                        m->appendFiles(otuCountsFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.counts"));
+                        m->appendFiles(otuCountsFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("counts")));
                         m->mothurRemove(otuCountsFileName);
-                        m->appendFiles(groupFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.groups"));
+                        m->appendFiles(groupFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("group")));
                         m->mothurRemove(groupFileName);
-                        m->mothurRemove(theseFlowFileNames[g]);
                     }
+                    m->mothurRemove(theseFlowFileNames[g]);
                 }
                        }
             
+            numCompleted++;
                        m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
                }
                
         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
         
-        return 0;
+        return numCompleted;
         
     }catch(exception& e) {
             m->errorOut(e, "ShhherCommand", "driver");
@@ -2298,17 +2383,21 @@ int ShhherCommand::getFlowData(string filename, vector<string>& thisSeqNameVecto
                thisNameMap.clear();
                
                flowFile >> numFlowCells;
+        if (m->debug) { m->mothurOut("[DEBUG]: numFlowCells = " + toString(numFlowCells) + ".\n"); }
                int index = 0;//pcluster
                while(!flowFile.eof()){
                        
                        if (m->control_pressed) { break; }
                        
                        flowFile >> seqName >> currentNumFlowCells;
+            
                        thisLengths.push_back(currentNumFlowCells);
            
                        thisSeqNameVector.push_back(seqName);
                        thisNameMap[seqName] = index++;//pcluster
-
+            
+            if (m->debug) { m->mothurOut("[DEBUG]: seqName = " + seqName + " length = " + toString(currentNumFlowCells) + " index = " + toString(index) + "\n"); }
+            
                        for(int i=0;i<numFlowCells;i++){
                                flowFile >> intensity;
                                if(intensity > 9.99)    {       intensity = 9.99;       }
@@ -2543,7 +2632,7 @@ int ShhherCommand::cluster(string filename, string distFileName, string namesFil
                read->read(clusterNameMap);
         
                ListVector* list = read->getListVector();
-               SparseMatrix* matrix = read->getMatrix();
+               SparseDistanceMatrix* matrix = read->getDMatrix();
                
                delete read; 
                delete clusterNameMap; 
@@ -2596,6 +2685,8 @@ int ShhherCommand::getOTUData(int numSeqs, string fileName,  vector<int>& otuDat
                
                listFile >> label >> numOTUs;
         
+        if (m->debug) { m->mothurOut("[DEBUG]: Getting OTU Data...\n"); }
+        
                otuData.assign(numSeqs, 0);
                cumNumSeqs.assign(numOTUs, 0);
                nSeqsPerOTU.assign(numOTUs, 0);
@@ -2610,6 +2701,7 @@ int ShhherCommand::getOTUData(int numSeqs, string fileName,  vector<int>& otuDat
                for(int i=0;i<numOTUs;i++){
                        
                        if (m->control_pressed) { break; }
+            if (m->debug) { m->mothurOut("[DEBUG]: processing OTU " + toString(i) + ".\n"); }
             
                        listFile >> singleOTU;
                        
@@ -3319,7 +3411,9 @@ void ShhherCommand::writeClusters(string otuCountsFileName, int numOTUs, int num
                                                        //otuCountsFile << base;
                                                }
                                        }
-                                       otuCountsFile << newSeq.substr(4) << endl;
+                                       
+                    if (newSeq.length() >= 4) {  otuCountsFile << newSeq.substr(4) << endl;  }
+                    else {  otuCountsFile << "NNNN" << endl;  }
                                }
                                otuCountsFile << endl;
                        }