Merge remote-tracking branch 'mothur/master'

[mothur.git] / shhhercommand.cpp

diff --git a/shhhercommand.cpp b/shhhercommand.cpp
index 28e6615f3dffba8fb88c88cd9ca76e79e66faa06..32891c333ab5229aec50ed05ff5f03b6799d1396 100644 (file)
--- a/shhhercommand.cpp
+++ b/shhhercommand.cpp
@@ -67,6 +67,13 @@ ShhherCommand::ShhherCommand(){
 
 ShhherCommand::ShhherCommand(string option) {
        try {
+        
+#ifdef USE_MPI
+               MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
+               MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
+        
+               if(pid == 0){
+#endif
                abort = false; calledHelp = false;   
                
                //allow user to run help
@@ -314,6 +321,9 @@ ShhherCommand::ShhherCommand(string option) {
                        }
                        
                }
+#ifdef USE_MPI
+               }                               
+#endif
        }
        catch(exception& e) {
                m->errorOut(e, "ShhherCommand", "ShhherCommand");
@@ -328,9 +338,7 @@ int ShhherCommand::execute(){
                
                int tag = 1976;
                MPI_Status status; 
-               MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-               MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
-        
+                       
                if(pid == 0){
 
                        for(int i=1;i<ncpus;i++){
@@ -344,6 +352,22 @@ int ShhherCommand::execute(){
                        getSingleLookUp();      if (m->control_pressed) { return 0; }
                        getJointLookUp();       if (m->control_pressed) { return 0; }
                        
+            vector<string> flowFileVector;
+                       if(flowFilesFileName != "not found"){
+                               string fName;
+                
+                               ifstream flowFilesFile;
+                               m->openInputFile(flowFilesFileName, flowFilesFile);
+                               while(flowFilesFile){
+                                       fName = m->getline(flowFilesFile);
+                                       flowFileVector.push_back(fName);
+                                       m->gobble(flowFilesFile);
+                               }
+                       }
+                       else{
+                               flowFileVector.push_back(flowFileName);
+                       }
+            
                        int numFiles = flowFileVector.size();
 
                        for(int i=1;i<ncpus;i++){
@@ -410,6 +434,8 @@ int ShhherCommand::execute(){
                                
                                if (m->control_pressed) { break; }
                                
+                
+                
                                getOTUData(listFileName);
 
                                m->mothurRemove(distFileName);
@@ -422,6 +448,7 @@ int ShhherCommand::execute(){
 
                                if (m->control_pressed) { break; }
                                
+            
                                for(int i=1;i<ncpus;i++){
                                        MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
                                        MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
@@ -719,6 +746,36 @@ int ShhherCommand::execute(){
        }
 }
 /**************************************************************************************************/
+string ShhherCommand::createNamesFile(){
+       try{
+               
+               vector<string> duplicateNames(numUniques, "");
+               for(int i=0;i<numSeqs;i++){
+                       duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
+               }
+               
+               string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
+               
+               ofstream nameFile;
+               m->openOutputFile(nameFileName, nameFile);
+               
+               for(int i=0;i<numUniques;i++){
+                       
+                       if (m->control_pressed) { break; }
+                       
+            //                 nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
+                       nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
+               }
+               
+               nameFile.close();
+               return  nameFileName;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ShhherCommand", "createNamesFile");
+               exit(1);
+       }
+}
+/**************************************************************************************************/
 
 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
        try{            
@@ -2346,7 +2403,8 @@ int ShhherCommand::getUniques(int numSeqs, int numFlowCells, vector<short>& uniq
                }
                uniqueFlowDataIntI.resize(numFlowCells * numUniques);
                uniqueLengths.resize(numUniques);       
-               
+               uniqueFlowgrams.resize(numFlowCells * numUniques);
+        
                flowDataPrI.resize(numSeqs * numFlowCells, 0);
                for(int i=0;i<flowDataPrI.size();i++)   {       if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]);          }
         
@@ -2397,7 +2455,7 @@ int ShhherCommand::cluster(string filename, string distFileName, string namesFil
                NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
                clusterNameMap->readMap();
                read->read(clusterNameMap);
-               
+        
                ListVector* list = read->getListVector();
                SparseMatrix* matrix = read->getMatrix();
                
@@ -3012,11 +3070,12 @@ void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string filenam
                        
                        if(otuCounts[i] > 0){
                                qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
-                               
+                       
                                int j=4;        //need to get past the first four bases
                                while(qualities[i][j] != -1){
-                                       qualityFile << qualities[i][j] << ' ';
-                                       j++;
+                    qualityFile << qualities[i][j] << ' ';
+                    if (j > qualities[i].size()) { break; }
+                    j++;
                                }
                                qualityFile << endl;
                        }