X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=seqsummarycommand.cpp;h=620de95c333cdee1360157f7aa7068f3ce1f1f62;hb=a150db3c2bfd7b76420048a1e94ebe397f2c6045;hp=7362c1b6cfe9e8b11c03ff4c69d3e32b83531942;hpb=284fd95c611ccc3b1a7875c4dacfca06d1f50ed6;p=mothur.git diff --git a/seqsummarycommand.cpp b/seqsummarycommand.cpp index 7362c1b..620de95 100644 --- a/seqsummarycommand.cpp +++ b/seqsummarycommand.cpp @@ -10,19 +10,67 @@ #include "seqsummarycommand.h" #include "sequence.hpp" +//********************************************************************************************************************** +vector 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); + + vector myArray; + for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } + return myArray; + } + catch(exception& e) { + m->errorOut(e, "SeqSummaryCommand", "setParameters"); + exit(1); + } +} +//********************************************************************************************************************** +string SeqSummaryCommand::getHelpString(){ + try { + string helpString = ""; + helpString += "The summary.seqs command reads a fastafile and summarizes the sequences.\n"; + helpString += "The summary.seqs command parameters are fasta, name and processors, fasta is required, unless you have a valid current fasta file.\n"; + 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"; + return helpString; + } + catch(exception& e) { + m->errorOut(e, "SeqSummaryCommand", "getHelpString"); + exit(1); + } +} + +//********************************************************************************************************************** +SeqSummaryCommand::SeqSummaryCommand(){ + try { + abort = true; calledHelp = true; + setParameters(); + vector tempOutNames; + outputTypes["summary"] = tempOutNames; + } + catch(exception& e) { + m->errorOut(e, "SeqSummaryCommand", "SeqSummaryCommand"); + exit(1); + } +} //*************************************************************************************************************** SeqSummaryCommand::SeqSummaryCommand(string option) { try { - abort = false; + abort = false; calledHelp = false; //allow user to run help - if(option == "help") { help(); abort = true; } + if(option == "help") { help(); abort = true; calledHelp = true; } + else if(option == "citation") { citation(); abort = true; calledHelp = true;} else { - //valid paramters for this command - string Array[] = {"fasta","processors","outputdir","inputdir"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + vector myArray = setParameters(); OptionParser parser(option); map parameters = parser.getParameters(); @@ -43,25 +91,46 @@ SeqSummaryCommand::SeqSummaryCommand(string option) { it = parameters.find("fasta"); //user has given a template file if(it != parameters.end()){ - path = hasPath(it->second); + path = m->hasPath(it->second); //if the user has not given a path then, add inputdir. else leave path alone. if (path == "") { parameters["fasta"] = inputDir + it->second; } } + + it = parameters.find("name"); + //user has given a template file + if(it != parameters.end()){ + path = m->hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["name"] = inputDir + it->second; } + } } + //initialize outputTypes + vector tempOutNames; + outputTypes["summary"] = tempOutNames; + //check for required parameters fastafile = validParameter.validFile(parameters, "fasta", true); if (fastafile == "not open") { abort = true; } - else if (fastafile == "not found") { fastafile = ""; m->mothurOut("fasta is a required parameter for the summary.seqs command."); m->mothurOutEndLine(); abort = true; } + else if (fastafile == "not found") { + 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; } + } + + namefile = validParameter.validFile(parameters, "name", true); + if (namefile == "not open") { namefile = ""; abort = true; } + else if (namefile == "not found") { 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"){ outputDir = ""; - outputDir += hasPath(fastafile); //if user entered a file with a path then preserve it + outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it } - string temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; } - convert(temp, processors); + string temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } + m->setProcessors(temp); + convert(temp, processors); } @@ -71,34 +140,17 @@ SeqSummaryCommand::SeqSummaryCommand(string option) { exit(1); } } -//********************************************************************************************************************** - -void SeqSummaryCommand::help(){ - try { - m->mothurOut("The summary.seqs command reads a fastafile and summarizes the sequences.\n"); - m->mothurOut("The summary.seqs command parameters are fasta and processors, fasta is required.\n"); - m->mothurOut("The summary.seqs command should be in the following format: \n"); - m->mothurOut("summary.seqs(fasta=yourFastaFile, processors=2) \n"); - m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n\n"); - } - catch(exception& e) { - m->errorOut(e, "SeqSummaryCommand", "help"); - exit(1); - } -} - -//*************************************************************************************************************** - -SeqSummaryCommand::~SeqSummaryCommand(){ /* do nothing */ } - //*************************************************************************************************************** int SeqSummaryCommand::execute(){ try{ - if (abort == true) { return 0; } + if (abort == true) { if (calledHelp) { return 0; } return 2; } - string summaryFile = outputDir + getSimpleName(fastafile) + ".summary"; + //set current fasta to fastafile + m->setFastaFile(fastafile); + + string summaryFile = outputDir + m->getSimpleName(fastafile) + ".summary"; int numSeqs = 0; @@ -107,6 +159,10 @@ int SeqSummaryCommand::execute(){ vector seqLength; vector ambigBases; vector longHomoPolymer; + + if (namefile != "") { nameMap = m->readNames(namefile); } + + if (m->control_pressed) { return 0; } #ifdef USE_MPI int pid, numSeqsPerProcessor; @@ -135,7 +191,7 @@ int SeqSummaryCommand::execute(){ if (pid == 0) { //you are the root process //print header - string outputString = "seqname\tstart\tend\tnbases\tambigs\tpolymer\n"; + string outputString = "seqname\tstart\tend\tnbases\tambigs\tpolymer\tnumSeqs\n"; int length = outputString.length(); char* buf2 = new char[length]; memcpy(buf2, outputString.c_str(), length); @@ -143,7 +199,7 @@ int SeqSummaryCommand::execute(){ MPI_File_write_shared(outMPI, buf2, length, MPI_CHAR, &statusOut); delete buf2; - MPIPos = setFilePosFasta(fastafile, numSeqs); //fills MPIPos, returns numSeqs + MPIPos = m->setFilePosFasta(fastafile, numSeqs); //fills MPIPos, returns numSeqs for(int i = 1; i < processors; i++) { MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD); @@ -224,7 +280,7 @@ int SeqSummaryCommand::execute(){ MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case #else - vector positions = divideFile(fastafile, processors); + vector positions = m->divideFile(fastafile, processors); for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(new linePair(positions[i], positions[(i+1)])); @@ -239,7 +295,7 @@ int SeqSummaryCommand::execute(){ rename((summaryFile + toString(processIDS[0]) + ".temp").c_str(), summaryFile.c_str()); //append files for(int i=1;iappendFiles((summaryFile + toString(processIDS[i]) + ".temp"), summaryFile); remove((summaryFile + toString(processIDS[i]) + ".temp").c_str()); } } @@ -260,13 +316,14 @@ int SeqSummaryCommand::execute(){ sort(seqLength.begin(), seqLength.end()); sort(ambigBases.begin(), ambigBases.end()); sort(longHomoPolymer.begin(), longHomoPolymer.end()); - - int ptile0_25 = int(numSeqs * 0.025); - int ptile25 = int(numSeqs * 0.250); - int ptile50 = int(numSeqs * 0.500); - int ptile75 = int(numSeqs * 0.750); - int ptile97_5 = int(numSeqs * 0.975); - int ptile100 = numSeqs - 1; + int size = startPosition.size(); + + int ptile0_25 = int(size * 0.025); + int ptile25 = int(size * 0.250); + int ptile50 = int(size * 0.500); + int ptile75 = int(size * 0.750); + int ptile97_5 = int(size * 0.975); + int ptile100 = size - 1; //to compensate for blank sequences that would result in startPosition and endPostion equalling -1 if (startPosition[0] == -1) { startPosition[0] = 0; } @@ -283,13 +340,14 @@ int SeqSummaryCommand::execute(){ 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("# of Seqs:\t" + toString(numSeqs)); 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; } m->mothurOutEndLine(); m->mothurOut("Output File Name: "); m->mothurOutEndLine(); - m->mothurOut(summaryFile); m->mothurOutEndLine(); + m->mothurOut(summaryFile); m->mothurOutEndLine(); outputNames.push_back(summaryFile); outputTypes["summary"].push_back(summaryFile); m->mothurOutEndLine(); #ifdef USE_MPI @@ -308,15 +366,15 @@ int SeqSummaryCommand::driverCreateSummary(vector& startPosition, vectoropenOutputFile(sumFile, outSummary); //print header if you are process 0 if (filePos->start == 0) { - outSummary << "seqname\tstart\tend\tnbases\tambigs\tpolymer" << endl; + outSummary << "seqname\tstart\tend\tnbases\tambigs\tpolymer\tnumSeqs" << endl; } ifstream in; - openInputFile(filename, in); + m->openInputFile(filename, in); in.seekg(filePos->start); @@ -327,30 +385,47 @@ int SeqSummaryCommand::driverCreateSummary(vector& startPosition, vectorcontrol_pressed) { in.close(); outSummary.close(); return 1; } - Sequence current(in); gobble(in); + Sequence current(in); m->gobble(in); if (current.getName() != "") { - startPosition.push_back(current.getStartPos()); - endPosition.push_back(current.getEndPos()); - seqLength.push_back(current.getNumBases()); - ambigBases.push_back(current.getAmbigBases()); - longHomoPolymer.push_back(current.getLongHomoPolymer()); + int num = 1; + if (namefile != "") { + //make sure this sequence is in the namefile, else error + map::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; } + else { num = it->second; } + } + + //for each sequence this sequence represents + for (int i = 0; i < num; i++) { + startPosition.push_back(current.getStartPos()); + endPosition.push_back(current.getEndPos()); + seqLength.push_back(current.getNumBases()); + ambigBases.push_back(current.getAmbigBases()); + longHomoPolymer.push_back(current.getLongHomoPolymer()); + } + + count++; outSummary << current.getName() << '\t'; outSummary << current.getStartPos() << '\t' << current.getEndPos() << '\t'; outSummary << current.getNumBases() << '\t' << current.getAmbigBases() << '\t'; - outSummary << current.getLongHomoPolymer() << endl; - count++; + outSummary << current.getLongHomoPolymer() << '\t' << num << endl; } - unsigned long int pos = in.tellg(); - if ((pos == -1) || (pos >= filePos->end)) { break; } + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + unsigned long int pos = in.tellg(); + if ((pos == -1) || (pos >= filePos->end)) { break; } + #else + if (in.eof()) { break; } + #endif //report progress - if((count) % 100 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } + //if((count) % 100 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } } //report progress - if((count) % 100 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } + //if((count) % 100 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } in.close(); @@ -370,7 +445,6 @@ int SeqSummaryCommand::MPICreateSummary(int start, int num, vector& startPo MPI_Status status; MPI_Comm_rank(MPI_COMM_WORLD, &pid); - for(int i=0;icontrol_pressed) { return 0; } @@ -389,14 +463,27 @@ int SeqSummaryCommand::MPICreateSummary(int start, int num, vector& startPo Sequence current(iss); if (current.getName() != "") { - startPosition.push_back(current.getStartPos()); - endPosition.push_back(current.getEndPos()); - seqLength.push_back(current.getNumBases()); - ambigBases.push_back(current.getAmbigBases()); - longHomoPolymer.push_back(current.getLongHomoPolymer()); + + int num = 1; + if (namefile != "") { + //make sure this sequence is in the namefile, else error + map::iterator it = nameMap.find(current.getName()); + + if (it == nameMap.end()) { cout << "[ERROR]: " << current.getName() << " is not in your namefile, please correct." << endl; m->control_pressed = true; } + else { num = it->second; } + } + + //for each sequence this sequence represents + for (int i = 0; i < num; i++) { + startPosition.push_back(current.getStartPos()); + endPosition.push_back(current.getEndPos()); + seqLength.push_back(current.getNumBases()); + ambigBases.push_back(current.getAmbigBases()); + longHomoPolymer.push_back(current.getLongHomoPolymer()); + } string outputString = current.getName() + "\t" + toString(current.getStartPos()) + "\t" + toString(current.getEndPos()) + "\t"; - outputString += toString(current.getNumBases()) + "\t" + toString(current.getAmbigBases()) + "\t" + toString(current.getLongHomoPolymer()) + "\n"; + outputString += toString(current.getNumBases()) + "\t" + toString(current.getAmbigBases()) + "\t" + toString(current.getLongHomoPolymer()) + "\t" + toString(num) + "\n"; //output to file length = outputString.length(); @@ -436,10 +523,11 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector& startPosition, //pass numSeqs to parent ofstream out; - string tempFile = toString(getpid()) + ".temp"; - openOutputFile(tempFile, out); + string tempFile = fastafile + toString(getpid()) + ".num.temp"; + m->openOutputFile(tempFile, out); out << num << endl; + out << startPosition.size() << endl; for (int k = 0; k < startPosition.size(); k++) { out << startPosition[k] << '\t'; } out << endl; for (int k = 0; k < endPosition.size(); k++) { out << endPosition[k] << '\t'; } out << endl; for (int k = 0; k < seqLength.size(); k++) { out << seqLength[k] << '\t'; } out << endl; @@ -449,7 +537,11 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector& startPosition, out.close(); exit(0); - }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); } + }else { + m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); + for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } + exit(0); + } } //force parent to wait until all the processes are done @@ -460,17 +552,18 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector& startPosition, //parent reads in and combine Filter info for (int i = 0; i < processIDS.size(); i++) { - string tempFilename = toString(processIDS[i]) + ".temp"; + string tempFilename = fastafile + toString(processIDS[i]) + ".num.temp"; ifstream in; - openInputFile(tempFilename, in); + m->openInputFile(tempFilename, in); int temp, tempNum; - in >> tempNum; gobble(in); num += tempNum; - for (int k = 0; k < tempNum; k++) { in >> temp; startPosition.push_back(temp); } gobble(in); - for (int k = 0; k < tempNum; k++) { in >> temp; endPosition.push_back(temp); } gobble(in); - for (int k = 0; k < tempNum; k++) { in >> temp; seqLength.push_back(temp); } gobble(in); - for (int k = 0; k < tempNum; k++) { in >> temp; ambigBases.push_back(temp); } gobble(in); - for (int k = 0; k < tempNum; k++) { in >> temp; longHomoPolymer.push_back(temp); } gobble(in); + in >> tempNum; m->gobble(in); num += tempNum; + in >> tempNum; m->gobble(in); + for (int k = 0; k < tempNum; k++) { in >> temp; startPosition.push_back(temp); } m->gobble(in); + for (int k = 0; k < tempNum; k++) { in >> temp; endPosition.push_back(temp); } m->gobble(in); + for (int k = 0; k < tempNum; k++) { in >> temp; seqLength.push_back(temp); } m->gobble(in); + for (int k = 0; k < tempNum; k++) { in >> temp; ambigBases.push_back(temp); } m->gobble(in); + for (int k = 0; k < tempNum; k++) { in >> temp; longHomoPolymer.push_back(temp); } m->gobble(in); in.close(); remove(tempFilename.c_str()); @@ -484,6 +577,7 @@ int SeqSummaryCommand::createProcessesCreateSummary(vector& startPosition, exit(1); } } -//*************************************************************************************************************** +/**********************************************************************************************************************/ +