X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=trimflowscommand.cpp;h=33349decab52e90be19c928b44b9b07c8f6c710a;hp=960a59d3950e92590063f348f890c267a3123e12;hb=b206f634aae1b4ce13978d203247fb64757d5482;hpb=902e0fcab76e75009ac43d3f4537e08182628d6f diff --git a/trimflowscommand.cpp b/trimflowscommand.cpp index 960a59d..33349de 100644 --- a/trimflowscommand.cpp +++ b/trimflowscommand.cpp @@ -10,87 +10,86 @@ #include "trimflowscommand.h" #include "needlemanoverlap.hpp" -//********************************************************************************************************************** -vector TrimFlowsCommand::getValidParameters(){ +//********************************************************************************************************************** +vector TrimFlowsCommand::setParameters(){ try { - string Array[] = {"flow", "maxflows", "minflows", - "fasta", "minlength", "maxlength", "maxhomop", "signal", "noise" - "oligos", "pdiffs", "bdiffs", "tdiffs", - "allfiles", "processors", - "outputdir","inputdir" + CommandParameter pflow("flow", "InputTypes", "", "", "none", "none", "none","flow-file",false,true,true); parameters.push_back(pflow); + CommandParameter poligos("oligos", "InputTypes", "", "", "none", "none", "none","",false,false,true); parameters.push_back(poligos); + CommandParameter pmaxhomop("maxhomop", "Number", "", "9", "", "", "","",false,false); parameters.push_back(pmaxhomop); + CommandParameter pmaxflows("maxflows", "Number", "", "450", "", "", "","",false,false); parameters.push_back(pmaxflows); + CommandParameter pminflows("minflows", "Number", "", "450", "", "", "","",false,false); parameters.push_back(pminflows); + CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(ppdiffs); + CommandParameter pbdiffs("bdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(pbdiffs); + CommandParameter pldiffs("ldiffs", "Number", "", "0", "", "", "","",false,false); parameters.push_back(pldiffs); + CommandParameter psdiffs("sdiffs", "Number", "", "0", "", "", "","",false,false); parameters.push_back(psdiffs); + CommandParameter ptdiffs("tdiffs", "Number", "", "0", "", "", "","",false,false); parameters.push_back(ptdiffs); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); 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", "Multiple", "A-B-I", "A", "", "", "","",false,false, true); 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))); + 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", "getValidParameters"); + m->errorOut(e, "TrimFlowsCommand", "setParameters"); exit(1); } } - //********************************************************************************************************************** - -vector TrimFlowsCommand::getRequiredParameters(){ +string TrimFlowsCommand::getHelpString(){ try { - string Array[] = {"flow"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); - return myArray; + string helpString = ""; + helpString += "The trim.flows command reads a flowgram file and creates .....\n"; + helpString += "The order parameter options are A, B or I. Default=A. A = TACG and B = TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC and I = TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC.\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", "getRequiredParameters"); + m->errorOut(e, "TrimFlowsCommand", "getHelpString"); exit(1); } } - //********************************************************************************************************************** - -vector TrimFlowsCommand::getRequiredFiles(){ - try { - vector myArray; - return myArray; - } - catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "getRequiredFiles"); - exit(1); - } +string TrimFlowsCommand::getOutputPattern(string type) { + try { + string pattern = ""; + + if (type == "flow") { pattern = "[filename],[tag],flow"; } + else if (type == "fasta") { pattern = "[filename],flow.fasta"; } + else if (type == "file") { pattern = "[filename],flow.files"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } + + return pattern; + } + catch(exception& e) { + m->errorOut(e, "TrimFlowsCommand", "getOutputPattern"); + exit(1); + } } - //********************************************************************************************************************** TrimFlowsCommand::TrimFlowsCommand(){ try { abort = true; calledHelp = true; + setParameters(); vector tempOutNames; outputTypes["flow"] = tempOutNames; outputTypes["fasta"] = tempOutNames; + outputTypes["file"] = tempOutNames; } catch(exception& e) { m->errorOut(e, "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) { @@ -101,20 +100,11 @@ TrimFlowsCommand::TrimFlowsCommand(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 { - //valid paramters for this command - string AlignArray[] = {"flow", "maxflows", "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(); @@ -131,6 +121,7 @@ TrimFlowsCommand::TrimFlowsCommand(string option) { vector tempOutNames; outputTypes["flow"] = tempOutNames; outputTypes["fasta"] = tempOutNames; + outputTypes["file"] = tempOutNames; //if the user changes the input directory command factory will send this info to us in the output parameter string inputDir = validParameter.validFile(parameters, "inputdir", false); @@ -153,21 +144,19 @@ TrimFlowsCommand::TrimFlowsCommand(string option) { if (path == "") { parameters["oligos"] = inputDir + it->second; } } - -// it = parameters.find("group"); -// //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["group"] = inputDir + it->second; } -// } } //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"){ @@ -180,62 +169,77 @@ 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"; } - convert(temp, minFlows); + temp = validParameter.validFile(parameters, "minflows", false); if (temp == "not found") { temp = "450"; } + m->mothurConvert(temp, minFlows); - temp = validParameter.validFile(parameters, "maxflows", false); if (temp == "not found") { temp = "720"; } - convert(temp, maxFlows); + temp = validParameter.validFile(parameters, "maxflows", false); if (temp == "not found") { temp = "450"; } + m->mothurConvert(temp, maxFlows); temp = validParameter.validFile(parameters, "oligos", true); if (temp == "not found") { oligoFileName = ""; } else if(temp == "not open") { abort = true; } - else { oligoFileName = temp; } + else { oligoFileName = temp; m->setOligosFile(oligoFileName); } temp = validParameter.validFile(parameters, "fasta", false); if (temp == "not found"){ fasta = 0; } else if(m->isTrue(temp)) { fasta = 1; } temp = validParameter.validFile(parameters, "maxhomop", false); if (temp == "not found"){ temp = "9"; } - convert(temp, maxHomoP); + m->mothurConvert(temp, maxHomoP); temp = validParameter.validFile(parameters, "signal", false); if (temp == "not found"){ temp = "0.50"; } - convert(temp, signal); + m->mothurConvert(temp, signal); temp = validParameter.validFile(parameters, "noise", false); if (temp == "not found"){ temp = "0.70"; } - convert(temp, noise); - - temp = validParameter.validFile(parameters, "minlength", false); if (temp == "not found"){ temp = "0"; } - convert(temp, minLength); - - temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found"){ temp = "0"; } - convert(temp, maxLength); - + m->mothurConvert(temp, noise); + temp = validParameter.validFile(parameters, "bdiffs", false); if (temp == "not found"){ temp = "0"; } - convert(temp, bdiffs); + m->mothurConvert(temp, bdiffs); temp = validParameter.validFile(parameters, "pdiffs", false); if (temp == "not found"){ temp = "0"; } - convert(temp, pdiffs); + m->mothurConvert(temp, pdiffs); - temp = validParameter.validFile(parameters, "tdiffs", false); - if (temp == "not found"){ int tempTotal = pdiffs + bdiffs; temp = toString(tempTotal); } - convert(temp, tdiffs); - if(tdiffs == 0){ tdiffs = bdiffs + pdiffs; } + temp = validParameter.validFile(parameters, "ldiffs", false); if (temp == "not found") { temp = "0"; } + m->mothurConvert(temp, ldiffs); + + temp = validParameter.validFile(parameters, "sdiffs", false); if (temp == "not found") { temp = "0"; } + m->mothurConvert(temp, sdiffs); - temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found") { temp = "T"; } - allFiles = m->isTrue(temp); + temp = validParameter.validFile(parameters, "tdiffs", false); if (temp == "not found") { int tempTotal = pdiffs + bdiffs + ldiffs + sdiffs; temp = toString(tempTotal); } + m->mothurConvert(temp, tdiffs); - temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found") { temp = "1"; } - convert(temp, processors); + if(tdiffs == 0){ tdiffs = bdiffs + pdiffs + ldiffs + sdiffs; } + - if(oligoFileName == ""){ allFiles = 0; } + temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } + m->setProcessors(temp); + m->mothurConvert(temp, processors); + + temp = validParameter.validFile(parameters, "order", false); if (temp == "not found"){ temp = "A"; } + if (temp.length() > 1) { m->mothurOut("[ERROR]: " + temp + " is not a valid option for order. order options are A, B, or I. A = TACG, B = TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC, and I = TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC.\n"); abort=true; + } + else { + if (toupper(temp[0]) == 'A') { flowOrder = "TACG"; } + else if(toupper(temp[0]) == 'B'){ + flowOrder = "TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC"; } + else if(toupper(temp[0]) == 'I'){ + flowOrder = "TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC"; } + else { + m->mothurOut("[ERROR]: " + temp + " is not a valid option for order. order options are A, B, or I. A = TACG, B = TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC, and I = TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC.\n"); abort=true; + } + } + + if(oligoFileName == "") { allFiles = 0; } + else { allFiles = 1; } numFPrimers = 0; numRPrimers = 0; + numLinkers = 0; + numSpacers = 0; } - } catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "TrimSeqsCommand"); + m->errorOut(e, "TrimFlowsCommand", "TrimFlowsCommand"); exit(1); } } @@ -247,36 +251,61 @@ int TrimFlowsCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } - string trimFlowFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "trim.flow"; + map variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName)); + string fastaFileName = getOutputFileName("fasta",variables); + if(fasta){ outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName); } + + variables["[tag]"] = "trim"; + string trimFlowFileName = getOutputFileName("flow",variables); outputNames.push_back(trimFlowFileName); outputTypes["flow"].push_back(trimFlowFileName); - string scrapFlowFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "scrap.flow"; + variables["[tag]"] = "scrap"; + string scrapFlowFileName = getOutputFileName("flow",variables); outputNames.push_back(scrapFlowFileName); outputTypes["flow"].push_back(scrapFlowFileName); - string fastaFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.fasta"; - if(fasta){ - outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName); - } - vector flowFilePos = getFlowFileBreaks(); + + vector flowFilePos; + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + flowFilePos = getFlowFileBreaks(); for (int i = 0; i < (flowFilePos.size()-1); i++) { lines.push_back(new linePair(flowFilePos[i], flowFilePos[(i+1)])); } - + #else + ifstream in; m->openInputFile(flowFileName, in); in >> numFlows; in.close(); + ///////////////////////////////////////// until I fix multiple processors for windows ////////////////// + processors = 1; + ///////////////////////////////////////// until I fix multiple processors for windows ////////////////// + if (processors == 1) { + lines.push_back(new linePair(0, 1000)); + }else { + int numFlowLines; + flowFilePos = m->setFilePosEachLine(flowFileName, numFlowLines); + flowFilePos.erase(flowFilePos.begin() + 1); numFlowLines--; + + //figure out how many sequences you have to process + int numSeqsPerProcessor = numFlowLines / processors; + cout << numSeqsPerProcessor << '\t' << numFlowLines << endl; + for (int i = 0; i < processors; i++) { + int startIndex = i * numSeqsPerProcessor; + if(i == (processors - 1)){ numSeqsPerProcessor = numFlowLines - i * numSeqsPerProcessor; } + lines.push_back(new linePair(flowFilePos[startIndex], numSeqsPerProcessor)); + cout << flowFilePos[startIndex] << '\t' << numSeqsPerProcessor << endl; + } + } + #endif + vector > barcodePrimerComboFileNames; if(oligoFileName != ""){ getOligos(barcodePrimerComboFileNames); } -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) if(processors == 1){ driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]); }else{ createProcessesCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames); } -#else - driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]); -#endif if (m->control_pressed) { return 0; } @@ -284,45 +313,52 @@ int TrimFlowsCommand::execute(){ ofstream output; if(allFiles){ - - flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files"; + set namesAlreadyProcessed; + flowFilesFileName = getOutputFileName("file",variables); m->openOutputFile(flowFilesFileName, output); for(int i=0;imothurRemove(barcodePrimerComboFileNames[i][j]); + } + else{ + output << m->getFullPathName(barcodePrimerComboFileNames[i][j]) << endl; + outputNames.push_back(barcodePrimerComboFileNames[i][j]); + outputTypes["flow"].push_back(barcodePrimerComboFileNames[i][j]); + } + namesAlreadyProcessed.insert(barcodePrimerComboFileNames[i][j]); + } } } } output.close(); } else{ - flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files"; + flowFilesFileName = getOutputFileName("file",variables); m->openOutputFile(flowFilesFileName, output); - output << trimFlowFileName << endl; + output << m->getFullPathName(trimFlowFileName) << endl; output.close(); } - outputTypes["flow.files"].push_back(flowFilesFileName); - + outputTypes["file"].push_back(flowFilesFileName); + outputNames.push_back(flowFilesFileName); + m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } @@ -338,10 +374,9 @@ int TrimFlowsCommand::execute(){ //*************************************************************************************************************** -int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileName, string scrapFlowFileName, string fastaFileName, vector > barcodePrimerComboFileNames, linePair* line){ +int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileName, string scrapFlowFileName, string fastaFileName, vector > thisBarcodePrimerComboFileNames, linePair* line){ try { - ofstream trimFlowFile; m->openOutputFile(trimFlowFileName, trimFlowFile); trimFlowFile.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint); @@ -349,24 +384,6 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN ofstream scrapFlowFile; m->openOutputFile(scrapFlowFileName, scrapFlowFile); scrapFlowFile.setf(ios::fixed, ios::floatfield); scrapFlowFile.setf(ios::showpoint); - - if(line->start == 4){ - 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); ofstream fastaFile; if(fasta){ m->openOutputFile(fastaFileName, fastaFile); } @@ -375,44 +392,81 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN m->openInputFile(flowFileName, flowFile); flowFile.seekg(line->start); - + + if(line->start == 0){ + flowFile >> numFlows; m->gobble(flowFile); + scrapFlowFile << numFlows << endl; + trimFlowFile << maxFlows << endl; + if(allFiles){ + for(int i=0;iopenOutputFile(thisBarcodePrimerComboFileNames[i][j], temp); + temp << maxFlows << endl; + temp.close(); + } + } + } + } + } + + FlowData flowData(numFlows, signal, noise, maxHomoP, flowOrder); + //cout << " driver flowdata address " << &flowData << &flowFile << endl; int count = 0; bool moreSeqs = 1; - + + TrimOligos trimOligos(pdiffs, bdiffs, ldiffs, sdiffs, primers, barcodes, revPrimer, linker, spacer); + while(moreSeqs) { + + if (m->control_pressed) { break; } int success = 1; int currentSeqDiffs = 0; string trashCode = ""; - flowData.getNext(flowFile); + flowData.getNext(flowFile); flowData.capFlows(maxFlows); Sequence currSeq = flowData.getSequence(); + //cout << currSeq.getName() << '\t' << currSeq.getUnaligned() << endl; if(!flowData.hasMinFlows(minFlows)){ //screen to see if sequence is of a minimum number of flows success = 0; trashCode += 'l'; } - - if(minLength > 0 || maxLength > 0){ //screen to see if sequence is above and below a specific number of bases - int seqLength = currSeq.getNumBases(); - if(seqLength < minLength || seqLength > maxLength){ - success = 0; - trashCode += 'l'; - } + if(!flowData.hasGoodHomoP()){ //screen to see if sequence meets the maximum homopolymer limit + success = 0; + trashCode += 'h'; } - + int primerIndex = 0; int barcodeIndex = 0; + if(numLinkers != 0){ + success = trimOligos.stripLinker(currSeq); + if(success > ldiffs) { trashCode += 'k'; } + else{ currentSeqDiffs += success; } + + } + + if (m->debug) { m->mothurOut("[DEBUG]: " + currSeq.getName() + " " + currSeq.getUnaligned() + "\n"); } + if(barcodes.size() != 0){ - success = stripBarcode(currSeq, barcodeIndex); + success = trimOligos.stripBarcode(currSeq, barcodeIndex); if(success > bdiffs) { trashCode += 'b'; } else{ currentSeqDiffs += success; } } + if(numSpacers != 0){ + success = trimOligos.stripSpacer(currSeq); + if(success > sdiffs) { trashCode += 's'; } + else{ currentSeqDiffs += success; } + + } + if(numFPrimers != 0){ - success = stripForward(currSeq, primerIndex); + success = trimOligos.stripForward(currSeq, primerIndex); if(success > pdiffs) { trashCode += 'f'; } else{ currentSeqDiffs += success; } } @@ -420,38 +474,52 @@ int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileN if (currentSeqDiffs > tdiffs) { trashCode += 't'; } if(numRPrimers != 0){ - success = stripReverse(currSeq); + success = trimOligos.stripReverse(currSeq); if(!success) { trashCode += 'r'; } } - - if(trashCode.length() == 0){ - - flowData.printFlows(trimFlowFile); - if(fasta) { currSeq.printSequence(fastaFile); } - - if(allFiles){ -// string fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + barcodePrimerCombos[barcodeIndex][primerIndex] + ".flow"; - ofstream output; - m->openOutputFileAppend(barcodePrimerComboFileNames[barcodeIndex][primerIndex], output); - output.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint); - - flowData.printFlows(output); - output.close(); - } - + if(trashCode.length() == 0){ + string thisGroup = ""; + if(barcodes.size() != 0){ + thisGroup = barcodeNameVector[barcodeIndex]; + if (primers.size() != 0) { + if (primerNameVector[primerIndex] != "") { + if(thisGroup != "") { + thisGroup += "." + primerNameVector[primerIndex]; + }else { + thisGroup = primerNameVector[primerIndex]; + } + } + } + } + + int pos = thisGroup.find("ignore"); + if (pos == string::npos) { + flowData.printFlows(trimFlowFile); + + if(fasta) { currSeq.printSequence(fastaFile); } + + if(allFiles){ + ofstream output; + m->openOutputFileAppend(thisBarcodePrimerComboFileNames[barcodeIndex][primerIndex], output); + output.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint); + + flowData.printFlows(output); + output.close(); + } + } } else{ flowData.printFlows(scrapFlowFile, trashCode); } count++; - + //cout << "driver" << '\t' << currSeq.getName() << endl; //report progress if((count) % 10000 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - unsigned long int pos = flowFile.tellg(); +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + unsigned long long pos = flowFile.tellg(); if ((pos == -1) || (pos >= line->end)) { break; } #else @@ -489,13 +557,16 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ while(!oligosFile.eof()){ - oligosFile >> type; m->gobble(oligosFile); //get the first column value of the row - is it a comment or a feature we are interested in? - + oligosFile >> type; //get the first column value of the row - is it a comment or a feature we are interested in? + + if (m->debug) { m->mothurOut("[DEBUG]: type = " + type + ".\n"); } + if(type[0] == '#'){ //igore the line because there's a comment - while (!oligosFile.eof()) { char c = oligosFile.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there + while (!oligosFile.eof()) { char c = oligosFile.get(); if (c == 10 || c == 13){ break; } } + m->gobble(oligosFile);// get rest of line if there's any crap there } else{ //there's a feature we're interested in - + m->gobble(oligosFile); for(int i=0;i> oligo; //get the DNA sequence for the feature @@ -504,13 +575,15 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ oligo[i] = toupper(oligo[i]); if(oligo[i] == 'U') { oligo[i] = 'T'; } } - + + if (m->debug) { m->mothurOut("[DEBUG]: oligos = " + oligo + ".\n"); } + if(type == "FORWARD"){ //if the feature is a forward primer... group = ""; while (!oligosFile.eof()) { // get rest of line in case there is a primer name = will have the name of the primer char c = oligosFile.get(); - if (c == 10 || c == 13){ break; } + if (c == 10 || c == 13 || c == -1){ break; } else if (c == 32 || c == 9){;} //space or tab else { group += c; } } @@ -524,9 +597,9 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ } else if(type == "REVERSE"){ - Sequence oligoRC("reverse", oligo); - oligoRC.reverseComplement(); - revPrimer.push_back(oligoRC.getUnaligned()); + string oligoRC = reverseOligo(oligo); + revPrimer.push_back(oligoRC); + if (m->debug) { m->mothurOut("[DEBUG]: reverse oligos = " + oligoRC + ".\n"); } } else if(type == "BARCODE"){ oligosFile >> group; @@ -534,9 +607,15 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ //check for repeat barcodes map::iterator itBar = barcodes.find(oligo); if (itBar != barcodes.end()) { m->mothurOut("barcode " + oligo + " is in your oligos file already."); m->mothurOutEndLine(); } - + + if (m->debug) { m->mothurOut("[DEBUG]: group = " + group + ".\n"); } + barcodes[oligo]=indexBarcode; indexBarcode++; barcodeNameVector.push_back(group); + }else if(type == "LINKER"){ + linker.push_back(oligo); + }else if(type == "SPACER"){ + spacer.push_back(oligo); } else{ m->mothurOut(type + " is not recognized as a valid type. Choices are forward, reverse, and barcode. Ignoring " + oligo + "."); m->mothurOutEndLine(); @@ -559,8 +638,7 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ primers[""] = 0; primerNameVector.push_back(""); } - - + outFlowFileNames.resize(barcodeNameVector.size()); for(int i=0;i >& outFlowFileNames){ string primerName = primerNameVector[itPrimer->second]; string barcodeName = barcodeNameVector[itBar->second]; - - string comboGroupName = ""; - string fileName = ""; - - if(primerName == ""){ - comboGroupName = barcodeNameVector[itBar->second]; - fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow"; - } - else{ - if(barcodeName == ""){ - comboGroupName = primerNameVector[itPrimer->second]; - } - else{ - comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second]; - } - fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow"; - } - - outputNames.push_back(fileName); - outputTypes["flow"].push_back(fileName); - outFlowFileNames[itBar->second][itPrimer->second] = fileName; - - ofstream temp; - m->openOutputFile(fileName, temp); - temp.close(); + + if ((primerName == "ignore") || (barcodeName == "ignore")) { } //do nothing + else { + string comboGroupName = ""; + string fileName = ""; + + map variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName)); + + if(primerName == ""){ + comboGroupName = barcodeNameVector[itBar->second]; + variables["[tag]"] = comboGroupName; + fileName = getOutputFileName("flow", variables); + } + else{ + if(barcodeName == ""){ + comboGroupName = primerNameVector[itPrimer->second]; + } + else{ + comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second]; + } + variables["[tag]"] = comboGroupName; + fileName = getOutputFileName("flow", variables); + } + + outFlowFileNames[itBar->second][itPrimer->second] = fileName; + + ofstream temp; + m->openOutputFile(fileName, temp); + temp.close(); + } } } } numFPrimers = primers.size(); numRPrimers = revPrimer.size(); + numLinkers = linker.size(); + numSpacers = spacer.size(); } catch(exception& e) { @@ -611,352 +697,57 @@ void TrimFlowsCommand::getOligos(vector >& outFlowFileNames){ exit(1); } } - -//*************************************************************************************************************** - -int TrimFlowsCommand::stripBarcode(Sequence& seq, int& group){ - try { - - string rawSequence = seq.getUnaligned(); - int success = bdiffs + 1; //guilty until proven innocent - - //can you find the barcode - for(map::iterator it=barcodes.begin();it!=barcodes.end();it++){ - string oligo = it->first; - if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length - success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out - break; - } - - if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){ - group = it->second; - seq.setUnaligned(rawSequence.substr(oligo.length())); - - success = 0; - break; - } - } - - //if you found the barcode or if you don't want to allow for diffs - if ((bdiffs == 0) || (success == 0)) { return success; } - - else { //try aligning and see if you can find it - - int maxLength = 0; - - Alignment* alignment; - if (barcodes.size() > 0) { - map::iterator it=barcodes.begin(); - - for(it;it!=barcodes.end();it++){ - if(it->first.length() > maxLength){ - maxLength = it->first.length(); - } - } - alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1)); - - }else{ alignment = NULL; } - - //can you find the barcode - int minDiff = 1e6; - int minCount = 1; - int minGroup = -1; - int minPos = 0; - - for(map::iterator it=barcodes.begin();it!=barcodes.end();it++){ - string oligo = it->first; - // int length = oligo.length(); - - if(rawSequence.length() < maxLength){ //let's just assume that the barcodes are the same length - success = bdiffs + 10; - break; - } - - //use needleman to align first barcode.length()+numdiffs of sequence to each barcode - alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs)); - oligo = alignment->getSeqAAln(); - string temp = alignment->getSeqBAln(); - - int alnLength = oligo.length(); - - for(int i=oligo.length()-1;i>=0;i--){ - if(oligo[i] != '-'){ alnLength = i+1; break; } - } - oligo = oligo.substr(0,alnLength); - temp = temp.substr(0,alnLength); - - int numDiff = countDiffs(oligo, temp); - - if(numDiff < minDiff){ - minDiff = numDiff; - minCount = 1; - minGroup = it->second; - minPos = 0; - for(int i=0;i bdiffs) { success = minDiff; } //no good matches - else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes - else{ //use the best match - group = minGroup; - seq.setUnaligned(rawSequence.substr(minPos)); - success = minDiff; - } - - if (alignment != NULL) { delete alignment; } - - } - - return success; - - } - catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "stripBarcode"); - exit(1); - } - -} - -//*************************************************************************************************************** - -int TrimFlowsCommand::stripForward(Sequence& seq, int& group){ - try { - - 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; - if(rawSequence.length() < oligo.length()){ //let's just assume that the primers are the same length - success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out - break; - } - - if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){ - group = it->second; - seq.setUnaligned(rawSequence.substr(oligo.length())); - success = 0; - break; - } - } - - //if you found the barcode or if you don't want to allow for diffs - if ((pdiffs == 0) || (success == 0)) { return success; } - - else { //try aligning and see if you can find it - - int maxLength = 0; - - Alignment* alignment; - if (primers.size() > 0) { - map::iterator it=primers.begin(); - - for(it;it!=primers.end();it++){ - if(it->first.length() > maxLength){ - maxLength = it->first.length(); - } - } - alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1)); - - }else{ alignment = NULL; } - - //can you find the barcode - int minDiff = 1e6; - int minCount = 1; - int minGroup = -1; - int minPos = 0; - - for(map::iterator it=primers.begin();it!=primers.end();it++){ - string oligo = it->first; - // int length = oligo.length(); - - if(rawSequence.length() < maxLength){ - success = pdiffs + 100; - break; - } - - //use needleman to align first barcode.length()+numdiffs of sequence to each barcode - alignment->align(oligo, rawSequence.substr(0,oligo.length()+pdiffs)); - oligo = alignment->getSeqAAln(); - string temp = alignment->getSeqBAln(); - - int alnLength = oligo.length(); - - for(int i=oligo.length()-1;i>=0;i--){ - if(oligo[i] != '-'){ alnLength = i+1; break; } - } - oligo = oligo.substr(0,alnLength); - temp = temp.substr(0,alnLength); - - int numDiff = countDiffs(oligo, temp); - - if(numDiff < minDiff){ - minDiff = numDiff; - minCount = 1; - minGroup = it->second; - minPos = 0; - for(int i=0;i pdiffs) { success = minDiff; } //no good matches - else if(minCount > 1) { success = pdiffs + 10; } //can't tell the difference between multiple primers - else{ //use the best match - group = minGroup; - seq.setUnaligned(rawSequence.substr(minPos)); - success = minDiff; - } - - if (alignment != NULL) { delete alignment; } - - } - - return success; - - } - catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "stripForward"); - exit(1); - } -} - -//*************************************************************************************************************** - -bool TrimFlowsCommand::stripReverse(Sequence& seq){ - try { - - string rawSequence = seq.getUnaligned(); - bool success = 0; //guilty until proven innocent - - for(int i=0;ierrorOut(e, "TrimFlowsCommand", "stripReverse"); - exit(1); - } -} - - -//*************************************************************************************************************** - -bool TrimFlowsCommand::compareDNASeq(string oligo, string seq){ +//********************************************************************/ +string TrimFlowsCommand::reverseOligo(string oligo){ try { - bool success = 1; - int length = oligo.length(); - - for(int i=0;i=0;i--){ + + if(oligo[i] == 'A') { reverse += 'T'; } + else if(oligo[i] == 'T'){ reverse += 'A'; } + else if(oligo[i] == 'U'){ reverse += 'A'; } + + else if(oligo[i] == 'G'){ reverse += 'C'; } + else if(oligo[i] == 'C'){ reverse += 'G'; } + + else if(oligo[i] == 'R'){ reverse += 'Y'; } + else if(oligo[i] == 'Y'){ reverse += 'R'; } + + else if(oligo[i] == 'M'){ reverse += 'K'; } + else if(oligo[i] == 'K'){ reverse += 'M'; } + + else if(oligo[i] == 'W'){ reverse += 'W'; } + else if(oligo[i] == 'S'){ reverse += 'S'; } + + else if(oligo[i] == 'B'){ reverse += 'V'; } + else if(oligo[i] == 'V'){ reverse += 'B'; } + + else if(oligo[i] == 'D'){ reverse += 'H'; } + else if(oligo[i] == 'H'){ reverse += 'D'; } + + else { reverse += 'N'; } + } + + + return reverse; + } catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "compareDNASeq"); + m->errorOut(e, "TrimFlowsCommand", "reverseOligo"); exit(1); } - -} - -//*************************************************************************************************************** - -int TrimFlowsCommand::countDiffs(string oligo, string seq){ - try { - - int length = oligo.length(); - int countDiffs = 0; - - for(int i=0;ierrorOut(e, "TrimFlowsCommand", "countDiffs"); - exit(1); - } - } /**************************************************************************************************/ - -vector TrimFlowsCommand::getFlowFileBreaks() { +vector TrimFlowsCommand::getFlowFileBreaks() { try{ - vector filePos; + vector filePos; filePos.push_back(0); FILE * pFile; - unsigned long int size; + unsigned long long size; //get num bytes in file pFile = fopen (flowFileName.c_str(),"rb"); @@ -968,7 +759,7 @@ vector TrimFlowsCommand::getFlowFileBreaks() { } //estimate file breaks - unsigned long int chunkSize = 0; + unsigned long long chunkSize = 0; chunkSize = size / processors; //file too small to divide by processors @@ -976,7 +767,7 @@ vector TrimFlowsCommand::getFlowFileBreaks() { //for each process seekg to closest file break and search for next '>' char. make that the filebreak for (int i = 0; i < processors; i++) { - unsigned long int spot = (i+1) * chunkSize; + unsigned long long spot = (i+1) * chunkSize; ifstream in; m->openInputFile(flowFileName, in); @@ -985,7 +776,7 @@ vector TrimFlowsCommand::getFlowFileBreaks() { string dummy = m->getline(in); //there was not another sequence before the end of the file - unsigned long int sanityPos = in.tellg(); + unsigned long long sanityPos = in.tellg(); // if (sanityPos == -1) { break; } // else { filePos.push_back(newSpot); } @@ -1007,8 +798,8 @@ vector TrimFlowsCommand::getFlowFileBreaks() { m->openInputFile(flowFileName, in); in >> numFlows; m->gobble(in); - unsigned long int spot = in.tellg(); - filePos[0] = spot; + //unsigned long long spot = in.tellg(); + //filePos[0] = spot; in.close(); processors = (filePos.size() - 1); @@ -1026,10 +817,11 @@ vector TrimFlowsCommand::getFlowFileBreaks() { int TrimFlowsCommand::createProcessesCreateTrim(string flowFileName, string trimFlowFileName, string scrapFlowFileName, string fastaFileName, vector > barcodePrimerComboFileNames){ try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 1; - int exitCommand = 1; processIDS.clear(); + int exitCommand = 1; + +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + int process = 1; //loop through and create all the processes you want while (process != processors) { @@ -1044,11 +836,12 @@ int TrimFlowsCommand::createProcessesCreateTrim(string flowFileName, string trim if(allFiles){ for(int i=0;iopenOutputFile(tempBarcodePrimerComboFileNames[i][j], temp); - temp.close(); - + if (tempBarcodePrimerComboFileNames[i][j] != "") { + tempBarcodePrimerComboFileNames[i][j] += toString(getpid()) + ".temp"; + ofstream temp; + m->openOutputFile(tempBarcodePrimerComboFileNames[i][j], temp); + temp.close(); + } } } } @@ -1086,7 +879,88 @@ int TrimFlowsCommand::createProcessesCreateTrim(string flowFileName, string trim int temp = processIDS[i]; wait(&temp); } +#else + ////////////////////////////////////////////////////////////////////////////////////////////////////// + //Windows version shared memory, so be careful when passing variables through the trimFlowData struct. + //Above fork() will clone, so memory is separate, but that's not the case with windows, + ////////////////////////////////////////////////////////////////////////////////////////////////////// + + vector pDataArray; + DWORD dwThreadIdArray[processors-1]; + HANDLE hThreadArray[processors-1]; + + //Create processor worker threads. + for( int i=0; i > tempBarcodePrimerComboFileNames = barcodePrimerComboFileNames; + if(allFiles){ + for(int i=0;iopenOutputFile(tempBarcodePrimerComboFileNames[i][j], temp); + temp.close(); + } + } + } + } + + trimFlowData* tempflow = new trimFlowData(flowFileName, (trimFlowFileName + extension), (scrapFlowFileName + extension), fastaFileName, flowOrder, tempBarcodePrimerComboFileNames, barcodes, primers, revPrimer, fasta, allFiles, lines[i]->start, lines[i]->end, m, signal, noise, numFlows, maxFlows, minFlows, maxHomoP, tdiffs, bdiffs, pdiffs, i); + pDataArray.push_back(tempflow); + + //MyTrimFlowThreadFunction 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, MyTrimFlowThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]); + } + //using the main process as a worker saves time and memory + ofstream temp; + m->openOutputFile(trimFlowFileName, temp); + temp.close(); + + m->openOutputFile(scrapFlowFileName, temp); + temp.close(); + + if(fasta){ + m->openOutputFile(fastaFileName, temp); + temp.close(); + } + + vector > tempBarcodePrimerComboFileNames = barcodePrimerComboFileNames; + if(allFiles){ + for(int i=0;iopenOutputFile(tempBarcodePrimerComboFileNames[i][j], temp); + temp.close(); + } + + } + } + } + + //do my part - do last piece because windows is looking for eof + int num = driverCreateTrim(flowFileName, (trimFlowFileName + toString(processors-1) + ".temp"), (scrapFlowFileName + toString(processors-1) + ".temp"), (fastaFileName + toString(processors-1) + ".temp"), tempBarcodePrimerComboFileNames, lines[processors-1]); + processIDS.push_back((processors-1)); + + //Wait until all threads have terminated. + 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; + CloseHandle(hThreadArray[i]); + delete pDataArray[i]; + } + + +#endif //append files m->mothurOutEndLine(); for(int i=0;imothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine(); m->appendFiles((trimFlowFileName + toString(processIDS[i]) + ".temp"), trimFlowFileName); - remove((trimFlowFileName + toString(processIDS[i]) + ".temp").c_str()); + m->mothurRemove((trimFlowFileName + toString(processIDS[i]) + ".temp")); // m->mothurOut("\tDone with trim.flow file"); m->mothurOutEndLine(); m->appendFiles((scrapFlowFileName + toString(processIDS[i]) + ".temp"), scrapFlowFileName); - remove((scrapFlowFileName + toString(processIDS[i]) + ".temp").c_str()); + m->mothurRemove((scrapFlowFileName + toString(processIDS[i]) + ".temp")); // m->mothurOut("\tDone with scrap.flow file"); m->mothurOutEndLine(); if(fasta){ m->appendFiles((fastaFileName + toString(processIDS[i]) + ".temp"), fastaFileName); - remove((fastaFileName + toString(processIDS[i]) + ".temp").c_str()); + m->mothurRemove((fastaFileName + toString(processIDS[i]) + ".temp")); // m->mothurOut("\tDone with flow.fasta file"); m->mothurOutEndLine(); } if(allFiles){ for (int j = 0; j < barcodePrimerComboFileNames.size(); j++) { for (int k = 0; k < barcodePrimerComboFileNames[0].size(); k++) { - m->appendFiles((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp"), barcodePrimerComboFileNames[j][k]); - remove((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp").c_str()); + if (barcodePrimerComboFileNames[j][k] != "") { + m->appendFiles((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp"), barcodePrimerComboFileNames[j][k]); + m->mothurRemove((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp")); + } } } } } return exitCommand; -#endif + } catch(exception& e) { m->errorOut(e, "TrimFlowsCommand", "createProcessesCreateTrim");