X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;ds=sidebyside;f=trimflowscommand.cpp;h=05ed31ebd880af06984fa87cba97349df3f5e8c4;hb=fd98ee6efb944d38bbd61fc36ea9fea2557e3830;hp=e4d0165d136da966addcf3f235a162cad9f5baa1;hpb=5694c92fbf646fe01abc87bde2af59e14a9a56b6;p=mothur.git diff --git a/trimflowscommand.cpp b/trimflowscommand.cpp index e4d0165..05ed31e 100644 --- a/trimflowscommand.cpp +++ b/trimflowscommand.cpp @@ -11,52 +11,47 @@ #include "needlemanoverlap.hpp" //********************************************************************************************************************** - -vector TrimFlowsCommand::getValidParameters(){ +vector TrimFlowsCommand::setParameters(){ try { - string Array[] = {"flow", "totalflows", "minflows", - "fasta", "minlength", "maxlength", "maxhomop", "signal", "noise" - "oligos", "pdiffs", "bdiffs", "tdiffs", - "allfiles", "processors", - - - -// "group", - "outputdir","inputdir" + CommandParameter pflow("flow", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pflow); + CommandParameter poligos("oligos", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(poligos); + CommandParameter pmaxhomop("maxhomop", "Number", "", "9", "", "", "",false,false); parameters.push_back(pmaxhomop); + CommandParameter pmaxflows("maxflows", "Number", "", "720", "", "", "",false,false); parameters.push_back(pmaxflows); + CommandParameter pminflows("minflows", "Number", "", "360", "", "", "",false,false); parameters.push_back(pminflows); + CommandParameter pminlength("minlength", "Number", "", "0", "", "", "",false,false); parameters.push_back(pminlength); + CommandParameter pmaxlength("maxlength", "Number", "", "0", "", "", "",false,false); parameters.push_back(pmaxlength); + CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(ppdiffs); + CommandParameter pbdiffs("bdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(pbdiffs); + CommandParameter ptdiffs("tdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(ptdiffs); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); + CommandParameter psignal("signal", "Number", "", "0.50", "", "", "",false,false); parameters.push_back(psignal); + CommandParameter pnoise("noise", "Number", "", "0.70", "", "", "",false,false); parameters.push_back(pnoise); + CommandParameter pallfiles("allfiles", "Boolean", "", "t", "", "", "",false,false); parameters.push_back(pallfiles); + CommandParameter porder("order", "String", "", "", "", "", "",false,false); parameters.push_back(porder); + CommandParameter pfasta("fasta", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pfasta); + CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); + CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); - }; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); - return myArray; - } - catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "getValidParameters"); - exit(1); - } -} - -//********************************************************************************************************************** - -vector TrimFlowsCommand::getRequiredParameters(){ - try { - string Array[] = {"flow"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + vector myArray; + for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } return myArray; } catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "getRequiredParameters"); + m->errorOut(e, "TrimFlowsCommand", "setParameters"); exit(1); } } - //********************************************************************************************************************** - -vector TrimFlowsCommand::getRequiredFiles(){ +string TrimFlowsCommand::getHelpString(){ try { - vector myArray; - return myArray; + string helpString = ""; + helpString += "The trim.flows command reads a flowgram file and creates .....\n"; + helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n"; + helpString += "For more details please check out the wiki http://www.mothur.org/wiki/Trim.flows.\n"; + return helpString; } catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "getRequiredFiles"); + m->errorOut(e, "TrimFlowsCommand", "getHelpString"); exit(1); } } @@ -65,8 +60,8 @@ vector TrimFlowsCommand::getRequiredFiles(){ TrimFlowsCommand::TrimFlowsCommand(){ try { - abort = true; - //initialize outputTypes + abort = true; calledHelp = true; + setParameters(); vector tempOutNames; outputTypes["flow"] = tempOutNames; outputTypes["fasta"] = tempOutNames; @@ -76,50 +71,21 @@ TrimFlowsCommand::TrimFlowsCommand(){ exit(1); } } - -//*************************************************************************************************************** - -TrimFlowsCommand::~TrimFlowsCommand(){ /* do nothing */ } - -//*************************************************************************************************************** - -void TrimFlowsCommand::help(){ - try { - m->mothurOut("The trim.flows command reads a flowgram file and creates .....\n"); - m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n"); - m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/Trim.flows.\n\n"); - - } - catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "help"); - exit(1); - } -} - //********************************************************************************************************************** TrimFlowsCommand::TrimFlowsCommand(string option) { try { - abort = false; + abort = false; calledHelp = false; comboStarts = 0; //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 AlignArray[] = {"flow", "totalflows", "minflows", - "fasta", "minlength", "maxlength", "maxhomop", "signal", "noise", - "oligos", "pdiffs", "bdiffs", "tdiffs", - "allfiles", "processors", - - // "group", - "outputdir","inputdir" - - }; - - vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); + + vector myArray = setParameters(); OptionParser parser(option); map parameters = parser.getParameters(); @@ -171,8 +137,14 @@ TrimFlowsCommand::TrimFlowsCommand(string option) { //check for required parameters flowFileName = validParameter.validFile(parameters, "flow", true); - if (flowFileName == "not found") { m->mothurOut("flow is a required parameter for the trim.flows command."); m->mothurOutEndLine(); abort = true; } - else if (flowFileName == "not open") { abort = true; } + if (flowFileName == "not found") { + flowFileName = m->getFlowFile(); + if (flowFileName != "") { m->mothurOut("Using " + flowFileName + " as input file for the flow parameter."); m->mothurOutEndLine(); } + else { + m->mothurOut("No valid current flow file. You must provide a flow file."); m->mothurOutEndLine(); + abort = true; + } + }else if (flowFileName == "not open") { flowFileName = ""; abort = true; } //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"){ @@ -185,11 +157,11 @@ TrimFlowsCommand::TrimFlowsCommand(string option) { // ...at some point should added some additional type checking... string temp; - temp = validParameter.validFile(parameters, "minflows", false); if (temp == "not found") { temp = "360"; } + temp = validParameter.validFile(parameters, "minflows", false); if (temp == "not found") { temp = "360"; } convert(temp, minFlows); - temp = validParameter.validFile(parameters, "totalflows", false); if (temp == "not found") { temp = "720"; } - convert(temp, totalFlows); + temp = validParameter.validFile(parameters, "maxflows", false); if (temp == "not found") { temp = "720"; } + convert(temp, maxFlows); temp = validParameter.validFile(parameters, "oligos", true); @@ -226,12 +198,19 @@ TrimFlowsCommand::TrimFlowsCommand(string option) { convert(temp, tdiffs); if(tdiffs == 0){ tdiffs = bdiffs + pdiffs; } - temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found") { temp = "T"; } + temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found"){ temp = "T"; } allFiles = m->isTrue(temp); - temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found") { temp = "1"; } - convert(temp, processors); - + temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } + m->setProcessors(temp); + convert(temp, processors); + + flowOrder = validParameter.validFile(parameters, "order", false); + if (flowOrder == "not found"){ flowOrder = "TACG"; } + else if(flowOrder.length() != 4){ + m->mothurOut("The value of the order option must be four bases long\n"); + } + if(oligoFileName == ""){ allFiles = 0; } numFPrimers = 0; @@ -250,7 +229,7 @@ TrimFlowsCommand::TrimFlowsCommand(string option) { int TrimFlowsCommand::execute(){ try{ - if (abort == true) { return 0; } + if (abort == true) { if (calledHelp) { return 0; } return 2; } string trimFlowFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "trim.flow"; outputNames.push_back(trimFlowFileName); outputTypes["flow"].push_back(trimFlowFileName); @@ -285,10 +264,12 @@ int TrimFlowsCommand::execute(){ if (m->control_pressed) { return 0; } + string flowFilesFileName; + ofstream output; + if(allFiles){ - ofstream output; - string flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files"; + flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files"; m->openOutputFile(flowFilesFileName, output); for(int i=0;imothurOut("deleting: " + barcodePrimerComboFileNames[i][j] + '\n'); remove(barcodePrimerComboFileNames[i][j].c_str()); } else{ -// m->mothurOut("saving: " + barcodePrimerComboFileNames[i][j] + '\n'); output << barcodePrimerComboFileNames[i][j] << endl; } } } output.close(); } + else{ + flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files"; + m->openOutputFile(flowFilesFileName, output); + + output << trimFlowFileName << endl; + + output.close(); + } + outputTypes["flow.files"].push_back(flowFilesFileName); + outputNames.push_back(flowFileName); + + //set fasta file as new current fastafile + string current = ""; + itTypes = outputTypes.find("fasta"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); } + } m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); @@ -347,20 +343,22 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN scrapFlowFile.setf(ios::fixed, ios::floatfield); scrapFlowFile.setf(ios::showpoint); if(line->start == 4){ - scrapFlowFile << totalFlows << endl; - trimFlowFile << totalFlows << endl; - for(int i=0;iopenOutputFile(barcodePrimerComboFileNames[i][j], temp); - temp << totalFlows << endl; - temp.close(); - } - } + scrapFlowFile << maxFlows << endl; + trimFlowFile << maxFlows << endl; + if(allFiles){ + for(int i=0;iopenOutputFile(barcodePrimerComboFileNames[i][j], temp); + temp << maxFlows << endl; + temp.close(); + } + } + } } - FlowData flowData(numFlows, signal, noise, maxHomoP); + FlowData flowData(numFlows, signal, noise, maxHomoP, flowOrder); ofstream fastaFile; if(fasta){ m->openOutputFile(fastaFileName, fastaFile); } @@ -372,7 +370,7 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN int count = 0; bool moreSeqs = 1; - + while(moreSeqs) { int success = 1; @@ -380,7 +378,7 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN string trashCode = ""; flowData.getNext(flowFile); - flowData.capFlows(totalFlows); + flowData.capFlows(maxFlows); Sequence currSeq = flowData.getSequence(); if(!flowData.hasMinFlows(minFlows)){ //screen to see if sequence is of a minimum number of flows @@ -396,7 +394,8 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN } } - int primerIndex, barcodeIndex; + int primerIndex = 0; + int barcodeIndex = 0; if(barcodes.size() != 0){ success = stripBarcode(currSeq, barcodeIndex); @@ -416,25 +415,21 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN success = stripReverse(currSeq); if(!success) { trashCode += 'r'; } } - - + if(trashCode.length() == 0){ + flowData.printFlows(trimFlowFile); + + if(fasta) { currSeq.printSequence(fastaFile); } - if(fasta){ - currSeq.printSequence(fastaFile); - } - - if(barcodes.size() != 0 || primers.size() != 0){ -// string fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + barcodePrimerCombos[barcodeIndex][primerIndex] + ".flow"; + if(allFiles){ ofstream output; m->openOutputFileAppend(barcodePrimerComboFileNames[barcodeIndex][primerIndex], output); output.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint); flowData.printFlows(output); output.close(); - } - + } } else{ flowData.printFlows(scrapFlowFile, trashCode); @@ -542,13 +537,27 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ } oligosFile.close(); - //add in potential combos + if(barcodeNameVector.size() == 0 && primerNameVector[0] == ""){ allFiles = 0; } + + //add in potential combos + if(barcodeNameVector.size() == 0){ + barcodes[""] = 0; + barcodeNameVector.push_back(""); + } + + if(primerNameVector.size() == 0){ + primers[""] = 0; + primerNameVector.push_back(""); + } + + outFlowFileNames.resize(barcodeNameVector.size()); for(int i=0;i::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){ for(map::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){ @@ -563,7 +572,12 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow"; } else{ - comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second]; + if(barcodeName == ""){ + comboGroupName = primerNameVector[itPrimer->second]; + } + else{ + comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second]; + } fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow"; } @@ -709,7 +723,7 @@ int TrimFlowsCommand::stripForward(Sequence& seq, int& group){ string rawSequence = seq.getUnaligned(); int success = pdiffs + 1; //guilty until proven innocent - + //can you find the primer for(map::iterator it=primers.begin();it!=primers.end();it++){ string oligo = it->first; @@ -727,7 +741,7 @@ int TrimFlowsCommand::stripForward(Sequence& seq, int& group){ } //if you found the barcode or if you don't want to allow for diffs - if ((pdiffs == 0) || (success == 0)) { return success; } + if ((pdiffs == 0) || (success == 0)) { return success; } else { //try aligning and see if you can find it @@ -1017,16 +1031,17 @@ int TrimFlowsCommand::createProcessesCreateTrim(string flowFileName, string trim }else if (pid == 0){ vector > tempBarcodePrimerComboFileNames = barcodePrimerComboFileNames; - for(int i=0;iopenOutputFile(tempBarcodePrimerComboFileNames[i][j], temp); - temp.close(); - + if(allFiles){ + for(int i=0;iopenOutputFile(tempBarcodePrimerComboFileNames[i][j], temp); + temp.close(); + + } } } - driverCreateTrim(flowFileName, (trimFlowFileName + toString(getpid()) + ".temp"), (scrapFlowFileName + toString(getpid()) + ".temp"),