X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=shhhercommand.cpp;h=4a6e5ff2375ce40560db1aa706edbdacefac0e0d;hp=28e6615f3dffba8fb88c88cd9ca76e79e66faa06;hb=615301e57c25e241356a9c2380648d117709458d;hpb=43ed0accfbc2852849e104ff7eccdd2c42acd4ec diff --git a/shhhercommand.cpp b/shhhercommand.cpp index 28e6615..4a6e5ff 100644 --- a/shhhercommand.cpp +++ b/shhhercommand.cpp @@ -12,17 +12,17 @@ //********************************************************************************************************************** vector ShhherCommand::setParameters(){ try { - CommandParameter pflow("flow", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pflow); - CommandParameter pfile("file", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pfile); - CommandParameter plookup("lookup", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(plookup); - CommandParameter pcutoff("cutoff", "Number", "", "0.01", "", "", "",false,false); parameters.push_back(pcutoff); - CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); - CommandParameter pmaxiter("maxiter", "Number", "", "1000", "", "", "",false,false); parameters.push_back(pmaxiter); - CommandParameter psigma("sigma", "Number", "", "60", "", "", "",false,false); parameters.push_back(psigma); - CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "",false,false); parameters.push_back(pmindelta); - CommandParameter porder("order", "String", "", "", "", "", "",false,false); parameters.push_back(porder); - CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); - CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); + CommandParameter pflow("flow", "InputTypes", "", "", "none", "fileflow", "none","fasta-name-group-counts-qfile",false,false,true); parameters.push_back(pflow); + CommandParameter pfile("file", "InputTypes", "", "", "none", "fileflow", "none","fasta-name-group-counts-qfile",false,false,true); parameters.push_back(pfile); + CommandParameter plookup("lookup", "InputTypes", "", "", "none", "none", "none","",false,false,true); parameters.push_back(plookup); + CommandParameter pcutoff("cutoff", "Number", "", "0.01", "", "", "","",false,false); parameters.push_back(pcutoff); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors); + CommandParameter pmaxiter("maxiter", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(pmaxiter); + CommandParameter plarge("large", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(plarge); + CommandParameter psigma("sigma", "Number", "", "60", "", "", "","",false,false); parameters.push_back(psigma); + CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "","",false,false); parameters.push_back(pmindelta); + CommandParameter porder("order", "Multiple", "A-B-I", "A", "", "", "","",false,false, true); parameters.push_back(porder); 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); } @@ -38,6 +38,11 @@ string ShhherCommand::getHelpString(){ try { string helpString = ""; helpString += "The shhh.flows command reads a file containing flowgrams and creates a file of corrected sequences.\n"; + helpString += "The shhh.flows command parameters are flow, file, lookup, cutoff, processors, large, maxiter, sigma, mindelta and order.\n"; + helpString += "The flow parameter is used to input your flow file.\n"; + helpString += "The file parameter is used to input the *flow.files file created by trim.flows.\n"; + helpString += "The lookup parameter is used specify the lookup file you would like to use. http://www.mothur.org/wiki/Lookup_files.\n"; + helpString += "The order parameter options are A, B or I. Default=A. A = TACG and B = TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC and I = TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC.\n"; return helpString; } catch(exception& e) { @@ -46,6 +51,25 @@ string ShhherCommand::getHelpString(){ } } //********************************************************************************************************************** +string ShhherCommand::getOutputPattern(string type) { + try { + string pattern = ""; + + if (type == "fasta") { pattern = "[filename],shhh.fasta"; } + else if (type == "name") { pattern = "[filename],shhh.names"; } + else if (type == "group") { pattern = "[filename],shhh.groups"; } + else if (type == "counts") { pattern = "[filename],shhh.counts"; } + else if (type == "qfile") { pattern = "[filename],shhh.qual"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } + + return pattern; + } + catch(exception& e) { + m->errorOut(e, "ShhherCommand", "getOutputPattern"); + exit(1); + } +} +//********************************************************************************************************************** ShhherCommand::ShhherCommand(){ try { @@ -53,8 +77,12 @@ ShhherCommand::ShhherCommand(){ setParameters(); //initialize outputTypes -// vector tempOutNames; -// outputTypes["pn.dist"] = tempOutNames; + vector tempOutNames; + outputTypes["fasta"] = tempOutNames; + outputTypes["name"] = tempOutNames; + outputTypes["group"] = tempOutNames; + outputTypes["counts"] = tempOutNames; + outputTypes["qfile"] = tempOutNames; } catch(exception& e) { @@ -67,6 +95,13 @@ ShhherCommand::ShhherCommand(){ ShhherCommand::ShhherCommand(string option) { try { + +#ifdef USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + MPI_Comm_size(MPI_COMM_WORLD, &ncpus); + + if(pid == 0){ +#endif abort = false; calledHelp = false; //allow user to run help @@ -88,9 +123,13 @@ ShhherCommand::ShhherCommand(string option) { } //initialize outputTypes - vector tempOutNames; -// outputTypes["pn.dist"] = tempOutNames; - // outputTypes["fasta"] = tempOutNames; + vector tempOutNames; + outputTypes["fasta"] = tempOutNames; + outputTypes["name"] = tempOutNames; + outputTypes["group"] = tempOutNames; + outputTypes["counts"] = tempOutNames; + outputTypes["qfile"] = 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); @@ -121,13 +160,15 @@ ShhherCommand::ShhherCommand(string option) { if (path == "") { parameters["file"] = inputDir + it->second; } } } - - + + //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 = ""; } + //check for required parameters flowFileName = validParameter.validFile(parameters, "flow", true); flowFilesFileName = validParameter.validFile(parameters, "file", true); if (flowFileName == "not found" && flowFilesFileName == "not found") { - m->mothurOut("values for either flow or file must be provided for the shhh.seqs command."); + m->mothurOut("values for either flow or file must be provided for the shhh.flows command."); m->mothurOutEndLine(); abort = true; } @@ -139,13 +180,22 @@ ShhherCommand::ShhherCommand(string option) { } else{ ofstream temp; - - //flow.files = 9 character offset - compositeFASTAFileName = flowFilesFileName.substr(0, flowFilesFileName.length()-10) + "shhh.fasta"; + + string thisoutputDir = outputDir; + if (outputDir == "") { thisoutputDir = m->hasPath(flowFilesFileName); } //if user entered a file with a path then preserve it + + //we want to rip off .files, and also .flow if its there + string fileroot = m->getRootName(m->getSimpleName(flowFilesFileName)); + if (fileroot[fileroot.length()-1] == '.') { fileroot = fileroot.substr(0, fileroot.length()-1); } //rip off dot + string extension = m->getExtension(fileroot); + if (extension == ".flow") { fileroot = m->getRootName(fileroot); } + else { fileroot += "."; } //add back if needed + + compositeFASTAFileName = thisoutputDir + fileroot + "shhh.fasta"; m->openOutputFile(compositeFASTAFileName, temp); temp.close(); - compositeNamesFileName = flowFilesFileName.substr(0, flowFilesFileName.length()-10) + "shhh.names"; + compositeNamesFileName = thisoutputDir + fileroot + "shhh.names"; m->openOutputFile(compositeNamesFileName, temp); temp.close(); } @@ -207,23 +257,26 @@ ShhherCommand::ShhherCommand(string option) { if (flowFileVector.size() == 0) { m->mothurOut("[ERROR]: no valid files."); m->mothurOutEndLine(); abort = true; } } else{ + if (outputDir == "") { outputDir = m->hasPath(flowFileName); } flowFileVector.push_back(flowFileName); } - - - //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 += m->hasPath(flowFileName); //if user entered a file with a path then preserve it - } - - + //check for optional parameter and set defaults // ...at some point should added some additional type checking... string temp; temp = validParameter.validFile(parameters, "lookup", true); if (temp == "not found") { - lookupFileName = "LookUp_Titanium.pat"; + string path = m->argv; + string tempPath = path; + for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); } + path = path.substr(0, (tempPath.find_last_of('m'))); + +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + path += "lookupFiles/"; +#else + path += "lookupFiles\\"; +#endif + lookupFileName = m->getFullPathName(path) + "LookUp_Titanium.pat"; int ableToOpen; ifstream in; @@ -233,7 +286,7 @@ ShhherCommand::ShhherCommand(string option) { //if you can't open it, try input location if (ableToOpen == 1) { if (inputDir != "") { //default path is set - string tryPath = inputDir + lookupFileName; + string tryPath = inputDir + m->getSimpleName(lookupFileName); m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine(); ifstream in2; ableToOpen = m->openInputFile(tryPath, in2, "noerror"); @@ -281,7 +334,7 @@ ShhherCommand::ShhherCommand(string option) { for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); } exepath = exepath.substr(0, (tempPath.find_last_of('m'))); - string tryPath = m->getFullPathName(exepath) + lookupFileName; + string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName); m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine(); ifstream in2; int ableToOpen = m->openInputFile(tryPath, in2, "noerror"); @@ -303,17 +356,40 @@ ShhherCommand::ShhherCommand(string option) { temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; } m->mothurConvert(temp, maxIters); + + temp = validParameter.validFile(parameters, "large", false); if (temp == "not found"){ temp = "0"; } + m->mothurConvert(temp, largeSize); + if (largeSize != 0) { large = true; } + else { large = false; } + if (largeSize < 0) { m->mothurOut("The value of the large cannot be negative.\n"); } + +#ifdef USE_MPI + if (large) { m->mothurOut("The large parameter is not available with the MPI-Enabled version.\n"); large=false; } +#endif + temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; } m->mothurConvert(temp, sigma); - 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"); - } + 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; + } + } + } +#ifdef USE_MPI + } +#endif } catch(exception& e) { m->errorOut(e, "ShhherCommand", "ShhherCommand"); @@ -328,9 +404,7 @@ int ShhherCommand::execute(){ int tag = 1976; MPI_Status status; - MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are - MPI_Comm_size(MPI_COMM_WORLD, &ncpus); - + if(pid == 0){ for(int i=1;icontrol_pressed) { return 0; } getJointLookUp(); if (m->control_pressed) { return 0; } + vector flowFileVector; + if(flowFilesFileName != "not found"){ + string fName; + + ifstream flowFilesFile; + m->openInputFile(flowFilesFileName, flowFilesFile); + while(flowFilesFile){ + fName = m->getline(flowFilesFile); + flowFileVector.push_back(fName); + m->gobble(flowFilesFile); + } + } + else{ + flowFileVector.push_back(flowFileName); + } + int numFiles = flowFileVector.size(); for(int i=1;icontrol_pressed) { break; } + + getOTUData(listFileName); m->mothurRemove(distFileName); @@ -422,6 +514,7 @@ int ShhherCommand::execute(){ if (m->control_pressed) { break; } + for(int i=1;imothurOutEndLine(); @@ -719,6 +812,37 @@ int ShhherCommand::execute(){ } } /**************************************************************************************************/ +string ShhherCommand::createNamesFile(){ + try{ + + vector duplicateNames(numUniques, ""); + for(int i=0;i variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName)); + string nameFileName = getOutputFileName("name",variables); + + ofstream nameFile; + m->openOutputFile(nameFileName, nameFile); + + for(int i=0;icontrol_pressed) { break; } + + // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl; + nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl; + } + + nameFile.close(); + return nameFileName; + } + catch(exception& e) { + m->errorOut(e, "ShhherCommand", "createNamesFile"); + exit(1); + } +} +/**************************************************************************************************/ string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){ try{ @@ -746,7 +870,7 @@ string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){ } } if(i % 100 == 0){ - m->mothurOut(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n'); + m->mothurOutJustToScreen(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n'); } } @@ -755,7 +879,7 @@ string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){ if (m->control_pressed) { return fDistFileName; } - m->mothurOut(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n'); + m->mothurOutJustToScreen(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n'); ofstream distFile(fDistFileName.c_str()); distFile << outStream.str(); @@ -865,6 +989,8 @@ void ShhherCommand::initPyroCluster(){ try{ if (numOTUs < processors) { processors = 1; } + if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); } + dist.assign(numSeqs * numOTUs, 0); change.assign(numOTUs, 1); centroids.assign(numOTUs, -1); @@ -874,6 +1000,8 @@ void ShhherCommand::initPyroCluster(){ nSeqsBreaks.assign(processors+1, 0); nOTUsBreaks.assign(processors+1, 0); + if (m->debug) { m->mothurOut("[DEBUG]: made it through the memory allocation.\n"); } + nSeqsBreaks[0] = 0; for(int i=0;i> numFlowCells; + string numFlowTest; + flowFile >> numFlowTest; + + if (!m->isContainingOnlyDigits(numFlowTest)) { m->mothurOut("[ERROR]: expected a number and got " + numFlowTest + ", quitting. Did you use the flow parameter instead of the file parameter?"); m->mothurOutEndLine(); exit(1); } + else { convert(numFlowTest, numFlowCells); } + int index = 0;//pcluster while(!flowFile.eof()){ @@ -1267,17 +1400,17 @@ string ShhherCommand::cluster(string distFileName, string namesFileName){ try { ReadMatrix* read = new ReadColumnMatrix(distFileName); - read->setCutoff(cutoff); - - NameAssignment* clusterNameMap = new NameAssignment(namesFileName); - clusterNameMap->readMap(); - read->read(clusterNameMap); - - ListVector* list = read->getListVector(); - SparseMatrix* matrix = read->getMatrix(); + read->setCutoff(cutoff); + + NameAssignment* clusterNameMap = new NameAssignment(namesFileName); + clusterNameMap->readMap(); + read->read(clusterNameMap); - delete read; - delete clusterNameMap; + ListVector* list = read->getListVector(); + SparseDistanceMatrix* matrix = read->getDMatrix(); + + delete read; + delete clusterNameMap; RAbundVector* rabund = new RAbundVector(list->getRAbundVector()); @@ -1537,7 +1670,9 @@ void ShhherCommand::writeQualities(vector otuCounts){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); } - string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.qual"; + map variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName)); + string qualityFileName = getOutputFileName("qfile",variables); ofstream qualityFile; m->openOutputFile(qualityFileName, qualityFile); @@ -1629,7 +1764,7 @@ void ShhherCommand::writeQualities(vector otuCounts){ } } qualityFile.close(); - outputNames.push_back(qualityFileName); + outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName); } catch(exception& e) { @@ -1644,7 +1779,9 @@ void ShhherCommand::writeSequences(vector otuCounts){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); } - string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.fasta"; + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)); + string fastaFileName = getOutputFileName("fasta",variables); ofstream fastaFile; m->openOutputFile(fastaFileName, fastaFile); @@ -1663,7 +1800,7 @@ void ShhherCommand::writeSequences(vector otuCounts){ for(int j=0;j otuCounts){ } fastaFile.close(); - outputNames.push_back(fastaFileName); + outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName); if(compositeFASTAFileName != ""){ m->appendFiles(fastaFileName, compositeFASTAFileName); @@ -1692,7 +1829,9 @@ void ShhherCommand::writeNames(vector otuCounts){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); } - string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.names"; + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)); + string nameFileName = getOutputFileName("name",variables); ofstream nameFile; m->openOutputFile(nameFileName, nameFile); @@ -1711,7 +1850,7 @@ void ShhherCommand::writeNames(vector otuCounts){ } } nameFile.close(); - outputNames.push_back(nameFileName); + outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName); if(compositeNamesFileName != ""){ @@ -1730,17 +1869,22 @@ void ShhherCommand::writeGroups(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); } - string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)); - string groupFileName = fileRoot + "shhh.groups"; + string fileRoot = m->getRootName(m->getSimpleName(flowFileName)); + int pos = fileRoot.find_first_of('.'); + string fileGroup = fileRoot; + if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); } + map variables; + variables["[filename]"] = thisOutputDir + fileRoot; + string groupFileName = getOutputFileName("group",variables); ofstream groupFile; m->openOutputFile(groupFileName, groupFile); for(int i=0;icontrol_pressed) { break; } - groupFile << seqNameVector[i] << '\t' << fileRoot << endl; + groupFile << seqNameVector[i] << '\t' << fileGroup << endl; } groupFile.close(); - outputNames.push_back(groupFileName); + outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName); } catch(exception& e) { @@ -1755,7 +1899,9 @@ void ShhherCommand::writeClusters(vector otuCounts){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); } - string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.counts"; + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)); + string otuCountsFileName = getOutputFileName("counts",variables); ofstream otuCountsFile; m->openOutputFile(otuCountsFileName, otuCountsFile); @@ -1772,7 +1918,7 @@ void ShhherCommand::writeClusters(vector otuCounts){ otuCountsFile << "ideal\t"; for(int j=8;j otuCounts){ string newSeq = ""; for(int k=0;k otuCounts){ } } otuCountsFile.close(); - outputNames.push_back(otuCountsFileName); + outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName); } catch(exception& e) { @@ -1814,7 +1960,7 @@ void ShhherCommand::writeClusters(vector otuCounts){ int ShhherCommand::execute(){ try { - if (abort == true) { return 0; } + if (abort == true) { if (calledHelp) { return 0; } return 2; } getSingleLookUp(); if (m->control_pressed) { return 0; } getJointLookUp(); if (m->control_pressed) { return 0; } @@ -1824,15 +1970,15 @@ int ShhherCommand::execute(){ if (numFiles < processors) { processors = numFiles; } #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) - if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size()); } + if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName); } else { createProcesses(flowFileVector); } //each processor processes one file #else - driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size()); + driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName); #endif if(compositeFASTAFileName != ""){ - outputNames.push_back(compositeFASTAFileName); - outputNames.push_back(compositeNamesFileName); + outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName); + outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName); } m->mothurOutEndLine(); @@ -1848,6 +1994,40 @@ int ShhherCommand::execute(){ } } #endif +//******************************************************************************************************************** +//sorts biggest to smallest +inline bool compareFileSizes(string left, string right){ + + FILE * pFile; + long leftsize = 0; + + //get num bytes in file + string filename = left; + pFile = fopen (filename.c_str(),"rb"); + string error = "Error opening " + filename; + if (pFile==NULL) perror (error.c_str()); + else{ + fseek (pFile, 0, SEEK_END); + leftsize=ftell (pFile); + fclose (pFile); + } + + FILE * pFile2; + long rightsize = 0; + + //get num bytes in file + filename = right; + pFile2 = fopen (filename.c_str(),"rb"); + error = "Error opening " + filename; + if (pFile2==NULL) perror (error.c_str()); + else{ + fseek (pFile2, 0, SEEK_END); + rightsize=ftell (pFile2); + fclose (pFile2); + } + + return (leftsize > rightsize); +} /**************************************************************************************************/ int ShhherCommand::createProcesses(vector filenames){ @@ -1858,16 +2038,40 @@ int ShhherCommand::createProcesses(vector filenames){ //sanity check if (filenames.size() < processors) { processors = filenames.size(); } - + + //sort file names by size to divide load better + sort(filenames.begin(), filenames.end(), compareFileSizes); + + vector < vector > dividedFiles; //dividedFiles[1] = vector of filenames for process 1... + dividedFiles.resize(processors); + + //for each file, figure out which process will complete it + //want to divide the load intelligently so the big files are spread between processes + for (int i = 0; i < filenames.size(); i++) { + int processToAssign = (i+1) % processors; + if (processToAssign == 0) { processToAssign = processors; } + + dividedFiles[(processToAssign-1)].push_back(filenames[i]); + } + + //now lets reverse the order of ever other process, so we balance big files running with little ones + for (int i = 0; i < processors; i++) { + int remainder = ((i+1) % processors); + if (remainder) { reverse(dividedFiles[i].begin(), dividedFiles[i].end()); } + } + + //divide the groups between the processors - vector lines; + /*vector lines; + vector numFilesToComplete; int numFilesPerProcessor = filenames.size() / processors; for (int i = 0; i < processors; i++) { int startIndex = i * numFilesPerProcessor; int endIndex = (i+1) * numFilesPerProcessor; if(i == (processors - 1)){ endIndex = filenames.size(); } lines.push_back(linePair(startIndex, endIndex)); - } + numFilesToComplete.push_back((endIndex-startIndex)); + }*/ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) @@ -1879,7 +2083,15 @@ int ShhherCommand::createProcesses(vector filenames){ processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later process++; }else if (pid == 0){ - num = driver(filenames, compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp", lines[process].start, lines[process].end); + num = driver(dividedFiles[process], compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp"); + + //pass numSeqs to parent + ofstream out; + string tempFile = compositeFASTAFileName + toString(getpid()) + ".num.temp"; + m->openOutputFile(tempFile, out); + out << num << endl; + out.close(); + exit(0); }else { m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); @@ -1889,7 +2101,7 @@ int ShhherCommand::createProcesses(vector filenames){ } //do my part - driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[0].start, lines[0].end); + driver(dividedFiles[0], compositeFASTAFileName, compositeNamesFileName); //force parent to wait until all the processes are done for (int i=0;i filenames){ //Windows version shared memory, so be careful when passing variables through the shhhFlowsData 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]; @@ -1938,10 +2150,22 @@ int ShhherCommand::createProcesses(vector filenames){ CloseHandle(hThreadArray[i]); delete pDataArray[i]; } - + */ #endif for (int i=0;iopenInputFile(tempFile, in); + if (!in.eof()) { + int tempNum = 0; + in >> tempNum; + if (tempNum != dividedFiles[i+1].size()) { + m->mothurOut("[ERROR]: main process expected " + toString(processIDS[i]) + " to complete " + toString(dividedFiles[i+1].size()) + " files, and it only reported completing " + toString(tempNum) + ". This will cause file mismatches. The flow files may be too large to process with multiple processors. \n"); + } + } + in.close(); m->mothurRemove(tempFile); + if (compositeFASTAFileName != "") { m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName); m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName); @@ -1960,176 +2184,276 @@ int ShhherCommand::createProcesses(vector filenames){ } /**************************************************************************************************/ -int ShhherCommand::driver(vector filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){ +vector ShhherCommand::parseFlowFiles(string filename){ try { + vector files; + int count = 0; - for(int i=start;icontrol_pressed) { break; } - - string flowFileName = filenames[i]; - - m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n"); - m->mothurOut("Reading flowgrams...\n"); - - vector seqNameVector; - vector lengths; - vector flowDataIntI; - vector flowDataPrI; - map nameMap; - vector uniqueFlowgrams; - vector uniqueCount; - vector mapSeqToUnique; - vector mapUniqueToSeq; - vector uniqueLengths; - int numFlowCells; - - int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells); - - if (m->control_pressed) { break; } - - m->mothurOut("Identifying unique flowgrams...\n"); - int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI); - - if (m->control_pressed) { break; } - - m->mothurOut("Calculating distances between flowgrams...\n"); - string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist"; - unsigned long long begTime = time(NULL); - double begClock = clock(); + ifstream in; + m->openInputFile(filename, in); - flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI); - - m->mothurOutEndLine(); - m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n'); - + int thisNumFLows = 0; + in >> thisNumFLows; m->gobble(in); + + while (!in.eof()) { + if (m->control_pressed) { break; } - string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names"; - createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq); - - if (m->control_pressed) { break; } - - m->mothurOut("\nClustering flowgrams...\n"); - string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list"; - cluster(listFileName, distFileName, namesFileName); - - if (m->control_pressed) { break; } + ofstream out; + string outputFileName = filename + toString(count) + ".temp"; + m->openOutputFile(outputFileName, out); + out << thisNumFLows << endl; + files.push_back(outputFileName); + + int numLinesWrote = 0; + for (int i = 0; i < largeSize; i++) { + if (in.eof()) { break; } + string line = m->getline(in); m->gobble(in); + out << line << endl; + numLinesWrote++; + } + out.close(); - vector otuData; - vector cumNumSeqs; - vector nSeqsPerOTU; - vector > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices - vector > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI - vector seqNumber; //tMaster->anP: the sequence id number sorted by OTU - vector seqIndex; //tMaster->anI; the index that corresponds to seqNumber + if (numLinesWrote == 0) { m->mothurRemove(outputFileName); files.pop_back(); } + count++; + } + in.close(); + + if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } files.clear(); } + + m->mothurOut("\nDivided " + filename + " into " + toString(files.size()) + " files.\n\n"); + + return files; + } + catch(exception& e) { + m->errorOut(e, "ShhherCommand", "parseFlowFiles"); + exit(1); + } +} +/**************************************************************************************************/ - - int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap); +int ShhherCommand::driver(vector filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName){ + try { + + int numCompleted = 0; + + for(int i=0;icontrol_pressed) { break; } - m->mothurRemove(distFileName); - m->mothurRemove(namesFileName); - m->mothurRemove(listFileName); - - vector dist; //adDist - distance of sequences to centroids - vector change; //did the centroid sequence change? 0 = no; 1 = yes - vector centroids; //the representative flowgram for each cluster m - vector weight; - vector singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs) - vector nSeqsBreaks; - vector nOTUsBreaks; + vector theseFlowFileNames; theseFlowFileNames.push_back(filenames[i]); + if (large) { theseFlowFileNames = parseFlowFiles(filenames[i]); } - dist.assign(numSeqs * numOTUs, 0); - change.assign(numOTUs, 1); - centroids.assign(numOTUs, -1); - weight.assign(numOTUs, 0); - singleTau.assign(numSeqs, 1.0); + if (m->control_pressed) { break; } - nSeqsBreaks.assign(2, 0); - nOTUsBreaks.assign(2, 0); + double begClock = clock(); + unsigned long long begTime; - nSeqsBreaks[0] = 0; - nSeqsBreaks[1] = numSeqs; - nOTUsBreaks[1] = numOTUs; - - if (m->control_pressed) { break; } - - double maxDelta = 0; - int iter = 0; - - begClock = clock(); - begTime = time(NULL); + string fileNameForOutput = filenames[i]; - m->mothurOut("\nDenoising flowgrams...\n"); - m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n"); - - while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){ - - if (m->control_pressed) { break; } - - double cycClock = clock(); - unsigned long long cycTime = time(NULL); - fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); - - if (m->control_pressed) { break; } + for (int g = 0; g < theseFlowFileNames.size(); g++) { - calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber); - - if (m->control_pressed) { break; } + string flowFileName = theseFlowFileNames[g]; + m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n"); + m->mothurOut("Reading flowgrams...\n"); - maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight); + vector seqNameVector; + vector lengths; + vector flowDataIntI; + vector flowDataPrI; + map nameMap; + vector uniqueFlowgrams; + vector uniqueCount; + vector mapSeqToUnique; + vector mapUniqueToSeq; + vector uniqueLengths; + int numFlowCells; + + if (m->debug) { m->mothurOut("[DEBUG]: About to read flowgrams.\n"); } + int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells); if (m->control_pressed) { break; } - double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight); + m->mothurOut("Identifying unique flowgrams...\n"); + int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI); if (m->control_pressed) { break; } - checkCentroids(numOTUs, centroids, weight); - - if (m->control_pressed) { break; } - - calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths); - - if (m->control_pressed) { break; } - - iter++; - - m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n'); + m->mothurOut("Calculating distances between flowgrams...\n"); + string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist"; + begTime = time(NULL); + - } - - if (m->control_pressed) { break; } - - m->mothurOut("\nFinalizing...\n"); - fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); - - if (m->control_pressed) { break; } - - setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI); - - if (m->control_pressed) { break; } - - vector otuCounts(numOTUs, 0); - for(int i=0;icontrol_pressed) { break; } + flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI); + + m->mothurOutEndLine(); + m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n'); + + + string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names"; + createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq); + + if (m->control_pressed) { break; } + + m->mothurOut("\nClustering flowgrams...\n"); + string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list"; + cluster(listFileName, distFileName, namesFileName); + + if (m->control_pressed) { break; } + + vector otuData; + vector cumNumSeqs; + vector nSeqsPerOTU; + vector > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices + vector > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI + vector seqNumber; //tMaster->anP: the sequence id number sorted by OTU + vector seqIndex; //tMaster->anI; the index that corresponds to seqNumber + + + int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap); + + if (m->control_pressed) { break; } + + m->mothurRemove(distFileName); + m->mothurRemove(namesFileName); + m->mothurRemove(listFileName); + + vector dist; //adDist - distance of sequences to centroids + vector change; //did the centroid sequence change? 0 = no; 1 = yes + vector centroids; //the representative flowgram for each cluster m + vector weight; + vector singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs) + vector nSeqsBreaks; + vector nOTUsBreaks; + + if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); } + + dist.assign(numSeqs * numOTUs, 0); + change.assign(numOTUs, 1); + centroids.assign(numOTUs, -1); + weight.assign(numOTUs, 0); + singleTau.assign(numSeqs, 1.0); + + nSeqsBreaks.assign(2, 0); + nOTUsBreaks.assign(2, 0); + + nSeqsBreaks[0] = 0; + nSeqsBreaks[1] = numSeqs; + nOTUsBreaks[1] = numOTUs; + + if (m->debug) { m->mothurOut("[DEBUG]: done allocating memory, about to denoise.\n"); } + + if (m->control_pressed) { break; } + + double maxDelta = 0; + int iter = 0; + + begClock = clock(); + begTime = time(NULL); + + m->mothurOut("\nDenoising flowgrams...\n"); + m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n"); + + while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){ + + if (m->control_pressed) { break; } + + double cycClock = clock(); + unsigned long long cycTime = time(NULL); + fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); + + if (m->control_pressed) { break; } + + calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber); + + if (m->control_pressed) { break; } + + maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight); + + if (m->control_pressed) { break; } + + double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight); + + if (m->control_pressed) { break; } + + checkCentroids(numOTUs, centroids, weight); + + if (m->control_pressed) { break; } + + calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths); + + if (m->control_pressed) { break; } + + iter++; + + m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n'); + + } + + if (m->control_pressed) { break; } + + m->mothurOut("\nFinalizing...\n"); + fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); + + if (m->control_pressed) { break; } + + setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI); + + if (m->control_pressed) { break; } + + vector otuCounts(numOTUs, 0); + for(int j=0;jcontrol_pressed) { break; } + + if ((large) && (g == 0)) { flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; } + string thisOutputDir = outputDir; + if (outputDir == "") { thisOutputDir = m->hasPath(flowFileName); } + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)); + string qualityFileName = getOutputFileName("qfile",variables); + string fastaFileName = getOutputFileName("fasta",variables); + string nameFileName = getOutputFileName("name",variables); + string otuCountsFileName = getOutputFileName("counts",variables); + string fileRoot = m->getRootName(m->getSimpleName(flowFileName)); + int pos = fileRoot.find_first_of('.'); + string fileGroup = fileRoot; + if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); } + string groupFileName = getOutputFileName("group",variables); + + + writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; } + writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; } + writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; } + writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; } + writeGroups(groupFileName, fileGroup, numSeqs, seqNameVector); if (m->control_pressed) { break; } + + if (large) { + if (g > 0) { + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])); + m->appendFiles(qualityFileName, getOutputFileName("qfile",variables)); + m->mothurRemove(qualityFileName); + m->appendFiles(fastaFileName, getOutputFileName("fasta",variables)); + m->mothurRemove(fastaFileName); + m->appendFiles(nameFileName, getOutputFileName("name",variables)); + m->mothurRemove(nameFileName); + m->appendFiles(otuCountsFileName, getOutputFileName("counts",variables)); + m->mothurRemove(otuCountsFileName); + m->appendFiles(groupFileName, getOutputFileName("group",variables)); + m->mothurRemove(groupFileName); + } + m->mothurRemove(theseFlowFileNames[g]); + } + } - writeQualities(numOTUs, numFlowCells, flowFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; } - writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, flowFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; } - writeNames(thisCompositeNamesFileName, numOTUs, flowFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; } - writeClusters(flowFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; } - writeGroups(flowFileName, numSeqs, seqNameVector); if (m->control_pressed) { break; } - - m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n'); + numCompleted++; + m->mothurOut("Total time to process " + fileNameForOutput + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n'); } if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } - return 0; + return numCompleted; }catch(exception& e) { m->errorOut(e, "ShhherCommand", "driver"); @@ -2153,18 +2477,27 @@ int ShhherCommand::getFlowData(string filename, vector& thisSeqNameVecto thisFlowDataIntI.clear(); thisNameMap.clear(); - flowFile >> numFlowCells; + string numFlowTest; + flowFile >> numFlowTest; + + if (!m->isContainingOnlyDigits(numFlowTest)) { m->mothurOut("[ERROR]: expected a number and got " + numFlowTest + ", quitting. Did you use the flow parameter instead of the file parameter?"); m->mothurOutEndLine(); exit(1); } + else { convert(numFlowTest, numFlowCells); } + + if (m->debug) { m->mothurOut("[DEBUG]: numFlowCells = " + toString(numFlowCells) + ".\n"); } int index = 0;//pcluster while(!flowFile.eof()){ if (m->control_pressed) { break; } flowFile >> seqName >> currentNumFlowCells; + thisLengths.push_back(currentNumFlowCells); thisSeqNameVector.push_back(seqName); thisNameMap[seqName] = index++;//pcluster - + + if (m->debug) { m->mothurOut("[DEBUG]: seqName = " + seqName + " length = " + toString(currentNumFlowCells) + " index = " + toString(index) + "\n"); } + for(int i=0;i> intensity; if(intensity > 9.99) { intensity = 9.99; } @@ -2224,9 +2557,8 @@ int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int } } if(i % 100 == 0){ - m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime)); - m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)); - m->mothurOutEndLine(); + m->mothurOutJustToScreen(toString(i) + "\t" + toString(time(NULL) - begTime)); + m->mothurOutJustToScreen("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)+"\n"); } } @@ -2236,9 +2568,8 @@ int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int if (m->control_pressed) {} else { - m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime)); - m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)); - m->mothurOutEndLine(); + m->mothurOutJustToScreen(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime)); + m->mothurOutJustToScreen("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)+"\n"); } return 0; @@ -2397,16 +2728,17 @@ int ShhherCommand::cluster(string filename, string distFileName, string namesFil NameAssignment* clusterNameMap = new NameAssignment(namesFileName); clusterNameMap->readMap(); read->read(clusterNameMap); - + ListVector* list = read->getListVector(); - SparseMatrix* matrix = read->getMatrix(); + SparseDistanceMatrix* matrix = read->getDMatrix(); delete read; delete clusterNameMap; RAbundVector* rabund = new RAbundVector(list->getRAbundVector()); - Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest"); + float adjust = -1.0; + Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest", adjust); string tag = cluster->getTag(); double clusterCutoff = cutoff; @@ -2452,6 +2784,8 @@ int ShhherCommand::getOTUData(int numSeqs, string fileName, vector& otuDat listFile >> label >> numOTUs; + if (m->debug) { m->mothurOut("[DEBUG]: Getting OTU Data...\n"); } + otuData.assign(numSeqs, 0); cumNumSeqs.assign(numOTUs, 0); nSeqsPerOTU.assign(numOTUs, 0); @@ -2466,6 +2800,7 @@ int ShhherCommand::getOTUData(int numSeqs, string fileName, vector& otuDat for(int i=0;icontrol_pressed) { break; } + if (m->debug) { m->mothurOut("[DEBUG]: processing OTU " + toString(i) + ".\n"); } listFile >> singleOTU; @@ -2923,14 +3258,11 @@ void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector& seqNumber, ve } /**************************************************************************************************/ -void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string filename, vector otuCounts, vector& nSeqsPerOTU, vector& seqNumber, +void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string qualityFileName, vector otuCounts, vector& nSeqsPerOTU, vector& seqNumber, vector& singleTau, vector& flowDataIntI, vector& uniqueFlowgrams, vector& cumNumSeqs, vector& mapUniqueToSeq, vector& seqNameVector, vector& centroids, vector >& aaI){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir += m->hasPath(filename); } - string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.qual"; ofstream qualityFile; m->openOutputFile(qualityFileName, qualityFile); @@ -3012,17 +3344,18 @@ void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string filenam if(otuCounts[i] > 0){ qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl; - + int j=4; //need to get past the first four bases while(qualities[i][j] != -1){ - qualityFile << qualities[i][j] << ' '; - j++; + qualityFile << qualities[i][j] << ' '; + if (j > qualities[i].size()) { break; } + j++; } qualityFile << endl; } } qualityFile.close(); - outputNames.push_back(qualityFileName); + outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName); } catch(exception& e) { @@ -3033,11 +3366,9 @@ void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string filenam /**************************************************************************************************/ -void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTUs, int numFlowCells, string filename, vector otuCounts, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& centroids){ +void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTUs, int numFlowCells, string fastaFileName, vector otuCounts, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& centroids){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir += m->hasPath(filename); } - string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.fasta"; + ofstream fastaFile; m->openOutputFile(fastaFileName, fastaFile); @@ -3056,18 +3387,19 @@ void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTU for(int j=0;j= 4) { fastaFile << newSeq.substr(4) << endl; } + else { fastaFile << "NNNN" << endl; } } } fastaFile.close(); - outputNames.push_back(fastaFileName); + outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName); if(thisCompositeFASTAFileName != ""){ m->appendFiles(fastaFileName, thisCompositeFASTAFileName); @@ -3081,11 +3413,9 @@ void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTU /**************************************************************************************************/ -void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string filename, vector otuCounts, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU){ +void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string nameFileName, vector otuCounts, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir += m->hasPath(filename); } - string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.names"; + ofstream nameFile; m->openOutputFile(nameFileName, nameFile); @@ -3104,7 +3434,7 @@ void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, s } } nameFile.close(); - outputNames.push_back(nameFileName); + outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName); if(thisCompositeNamesFileName != ""){ @@ -3119,13 +3449,9 @@ void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, s /**************************************************************************************************/ -void ShhherCommand::writeGroups(string filename, int numSeqs, vector& seqNameVector){ +void ShhherCommand::writeGroups(string groupFileName, string fileRoot, int numSeqs, vector& seqNameVector){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir += m->hasPath(filename); } - string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(filename)); - string groupFileName = fileRoot + "shhh.groups"; - ofstream groupFile; + ofstream groupFile; m->openOutputFile(groupFileName, groupFile); for(int i=0;i& se groupFile << seqNameVector[i] << '\t' << fileRoot << endl; } groupFile.close(); - outputNames.push_back(groupFileName); + outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName); } catch(exception& e) { @@ -3144,11 +3470,8 @@ void ShhherCommand::writeGroups(string filename, int numSeqs, vector& se /**************************************************************************************************/ -void ShhherCommand::writeClusters(string filename, int numOTUs, int numFlowCells, vector otuCounts, vector& centroids, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU, vector& lengths, vector& flowDataIntI){ +void ShhherCommand::writeClusters(string otuCountsFileName, int numOTUs, int numFlowCells, vector otuCounts, vector& centroids, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU, vector& lengths, vector& flowDataIntI){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir += m->hasPath(filename); } - string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.counts"; ofstream otuCountsFile; m->openOutputFile(otuCountsFileName, otuCountsFile); @@ -3165,7 +3488,7 @@ void ShhherCommand::writeClusters(string filename, int numOTUs, int numFlowCells otuCountsFile << "ideal\t"; for(int j=8;j= 4) { otuCountsFile << newSeq.substr(4) << endl; } + else { otuCountsFile << "NNNN" << endl; } } otuCountsFile << endl; } } otuCountsFile.close(); - outputNames.push_back(otuCountsFileName); + outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName); } catch(exception& e) {