]> git.donarmstrong.com Git - mothur.git/blobdiff - seqsummarycommand.cpp
working on adding count file to mgcluster
[mothur.git] / seqsummarycommand.cpp
index ee8835386e1937bd8bc6d51ab2b0779f05f6a679..c328a041267d01e21ed2e99780f18de486571685 100644 (file)
@@ -45,6 +45,26 @@ string SeqSummaryCommand::getHelpString(){
                exit(1);
        }
 }
+//**********************************************************************************************************************
+string SeqSummaryCommand::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 == "summary")            {   outputFileName =  "summary";   }
+            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, "SeqSummaryCommand", "getOutputFileNameTag");
+               exit(1);
+       }
+}
 
 //**********************************************************************************************************************
 SeqSummaryCommand::SeqSummaryCommand(){        
@@ -131,9 +151,13 @@ SeqSummaryCommand::SeqSummaryCommand(string option)  {
                        
                        string temp = validParameter.validFile(parameters, "processors", false);        if (temp == "not found"){       temp = m->getProcessors();      }
                        m->setProcessors(temp);
-                       convert(temp, processors);
-
-
+                       m->mothurConvert(temp, processors);
+                       
+                       if (namefile == "") {
+                               vector<string> files; files.push_back(fastafile);
+                               parser.getNameFile(files);
+                       }
+                       
                }
        }
        catch(exception& e) {
@@ -151,7 +175,7 @@ int SeqSummaryCommand::execute(){
                //set current fasta to fastafile
                m->setFastaFile(fastafile);
                
-               string summaryFile = outputDir + m->getSimpleName(fastafile) + ".summary";
+               string summaryFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + getOutputFileNameTag("summary");
                                
                int numSeqs = 0;
                
@@ -170,7 +194,7 @@ int SeqSummaryCommand::execute(){
                                int tag = 2001;
                                int startTag = 1; int endTag = 2; int lengthTag = 3; int baseTag = 4; int lhomoTag = 5;
                                int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
-                               vector<unsigned long int> MPIPos;
+                               vector<unsigned long long> MPIPos;
                                
                                MPI_Status status; 
                                MPI_Status statusOut;
@@ -281,12 +305,13 @@ int SeqSummaryCommand::execute(){
                                
                                MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 #else
-                       vector<unsigned long int> positions; 
-                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+                       vector<unsigned long long> positions; 
+                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
                                positions = m->divideFile(fastafile, processors);
                                for (int i = 0; i < (positions.size()-1); i++) {        lines.push_back(new linePair(positions[i], positions[(i+1)]));  }
                        #else
                                positions = m->setFilePosFasta(fastafile, numSeqs); 
+                if (positions.size() < processors) { processors = positions.size(); }
                
                                //figure out how many sequences you have to process
                                int numSeqsPerProcessor = numSeqs / processors;
@@ -317,6 +342,19 @@ int SeqSummaryCommand::execute(){
                sort(ambigBases.begin(), ambigBases.end());
                sort(longHomoPolymer.begin(), longHomoPolymer.end());
                int size = startPosition.size();
+               
+               //find means
+               float meanStartPosition, meanEndPosition, meanSeqLength, meanAmbigBases, meanLongHomoPolymer;
+               meanStartPosition = 0; meanEndPosition = 0; meanSeqLength = 0; meanAmbigBases = 0; meanLongHomoPolymer = 0;
+               for (int i = 0; i < size; i++) {
+                       meanStartPosition += startPosition[i];
+                       meanEndPosition += endPosition[i];
+                       meanSeqLength += seqLength[i];
+                       meanAmbigBases += ambigBases[i];
+                       meanLongHomoPolymer += longHomoPolymer[i];
+               }
+               //this is an int divide so the remainder is lost
+               meanStartPosition /= (float) size; meanEndPosition /= (float) size; meanLongHomoPolymer /= (float) size; meanSeqLength /= (float) size; meanAmbigBases /= (float) size;
                                
                int ptile0_25   = int(size * 0.025);
                int ptile25             = int(size * 0.250);
@@ -332,14 +370,16 @@ int SeqSummaryCommand::execute(){
                if (m->control_pressed) {  m->mothurRemove(summaryFile); return 0; }
                
                m->mothurOutEndLine();
-               m->mothurOut("\t\tStart\tEnd\tNBases\tAmbigs\tPolymer"); m->mothurOutEndLine();
-               m->mothurOut("Minimum:\t" + toString(startPosition[0]) + "\t" + toString(endPosition[0]) + "\t" + toString(seqLength[0]) + "\t" + toString(ambigBases[0]) + "\t" + toString(longHomoPolymer[0])); m->mothurOutEndLine();
-               m->mothurOut("2.5%-tile:\t" + toString(startPosition[ptile0_25]) + "\t" + toString(endPosition[ptile0_25]) + "\t" + toString(seqLength[ptile0_25]) + "\t" + toString(ambigBases[ptile0_25]) + "\t"+ toString(longHomoPolymer[ptile0_25])); m->mothurOutEndLine();
-               m->mothurOut("25%-tile:\t" + toString(startPosition[ptile25]) + "\t" + toString(endPosition[ptile25]) + "\t" + toString(seqLength[ptile25]) + "\t" + toString(ambigBases[ptile25]) + "\t" + toString(longHomoPolymer[ptile25])); m->mothurOutEndLine();
-               m->mothurOut("Median: \t" + toString(startPosition[ptile50]) + "\t" + toString(endPosition[ptile50]) + "\t" + toString(seqLength[ptile50]) + "\t" + toString(ambigBases[ptile50]) + "\t" + toString(longHomoPolymer[ptile50])); m->mothurOutEndLine();
-               m->mothurOut("75%-tile:\t" + toString(startPosition[ptile75]) + "\t" + toString(endPosition[ptile75]) + "\t" + toString(seqLength[ptile75]) + "\t" + toString(ambigBases[ptile75]) + "\t" + toString(longHomoPolymer[ptile75])); m->mothurOutEndLine();
-               m->mothurOut("97.5%-tile:\t" + toString(startPosition[ptile97_5]) + "\t" + toString(endPosition[ptile97_5]) + "\t" + toString(seqLength[ptile97_5]) + "\t" + toString(ambigBases[ptile97_5]) + "\t" + toString(longHomoPolymer[ptile97_5])); m->mothurOutEndLine();
-               m->mothurOut("Maximum:\t" + toString(startPosition[ptile100]) + "\t" + toString(endPosition[ptile100]) + "\t" + toString(seqLength[ptile100]) + "\t" + toString(ambigBases[ptile100]) + "\t" + toString(longHomoPolymer[ptile100])); m->mothurOutEndLine();
+               m->mothurOut("\t\tStart\tEnd\tNBases\tAmbigs\tPolymer\tNumSeqs"); m->mothurOutEndLine();
+               m->mothurOut("Minimum:\t" + toString(startPosition[0]) + "\t" + toString(endPosition[0]) + "\t" + toString(seqLength[0]) + "\t" + toString(ambigBases[0]) + "\t" + toString(longHomoPolymer[0]) + "\t" + toString(1)); m->mothurOutEndLine();
+               m->mothurOut("2.5%-tile:\t" + toString(startPosition[ptile0_25]) + "\t" + toString(endPosition[ptile0_25]) + "\t" + toString(seqLength[ptile0_25]) + "\t" + toString(ambigBases[ptile0_25]) + "\t"+ toString(longHomoPolymer[ptile0_25]) + "\t" + toString(ptile0_25+1)); m->mothurOutEndLine();
+               m->mothurOut("25%-tile:\t" + toString(startPosition[ptile25]) + "\t" + toString(endPosition[ptile25]) + "\t" + toString(seqLength[ptile25]) + "\t" + toString(ambigBases[ptile25]) + "\t" + toString(longHomoPolymer[ptile25]) + "\t" + toString(ptile25+1)); m->mothurOutEndLine();
+               m->mothurOut("Median: \t" + toString(startPosition[ptile50]) + "\t" + toString(endPosition[ptile50]) + "\t" + toString(seqLength[ptile50]) + "\t" + toString(ambigBases[ptile50]) + "\t" + toString(longHomoPolymer[ptile50]) + "\t" + toString(ptile50+1)); m->mothurOutEndLine();
+               m->mothurOut("75%-tile:\t" + toString(startPosition[ptile75]) + "\t" + toString(endPosition[ptile75]) + "\t" + toString(seqLength[ptile75]) + "\t" + toString(ambigBases[ptile75]) + "\t" + toString(longHomoPolymer[ptile75]) + "\t" + toString(ptile75+1)); m->mothurOutEndLine();
+               m->mothurOut("97.5%-tile:\t" + toString(startPosition[ptile97_5]) + "\t" + toString(endPosition[ptile97_5]) + "\t" + toString(seqLength[ptile97_5]) + "\t" + toString(ambigBases[ptile97_5]) + "\t" + toString(longHomoPolymer[ptile97_5]) + "\t" + toString(ptile97_5+1)); m->mothurOutEndLine();
+               m->mothurOut("Maximum:\t" + toString(startPosition[ptile100]) + "\t" + toString(endPosition[ptile100]) + "\t" + toString(seqLength[ptile100]) + "\t" + toString(ambigBases[ptile100]) + "\t" + toString(longHomoPolymer[ptile100]) + "\t" + toString(ptile100+1)); m->mothurOutEndLine();
+               m->mothurOut("Mean:\t" + toString(meanStartPosition) + "\t" + toString(meanEndPosition) + "\t" + toString(meanSeqLength) + "\t" + toString(meanAmbigBases) + "\t" + toString(meanLongHomoPolymer)); m->mothurOutEndLine();
+
                if (namefile == "") {  m->mothurOut("# of Seqs:\t" + toString(numSeqs)); m->mothurOutEndLine(); }
                else { m->mothurOut("# of unique seqs:\t" + toString(numSeqs)); m->mothurOutEndLine(); m->mothurOut("total # of seqs:\t" + toString(startPosition.size())); m->mothurOutEndLine(); }
                
@@ -394,7 +434,7 @@ int SeqSummaryCommand::driverCreateSummary(vector<int>& startPosition, vector<in
                                        //make sure this sequence is in the namefile, else error 
                                        map<string, int>::iterator it = nameMap.find(current.getName());
                                        
-                                       if (it == nameMap.end()) { m->mothurOut("[ERROR]: " + current.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); m->control_pressed = true; }
+                                       if (it == nameMap.end()) { m->mothurOut("[ERROR]: '" + current.getName() + "' is not in your namefile, please correct."); m->mothurOutEndLine(); m->control_pressed = true; }
                                        else { num = it->second; }
                                }
                                
@@ -414,8 +454,8 @@ int SeqSummaryCommand::driverCreateSummary(vector<int>& startPosition, vector<in
                                outSummary << current.getLongHomoPolymer() << '\t' << num << endl;
                        }
                        
-                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-                               unsigned long int pos = in.tellg();
+                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+                               unsigned long long pos = in.tellg();
                                if ((pos == -1) || (pos >= filePos->end)) { break; }
                        #else
                                if (in.eof()) { break; }
@@ -438,7 +478,7 @@ int SeqSummaryCommand::driverCreateSummary(vector<int>& startPosition, vector<in
 }
 #ifdef USE_MPI
 /**************************************************************************************/
-int SeqSummaryCommand::MPICreateSummary(int start, int num, vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, MPI_File& inMPI, MPI_File& outMPI, vector<unsigned long int>& MPIPos) {       
+int SeqSummaryCommand::MPICreateSummary(int start, int num, vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, MPI_File& inMPI, MPI_File& outMPI, vector<unsigned long long>& MPIPos) {      
        try {
                
                int pid;
@@ -510,7 +550,7 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector<int>& startPosition,
                int num = 0;
                processIDS.clear();
                
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
                
                //loop through and create all the processes you want
                while (process != processors) {
@@ -584,34 +624,42 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector<int>& startPosition,
                //////////////////////////////////////////////////////////////////////////////////////////////////////
                
                vector<seqSumData*> pDataArray; 
-               DWORD   dwThreadIdArray[processors];
-               HANDLE  hThreadArray[processors]; 
+               DWORD   dwThreadIdArray[processors-1];
+               HANDLE  hThreadArray[processors-1]; 
                
                //Create processor worker threads.
-               for( int i=0; i<processors; i++ ){
-                       
-                       //cout << i << '\t' << lines[i]->start << '\t' << lines[i]->end << endl;
+               for( int i=0; i<processors-1; i++ ){
+            
+            string extension = "";
+            if (i != 0) { extension = toString(i) + ".temp"; processIDS.push_back(i); }
                        // Allocate memory for thread data.
-                       seqSumData* tempSum = new seqSumData(&startPosition, &endPosition, &seqLength, &ambigBases, &longHomoPolymer, filename, (sumFile + toString(i) + ".temp"), m, lines[i]->start, lines[i]->end, namefile, nameMap);
+                       seqSumData* tempSum = new seqSumData(filename, (sumFile+extension), m, lines[i]->start, lines[i]->end, namefile, nameMap);
                        pDataArray.push_back(tempSum);
-                       processIDS.push_back(i);
-                               
+                       
                        //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
                        //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
                        hThreadArray[i] = CreateThread(NULL, 0, MySeqSumThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);   
                }
-                       
+               
+        //do your part
+               num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, (sumFile+toString(processors-1)+".temp"), lines[processors-1]);
+        processIDS.push_back(processors-1);
+
                //Wait until all threads have terminated.
-               WaitForMultipleObjects(processors, hThreadArray, TRUE, INFINITE);
+               WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
                
                //Close all thread handles and free memory allocations.
                for(int i=0; i < pDataArray.size(); i++){
                        num += pDataArray[i]->count;
+            for (int k = 0; k < pDataArray[i]->startPosition.size(); k++) {    startPosition.push_back(pDataArray[i]->startPosition[k]);       }
+                       for (int k = 0; k < pDataArray[i]->endPosition.size(); k++) {   endPosition.push_back(pDataArray[i]->endPosition[k]);       }
+            for (int k = 0; k < pDataArray[i]->seqLength.size(); k++) {        seqLength.push_back(pDataArray[i]->seqLength[k]);       }
+            for (int k = 0; k < pDataArray[i]->ambigBases.size(); k++) {       ambigBases.push_back(pDataArray[i]->ambigBases[k]);       }
+            for (int k = 0; k < pDataArray[i]->longHomoPolymer.size(); k++) {  longHomoPolymer.push_back(pDataArray[i]->longHomoPolymer[k]);       }
                        CloseHandle(hThreadArray[i]);
                        delete pDataArray[i];
                }
-               
-               //rename((sumFile + toString(processIDS[0]) + ".temp").c_str(), sumFile.c_str());
+    
                //append files
                for(int i=0;i<processIDS.size();i++){
                        m->appendFiles((sumFile + toString(processIDS[i]) + ".temp"), sumFile);