]> git.donarmstrong.com Git - mothur.git/blobdiff - seqsummarycommand.cpp
fixed bug with shhh.flow from file path name in write functions, added "smart" featur...
[mothur.git] / seqsummarycommand.cpp
index 911e65fc9f8696834448ff9b3451b53e43f496b6..e8f73ca10875f8321f624adfb62ee02f7baf2770 100644 (file)
@@ -8,13 +8,14 @@
  */
 
 #include "seqsummarycommand.h"
-#include "sequence.hpp"
+
 
 //**********************************************************************************************************************
 vector<string> SeqSummaryCommand::setParameters(){     
        try {
                CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
                CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
+               CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
                CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
                CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
                
@@ -36,7 +37,7 @@ string SeqSummaryCommand::getHelpString(){
                helpString += "The name parameter allows you to enter a name file associated with your fasta file. \n";
                helpString += "The summary.seqs command should be in the following format: \n";
                helpString += "summary.seqs(fasta=yourFastaFile, processors=2) \n";
-               helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n\n";     
+               helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n";       
                return helpString;
        }
        catch(exception& e) {
@@ -66,6 +67,7 @@ SeqSummaryCommand::SeqSummaryCommand(string option)  {
                
                //allow user to run help
                if(option == "help") { help(); abort = true; calledHelp = true; }
+               else if(option == "citation") { citation(); abort = true; calledHelp = true;}
                
                else {
                        vector<string> myArray = setParameters();
@@ -114,11 +116,12 @@ SeqSummaryCommand::SeqSummaryCommand(string option)  {
                                fastafile = m->getFastaFile(); 
                                if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
                                else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
-                       }       
+                       }else { m->setFastaFile(fastafile); }   
                        
                        namefile = validParameter.validFile(parameters, "name", true);
                        if (namefile == "not open") { namefile = ""; abort = true; }
                        else if (namefile == "not found") { namefile = "";  }   
+                       else { m->setNameFile(namefile); }
                        
                        //if the user changes the output directory command factory will send this info to us in the output parameter 
                        outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  
@@ -128,9 +131,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) {
@@ -167,7 +174,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;
@@ -278,31 +285,30 @@ int SeqSummaryCommand::execute(){
                                
                                MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 #else
-                       vector<unsigned long int> positions = m->divideFile(fastafile, processors);
-                               
-                       for (int i = 0; i < (positions.size()-1); i++) {
-                               lines.push_back(new linePair(positions[i], positions[(i+1)]));
-                       }       
-
-               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-                               if(processors == 1){
-                                       numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, summaryFile, lines[0]);
-                               }else{
-                                       numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, summaryFile); 
-                                       
-                                       rename((summaryFile + toString(processIDS[0]) + ".temp").c_str(), summaryFile.c_str());
-                                       //append files
-                                       for(int i=1;i<processors;i++){
-                                               m->appendFiles((summaryFile + toString(processIDS[i]) + ".temp"), summaryFile);
-                                               remove((summaryFile + toString(processIDS[i]) + ".temp").c_str());
-                                       }
+                       vector<unsigned long long> positions; 
+                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+                               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); 
+               
+                               //figure out how many sequences you have to process
+                               int numSeqsPerProcessor = numSeqs / processors;
+                               for (int i = 0; i < processors; i++) {
+                                       int startIndex =  i * numSeqsPerProcessor;
+                                       if(i == (processors - 1)){      numSeqsPerProcessor = numSeqs - i * numSeqsPerProcessor;        }
+                                       lines.push_back(new linePair(positions[startIndex], numSeqsPerProcessor));
                                }
-                               
-                               if (m->control_pressed) {  return 0; }
-               #else
+                       #endif
+                       
+
+                       if(processors == 1){
                                numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, summaryFile, lines[0]);
-                               if (m->control_pressed) {  return 0; }
-               #endif
+                       }else{
+                               numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, summaryFile); 
+                       }
+                       
+                       if (m->control_pressed) {  return 0; }
 #endif
                        
                #ifdef USE_MPI
@@ -315,6 +321,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);
@@ -327,21 +346,23 @@ int SeqSummaryCommand::execute(){
                if (startPosition[0] == -1) {  startPosition[0] = 0;    }
                if (endPosition[0] == -1)       {  endPosition[0] = 0;          }
                
-               if (m->control_pressed) {  remove(summaryFile.c_str()); return 0; }
+               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(); }
                
-               if (m->control_pressed) {  remove(summaryFile.c_str()); return 0; }
+               if (m->control_pressed) {  m->mothurRemove(summaryFile); return 0; }
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Name: "); m->mothurOutEndLine();
@@ -413,7 +434,7 @@ int SeqSummaryCommand::driverCreateSummary(vector<int>& startPosition, vector<in
                        }
                        
                        #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-                               unsigned long int pos = in.tellg();
+                               unsigned long long pos = in.tellg();
                                if ((pos == -1) || (pos >= filePos->end)) { break; }
                        #else
                                if (in.eof()) { break; }
@@ -436,7 +457,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;
@@ -504,11 +525,12 @@ int SeqSummaryCommand::MPICreateSummary(int start, int num, vector<int>& startPo
 /**************************************************************************************************/
 int SeqSummaryCommand::createProcessesCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, string sumFile) {
        try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-               int process = 0;
+               int process = 1;
                int num = 0;
                processIDS.clear();
                
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+               
                //loop through and create all the processes you want
                while (process != processors) {
                        int pid = fork();
@@ -542,8 +564,11 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector<int>& startPosition,
                        }
                }
                
+               //do your part
+               num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, sumFile, lines[0]);
+
                //force parent to wait until all the processes are done
-               for (int i=0;i<processors;i++) { 
+               for (int i=0;i<processIDS.size();i++) { 
                        int temp = processIDS[i];
                        wait(&temp);
                }
@@ -564,11 +589,55 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector<int>& startPosition,
                        for (int k = 0; k < tempNum; k++)                       {               in >> temp; longHomoPolymer.push_back(temp);    }               m->gobble(in);
                                
                        in.close();
-                       remove(tempFilename.c_str());
+                       m->mothurRemove(tempFilename);
+                       
+                       m->appendFiles((sumFile + toString(processIDS[i]) + ".temp"), sumFile);
+                       m->mothurRemove((sumFile + toString(processIDS[i]) + ".temp"));
                }
                
-               return num;
+#else
+               //////////////////////////////////////////////////////////////////////////////////////////////////////
+               //Windows version shared memory, so be careful when passing variables through the seqSumData struct. 
+               //Above fork() will clone, so memory is separate, but that's not the case with windows, 
+               //Taking advantage of shared memory to allow both threads to add info to vectors.
+               //////////////////////////////////////////////////////////////////////////////////////////////////////
+               
+               vector<seqSumData*> pDataArray; 
+               DWORD   dwThreadIdArray[processors];
+               HANDLE  hThreadArray[processors]; 
+               
+               //Create processor worker threads.
+               for( int i=0; i<processors; i++ ){
+                       
+                       //cout << i << '\t' << lines[i]->start << '\t' << lines[i]->end << endl;
+                       // 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);
+                       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]);   
+               }
+                       
+               //Wait until all threads have terminated.
+               WaitForMultipleObjects(processors, hThreadArray, TRUE, INFINITE);
+               
+               //Close all thread handles and free memory allocations.
+               for(int i=0; i < pDataArray.size(); i++){
+                       num += pDataArray[i]->count;
+                       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);
+                       m->mothurRemove((sumFile + toString(processIDS[i]) + ".temp"));
+               }
 #endif         
+               return num;
        }
        catch(exception& e) {
                m->errorOut(e, "SeqSummaryCommand", "createProcessesCreateSummary");