X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=chimerauchimecommand.cpp;h=7ff3989dea527060934b417de22ef5078df7a074;hb=035f86272c776e1cccaa47021e26782e49cd41e7;hp=73e7ace9860eeb6410517497b5e0495e9fdb1b96;hpb=c47e480b743d1c242b8c527b6d12f992c68b8c2c;p=mothur.git diff --git a/chimerauchimecommand.cpp b/chimerauchimecommand.cpp index 73e7ace..7ff3989 100644 --- a/chimerauchimecommand.cpp +++ b/chimerauchimecommand.cpp @@ -9,18 +9,19 @@ #include "chimerauchimecommand.h" #include "deconvolutecommand.h" -#include "uc.h" +//#include "uc.h" #include "sequence.hpp" #include "referencedb.h" - +#include "systemcommand.h" //********************************************************************************************************************** vector ChimeraUchimeCommand::setParameters(){ try { CommandParameter ptemplate("reference", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptemplate); CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta); - CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname); - CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup); + CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none",false,false); parameters.push_back(pname); + CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none",false,false); parameters.push_back(pcount); + CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none",false,false); parameters.push_back(pgroup); CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); @@ -58,9 +59,10 @@ string ChimeraUchimeCommand::getHelpString(){ string helpString = ""; helpString += "The chimera.uchime command reads a fastafile and referencefile and outputs potentially chimeric sequences.\n"; helpString += "This command is a wrapper for uchime written by Robert C. Edgar.\n"; - helpString += "The chimera.uchime command parameters are fasta, name, reference, processors, abskew, chimealns, minh, mindiv, xn, dn, xa, chunks, minchunk, idsmoothwindow, minsmoothid, maxp, skipgaps, skipgaps2, minlen, maxlen, ucl and queryfact.\n"; + helpString += "The chimera.uchime command parameters are fasta, name, count, reference, processors, abskew, chimealns, minh, mindiv, xn, dn, xa, chunks, minchunk, idsmoothwindow, minsmoothid, maxp, skipgaps, skipgaps2, minlen, maxlen, ucl and queryfact.\n"; helpString += "The fasta parameter allows you to enter the fasta file containing your potentially chimeric sequences, and is required, unless you have a valid current fasta file. \n"; helpString += "The name parameter allows you to provide a name file, if you are using template=self. \n"; + helpString += "The count parameter allows you to provide a count file, if you are using template=self. \n"; helpString += "You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amazon.fasta \n"; helpString += "The group parameter allows you to provide a group file. The group file can be used with a namesfile and reference=self. When checking sequences, only sequences from the same group as the query sequence will be used as the reference. \n"; helpString += "The reference parameter allows you to enter a reference file containing known non-chimeric sequences, and is required. You may also set template=self, in this case the abundant sequences will be used as potential parents. \n"; @@ -98,6 +100,28 @@ string ChimeraUchimeCommand::getHelpString(){ } } //********************************************************************************************************************** +string ChimeraUchimeCommand::getOutputFileNameTag(string type, string inputName=""){ + try { + string outputFileName = ""; + map >::iterator it; + + //is this a type this command creates + it = outputTypes.find(type); + if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); } + else { + if (type == "chimera") { outputFileName = "uchime.chimeras"; } + else if (type == "accnos") { outputFileName = "uchime.accnos"; } + else if (type == "alns") { outputFileName = "uchime.alns"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; } + } + return outputFileName; + } + catch(exception& e) { + m->errorOut(e, "ChimeraUchimeCommand", "getOutputFileNameTag"); + exit(1); + } +} +//********************************************************************************************************************** ChimeraUchimeCommand::ChimeraUchimeCommand(){ try { abort = true; calledHelp = true; @@ -115,7 +139,7 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(){ //*************************************************************************************************************** ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { try { - abort = false; calledHelp = false; + abort = false; calledHelp = false; hasName=false; hasCount=false; ReferenceDB* rdb = ReferenceDB::getInstance(); //allow user to run help @@ -225,9 +249,8 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { //check for required parameters - bool hasName = true; namefile = validParameter.validFile(parameters, "name", false); - if (namefile == "not found") { namefile = ""; hasName = false; } + if (namefile == "not found") { namefile = ""; } else { m->splitAtDash(namefile, nameFileNames); @@ -294,12 +317,91 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { } } } + } + + if (nameFileNames.size() != 0) { hasName = true; } + + //check for required parameters + vector countfileNames; + countfile = validParameter.validFile(parameters, "count", false); + if (countfile == "not found") { + countfile = ""; + }else { + m->splitAtDash(countfile, countfileNames); - //make sure there is at least one valid file left - if (nameFileNames.size() == 0) { m->mothurOut("[ERROR]: no valid name files."); m->mothurOutEndLine(); abort = true; } + //go through files and make sure they are good, if not, then disregard them + for (int i = 0; i < countfileNames.size(); i++) { + + bool ignore = false; + if (countfileNames[i] == "current") { + countfileNames[i] = m->getCountTableFile(); + if (nameFileNames[i] != "") { m->mothurOut("Using " + countfileNames[i] + " as input file for the count parameter where you had given current."); m->mothurOutEndLine(); } + else { + m->mothurOut("You have no current count file, ignoring current."); m->mothurOutEndLine(); ignore=true; + //erase from file list + countfileNames.erase(countfileNames.begin()+i); + i--; + } + } + + if (!ignore) { + + if (inputDir != "") { + string path = m->hasPath(countfileNames[i]); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { countfileNames[i] = inputDir + countfileNames[i]; } + } + + int ableToOpen; + ifstream in; + + ableToOpen = m->openInputFile(countfileNames[i], in, "noerror"); + + //if you can't open it, try default location + if (ableToOpen == 1) { + if (m->getDefaultPath() != "") { //default path is set + string tryPath = m->getDefaultPath() + m->getSimpleName(countfileNames[i]); + m->mothurOut("Unable to open " + countfileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine(); + ifstream in2; + ableToOpen = m->openInputFile(tryPath, in2, "noerror"); + in2.close(); + countfileNames[i] = tryPath; + } + } + + if (ableToOpen == 1) { + if (m->getOutputDir() != "") { //default path is set + string tryPath = m->getOutputDir() + m->getSimpleName(countfileNames[i]); + m->mothurOut("Unable to open " + countfileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine(); + ifstream in2; + ableToOpen = m->openInputFile(tryPath, in2, "noerror"); + in2.close(); + countfileNames[i] = tryPath; + } + } + + in.close(); + + if (ableToOpen == 1) { + m->mothurOut("Unable to open " + countfileNames[i] + ". It will be disregarded."); m->mothurOutEndLine(); + //erase from file list + countfileNames.erase(countfileNames.begin()+i); + i--; + }else { + m->setCountTableFile(countfileNames[i]); + } + } + } } - - if (hasName && (nameFileNames.size() != fastaFileNames.size())) { m->mothurOut("[ERROR]: The number of namefiles does not match the number of fastafiles, please correct."); m->mothurOutEndLine(); abort=true; } + + if (countfileNames.size() != 0) { hasCount = true; } + + //make sure there is at least one valid file left + if (hasName && hasCount) { m->mothurOut("[ERROR]: You must enter ONLY ONE of the following: count or name."); m->mothurOutEndLine(); abort = true; } + + if (!hasName && hasCount) { nameFileNames = countfileNames; } + + if ((hasCount || hasName) && (nameFileNames.size() != fastaFileNames.size())) { m->mothurOut("[ERROR]: The number of name or count files does not match the number of fastafiles, please correct."); m->mothurOutEndLine(); abort=true; } bool hasGroup = true; groupfile = validParameter.validFile(parameters, "group", false); @@ -377,6 +479,10 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { if (hasGroup && (groupFileNames.size() != fastaFileNames.size())) { m->mothurOut("[ERROR]: The number of groupfiles does not match the number of fastafiles, please correct."); m->mothurOutEndLine(); abort=true; } + if (hasGroup && hasCount) { m->mothurOut("[ERROR]: You must enter ONLY ONE of the following: count or group."); m->mothurOutEndLine(); 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"){ outputDir = ""; } + //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 = ""; } @@ -405,6 +511,7 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { } } }else if (hasName) { templatefile = "self"; } + else if (hasCount) { templatefile = "self"; } else { if (rdb->getSavedReference() != "") { templatefile = rdb->getSavedReference(); @@ -418,7 +525,7 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { string temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } m->setProcessors(temp); - convert(temp, processors); + m->mothurConvert(temp, processors); abskew = validParameter.validFile(parameters, "abskew", false); if (abskew == "not found"){ useAbskew = false; abskew = "1.9"; }else{ useAbskew = true; } if (useAbskew && templatefile != "self") { m->mothurOut("The abskew parameter is only valid with template=self, ignoring."); m->mothurOutEndLine(); useAbskew = false; } @@ -453,8 +560,54 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { if (hasName && (templatefile != "self")) { m->mothurOut("You have provided a namefile and the reference parameter is not set to self. I am not sure what reference you are trying to use, aborting."); m->mothurOutEndLine(); abort=true; } if (hasGroup && (templatefile != "self")) { m->mothurOut("You have provided a group file and the reference parameter is not set to self. I am not sure what reference you are trying to use, aborting."); m->mothurOutEndLine(); abort=true; } + + //look for uchime exe + 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'))); + + string uchimeCommand; +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + uchimeCommand = path + "uchime"; // format the database, -o option gives us the ability + if (m->debug) { + m->mothurOut("[DEBUG]: Uchime location using \"which uchime\" = "); + Command* newCommand = new SystemCommand("which uchime"); m->mothurOutEndLine(); + newCommand->execute(); + delete newCommand; + m->mothurOut("[DEBUG]: Mothur's location using \"which mothur\" = "); + newCommand = new SystemCommand("which mothur"); m->mothurOutEndLine(); + newCommand->execute(); + delete newCommand; + } +#else + uchimeCommand = path + "uchime.exe"; +#endif + + //test to make sure uchime exists + ifstream in; + uchimeCommand = m->getFullPathName(uchimeCommand); + int ableToOpen = m->openInputFile(uchimeCommand, in, "no error"); in.close(); + if(ableToOpen == 1) { + m->mothurOut(uchimeCommand + " file does not exist. Checking path... \n"); + //check to see if uchime is in the path?? + + string uLocation = m->findProgramPath("uchime"); + + + ifstream in2; +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + ableToOpen = m->openInputFile(uLocation, in2, "no error"); in2.close(); +#else + ableToOpen = m->openInputFile((uLocation + ".exe"), in2, "no error"); in2.close(); +#endif - } + if(ableToOpen == 1) { m->mothurOut("[ERROR]: " + uLocation + " file does not exist. mothur requires the uchime executable."); m->mothurOutEndLine(); abort = true; } + else { m->mothurOut("Found uchime in your path, using " + uLocation + "\n");uchimeLocation = uLocation; } + }else { uchimeLocation = uchimeCommand; } + + uchimeLocation = m->getFullPathName(uchimeLocation); + } } catch(exception& e) { m->errorOut(e, "ChimeraSlayerCommand", "ChimeraSlayerCommand"); @@ -465,7 +618,8 @@ ChimeraUchimeCommand::ChimeraUchimeCommand(string option) { int ChimeraUchimeCommand::execute(){ try{ - if (abort == true) { if (calledHelp) { return 0; } return 2; } + + if (abort == true) { if (calledHelp) { return 0; } return 2; } m->mothurOut("\nuchime by Robert C. Edgar\nhttp://drive5.com/uchime\nThis code is donated to the public domain.\n\n"); @@ -476,28 +630,47 @@ int ChimeraUchimeCommand::execute(){ int start = time(NULL); string nameFile = ""; if (outputDir == "") { outputDir = m->hasPath(fastaFileNames[s]); }//if user entered a file with a path then preserve it - string outputFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "uchime.chimera"; - string accnosFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "uchime.accnos"; - string alnsFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "uchime.alns"; + string outputFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + getOutputFileNameTag("chimera"); + string accnosFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + getOutputFileNameTag("accnos"); + string alnsFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + getOutputFileNameTag("alns"); string newFasta = m->getRootName(fastaFileNames[s]) + "temp"; //you provided a groupfile string groupFile = ""; - if (groupFileNames.size() != 0) { groupFile = groupFileNames[s]; } + bool hasGroup = false; + if (groupFileNames.size() != 0) { groupFile = groupFileNames[s]; hasGroup = true; } + else if (hasCount) { + CountTable ct; + if (ct.testGroups(nameFileNames[s])) { hasGroup = true; } + } - if ((templatefile == "self") && (groupFile == "")) { //you want to run uchime with a reference template + if ((templatefile == "self") && (!hasGroup)) { //you want to run uchime with a template=self and no groups if (processors != 1) { m->mothurOut("When using template=self, mothur can only use 1 processor, continuing."); m->mothurOutEndLine(); processors = 1; } if (nameFileNames.size() != 0) { //you provided a namefile and we don't need to create one nameFile = nameFileNames[s]; }else { nameFile = getNamesFile(fastaFileNames[s]); } - + map seqs; readFasta(fastaFileNames[s], seqs); if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } //read namefile vector nameMapCount; - int error = m->readNames(nameFile, nameMapCount, seqs); if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } + int error; + if (hasCount) { + CountTable ct; + ct.readTable(nameFile); + for(map::iterator it = seqs.begin(); it != seqs.end(); it++) { + int num = ct.getNumSeqs(it->first); + if (num == 0) { error = 1; } + else { + seqPriorityNode temp(num, it->second, it->first); + nameMapCount.push_back(temp); + } + } + }else { + error = m->readNames(nameFile, nameMapCount, seqs); if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } + } if (error == 1) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } if (seqs.size() != nameMapCount.size()) { m->mothurOut( "The number of sequences in your fastafile does not match the number of sequences in your namefile, aborting."); m->mothurOutEndLine(); for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } @@ -507,14 +680,23 @@ int ChimeraUchimeCommand::execute(){ if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } - if (groupFile != "") { + if (hasGroup) { if (nameFileNames.size() != 0) { //you provided a namefile and we don't need to create one nameFile = nameFileNames[s]; }else { nameFile = getNamesFile(fastaFileNames[s]); } //Parse sequences by group - SequenceParser parser(groupFile, fastaFileNames[s], nameFile); - vector groups = parser.getNamesOfGroups(); + vector groups; + map uniqueNames; + if (hasCount) { + cparser = new SequenceCountParser(nameFile, fastaFileNames[s]); + groups = cparser->getNamesOfGroups(); + uniqueNames = cparser->getAllSeqsMap(); + }else{ + sparser = new SequenceParser(groupFile, fastaFileNames[s], nameFile); + groups = sparser->getNamesOfGroups(); + uniqueNames = sparser->getAllSeqsMap(); + } if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } @@ -525,15 +707,14 @@ int ChimeraUchimeCommand::execute(){ if (chimealns) { m->openOutputFile(alnsFileName, out2); out2.close(); } int totalSeqs = 0; - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors == 1) { totalSeqs = driverGroups(parser, outputFileName, newFasta, accnosFileName, alnsFileName, 0, groups.size(), groups); } - else { totalSeqs = createProcessesGroups(parser, outputFileName, newFasta, accnosFileName, alnsFileName, groups); } - #else - totalSeqs = driverGroups(parser, outputFileName, newFasta, accnosFileName, alnsFileName, 0, groups.size(), groups); - #endif + if(processors == 1) { totalSeqs = driverGroups(outputFileName, newFasta, accnosFileName, alnsFileName, 0, groups.size(), groups); } + else { totalSeqs = createProcessesGroups(outputFileName, newFasta, accnosFileName, alnsFileName, groups, nameFile, groupFile, fastaFileNames[s]); } + if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } + if (hasCount) { delete cparser; } + else { delete sparser; } - int totalChimeras = deconvoluteResults(parser, outputFileName, accnosFileName, alnsFileName); + int totalChimeras = deconvoluteResults(uniqueNames, outputFileName, accnosFileName, alnsFileName); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(totalSeqs) + " sequences. " + toString(totalChimeras) + " chimeras were found."); m->mothurOutEndLine(); m->mothurOut("The number of sequences checked may be larger than the number of unique sequences because some sequences are found in several samples."); m->mothurOutEndLine(); @@ -545,12 +726,19 @@ int ChimeraUchimeCommand::execute(){ int numSeqs = 0; int numChimeras = 0; - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + if(processors == 1){ numSeqs = driver(outputFileName, fastaFileNames[s], accnosFileName, alnsFileName, numChimeras); } else{ numSeqs = createProcesses(outputFileName, fastaFileNames[s], accnosFileName, alnsFileName, numChimeras); } - #else - numSeqs = driver(outputFileName, fastaFileNames[s], accnosFileName, alnsFileName, numChimeras); - #endif + + //add headings + ofstream out; + m->openOutputFile(outputFileName+".temp", out); + out << "Score\tQuery\tParentA\tParentB\tIdQM\tIdQA\tIdQB\tIdAB\tIdQT\tLY\tLN\tLA\tRY\tRN\tRA\tDiv\tYN\n"; + out.close(); + + m->appendFiles(outputFileName, outputFileName+".temp"); + m->mothurRemove(outputFileName); rename((outputFileName+".temp").c_str(), outputFileName.c_str()); + if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } //remove file made for uchime @@ -585,26 +773,66 @@ int ChimeraUchimeCommand::execute(){ } } //********************************************************************************************************************** -int ChimeraUchimeCommand::deconvoluteResults(SequenceParser& parser, string outputFileName, string accnosFileName, string alnsFileName){ +int ChimeraUchimeCommand::deconvoluteResults(map& uniqueNames, string outputFileName, string accnosFileName, string alnsFileName){ try { - map uniqueNames = parser.getAllSeqsMap(); map::iterator itUnique; int total = 0; + //edit accnos file + ifstream in2; + m->openInputFile(accnosFileName, in2); + + ofstream out2; + m->openOutputFile(accnosFileName+".temp", out2); + + string name; + set namesInFile; //this is so if a sequence is found to be chimera in several samples we dont write it to the results file more than once + set::iterator itNames; + set chimerasInFile; + set::iterator itChimeras; + + + while (!in2.eof()) { + if (m->control_pressed) { in2.close(); out2.close(); m->mothurRemove(outputFileName); m->mothurRemove((accnosFileName+".temp")); return 0; } + + in2 >> name; m->gobble(in2); + + //find unique name + itUnique = uniqueNames.find(name); + + if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing accnos results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; } + else { + itChimeras = chimerasInFile.find((itUnique->second)); + + if (itChimeras == chimerasInFile.end()) { + out2 << itUnique->second << endl; + chimerasInFile.insert((itUnique->second)); + total++; + } + } + } + in2.close(); + out2.close(); + + m->mothurRemove(accnosFileName); + rename((accnosFileName+".temp").c_str(), accnosFileName.c_str()); + + + //edit chimera file ifstream in; m->openInputFile(outputFileName, in); ofstream out; m->openOutputFile(outputFileName+".temp", out); out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); + out << "Score\tQuery\tParentA\tParentB\tIdQM\tIdQA\tIdQB\tIdAB\tIdQT\tLY\tLN\tLA\tRY\tRN\tRA\tDiv\tYN\n"; float temp1; - string name, rest, parent1, parent2; - set namesInFile; //this is so if a sequence is found to be chimera in several samples we dont write it to the results file more than once - set::iterator itNames; - + string parent1, parent2, temp2, temp3, temp4, temp5, temp6, temp7, temp8, temp9, temp10, temp11, temp12, temp13, flag; + name = ""; + namesInFile.clear(); //assumptions - in file each read will always look like - if uchime source is updated, revisit this code. - /* + /* 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0.000000 F11Fcsw_33372/ab=18/ * * * * * * * * * * * * * * N 0.018300 F11Fcsw_14980/ab=16/ F11Fcsw_1915/ab=35/ F11Fcsw_6032/ab=42/ 79.9 78.7 78.2 78.7 79.2 3 0 5 11 10 20 1.46 N */ @@ -613,11 +841,13 @@ int ChimeraUchimeCommand::deconvoluteResults(SequenceParser& parser, string outp if (m->control_pressed) { in.close(); out.close(); m->mothurRemove((outputFileName+".temp")); return 0; } + bool print = false; in >> temp1; m->gobble(in); in >> name; m->gobble(in); in >> parent1; m->gobble(in); in >> parent2; m->gobble(in); - rest = m->getline(in); m->gobble(in); + in >> temp2 >> temp3 >> temp4 >> temp5 >> temp6 >> temp7 >> temp8 >> temp9 >> temp10 >> temp11 >> temp12 >> temp13 >> flag; + m->gobble(in); //parse name - name will look like U68590/ab=1/ string restOfName = ""; @@ -632,46 +862,54 @@ int ChimeraUchimeCommand::deconvoluteResults(SequenceParser& parser, string outp if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; } else { - itNames = namesInFile.find((itUnique->second)); + name = itUnique->second; + //is this name already in the file + itNames = namesInFile.find((name)); - if (itNames == namesInFile.end()) { - out << temp1 << '\t' << itUnique->second << restOfName << '\t'; - namesInFile.insert((itUnique->second)); - - //parse parent1 names - if (parent1 != "*") { - restOfName = ""; - pos = parent1.find_first_of('/'); - if (pos != string::npos) { - restOfName = parent1.substr(pos); - parent1 = parent1.substr(0, pos); - } + if (itNames == namesInFile.end()) { //no not in file + if (flag == "N") { //are you really a no?? + //is this sequence really not chimeric?? + itChimeras = chimerasInFile.find(name); - itUnique = uniqueNames.find(parent1); - if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parentA "+ parent1 + "."); m->mothurOutEndLine(); m->control_pressed = true; } - else { - out << itUnique->second << restOfName << '\t'; - } - }else { out << parent1 << '\t'; } + //then you really are a no so print, otherwise skip + if (itChimeras == chimerasInFile.end()) { print = true; } + }else{ print = true; } + } + } + + if (print) { + out << temp1 << '\t' << name << restOfName << '\t'; + namesInFile.insert(name); + + //parse parent1 names + if (parent1 != "*") { + restOfName = ""; + pos = parent1.find_first_of('/'); + if (pos != string::npos) { + restOfName = parent1.substr(pos); + parent1 = parent1.substr(0, pos); + } - //parse parent2 names - if (parent2 != "*") { - restOfName = ""; - pos = parent2.find_first_of('/'); - if (pos != string::npos) { - restOfName = parent2.substr(pos); - parent2 = parent2.substr(0, pos); - } - - itUnique = uniqueNames.find(parent2); - if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parentB "+ parent2 + "."); m->mothurOutEndLine(); m->control_pressed = true; } - else { - out << itUnique->second << restOfName << '\t'; - } - }else { out << parent2 << '\t'; } + itUnique = uniqueNames.find(parent1); + if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parentA "+ parent1 + "."); m->mothurOutEndLine(); m->control_pressed = true; } + else { out << itUnique->second << restOfName << '\t'; } + }else { out << parent1 << '\t'; } + + //parse parent2 names + if (parent2 != "*") { + restOfName = ""; + pos = parent2.find_first_of('/'); + if (pos != string::npos) { + restOfName = parent2.substr(pos); + parent2 = parent2.substr(0, pos); + } - out << rest << endl; - } + itUnique = uniqueNames.find(parent2); + if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parentB "+ parent2 + "."); m->mothurOutEndLine(); m->control_pressed = true; } + else { out << itUnique->second << restOfName << '\t'; } + }else { out << parent2 << '\t'; } + + out << temp2 << '\t' << temp3 << '\t' << temp4 << '\t' << temp5 << '\t' << temp6 << '\t' << temp7 << '\t' << temp8 << '\t' << temp9 << '\t' << temp10 << '\t' << temp11 << '\t' << temp12 << temp13 << '\t' << flag << endl; } } in.close(); @@ -680,41 +918,7 @@ int ChimeraUchimeCommand::deconvoluteResults(SequenceParser& parser, string outp m->mothurRemove(outputFileName); rename((outputFileName+".temp").c_str(), outputFileName.c_str()); - //edit accnos file - ifstream in2; - m->openInputFile(accnosFileName, in2); - - ofstream out2; - m->openOutputFile(accnosFileName+".temp", out2); - - name = ""; - namesInFile.clear(); - - while (!in2.eof()) { - if (m->control_pressed) { in2.close(); out2.close(); m->mothurRemove(outputFileName); m->mothurRemove((accnosFileName+".temp")); return 0; } - - in2 >> name; m->gobble(in2); - - //find unique name - itUnique = uniqueNames.find(name); - - if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing accnos results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; } - else { - itNames = namesInFile.find((itUnique->second)); - - if (itNames == namesInFile.end()) { - out2 << itUnique->second << endl; - namesInFile.insert((itUnique->second)); - total++; - } - } - } - in2.close(); - out2.close(); - - m->mothurRemove(accnosFileName); - rename((accnosFileName+".temp").c_str(), accnosFileName.c_str()); - + //edit anls file //assumptions - in file each read will always look like - if uchime source is updated, revisit this code. /* @@ -888,14 +1092,15 @@ string ChimeraUchimeCommand::getNamesFile(string& inputFile){ string inputString = "fasta=" + inputFile; m->mothurOut("/******************************************/"); m->mothurOutEndLine(); m->mothurOut("Running command: unique.seqs(" + inputString + ")"); m->mothurOutEndLine(); - + m->mothurCalling = true; + Command* uniqueCommand = new DeconvoluteCommand(inputString); uniqueCommand->execute(); map > filenames = uniqueCommand->getOutputFiles(); delete uniqueCommand; - + m->mothurCalling = false; m->mothurOut("/******************************************/"); m->mothurOutEndLine(); nameFile = filenames["name"][0]; @@ -909,7 +1114,7 @@ string ChimeraUchimeCommand::getNamesFile(string& inputFile){ } } //********************************************************************************************************************** -int ChimeraUchimeCommand::driverGroups(SequenceParser& parser, string outputFName, string filename, string accnos, string alns, int start, int end, vector groups){ +int ChimeraUchimeCommand::driverGroups(string outputFName, string filename, string accnos, string alns, int start, int end, vector groups){ try { int totalSeqs = 0; @@ -917,8 +1122,10 @@ int ChimeraUchimeCommand::driverGroups(SequenceParser& parser, string outputFNam for (int i = start; i < end; i++) { int start = time(NULL); if (m->control_pressed) { return 0; } - - int error = parser.getSeqs(groups[i], filename, true); if ((error == 1) || m->control_pressed) { return 0; } + + int error; + if (hasCount) { error = cparser->getSeqs(groups[i], filename, true); if ((error == 1) || m->control_pressed) { return 0; } } + else { error = sparser->getSeqs(groups[i], filename, true); if ((error == 1) || m->control_pressed) { return 0; } } int numSeqs = driver((outputFName + groups[i]), filename, (accnos+ groups[i]), (alns+ groups[i]), numChimeras); totalSeqs += numSeqs; @@ -926,7 +1133,8 @@ int ChimeraUchimeCommand::driverGroups(SequenceParser& parser, string outputFNam if (m->control_pressed) { return 0; } //remove file made for uchime - m->mothurRemove(filename); + if (!m->debug) { m->mothurRemove(filename); } + else { m->mothurOut("[DEBUG]: saving file: " + filename + ".\n"); } //append files m->appendFiles((outputFName+groups[i]), outputFName); m->mothurRemove((outputFName+groups[i])); @@ -935,7 +1143,6 @@ int ChimeraUchimeCommand::driverGroups(SequenceParser& parser, string outputFNam m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences from group " + groups[i] + "."); m->mothurOutEndLine(); } - return totalSeqs; } @@ -949,23 +1156,35 @@ int ChimeraUchimeCommand::driverGroups(SequenceParser& parser, string outputFNam int ChimeraUchimeCommand::driver(string outputFName, string filename, string accnos, string alns, int& numChimeras){ try { + outputFName = m->getFullPathName(outputFName); + filename = m->getFullPathName(filename); + alns = m->getFullPathName(alns); + + //to allow for spaces in the path + outputFName = "\"" + outputFName + "\""; + filename = "\"" + filename + "\""; + alns = "\"" + alns + "\""; + vector cPara; - char* tempUchime = new char[8]; - strcpy(tempUchime, "./uchime "); + string uchimeCommand = uchimeLocation; +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + uchimeCommand += " "; +#else + uchimeCommand = "\"" + uchimeCommand + "\""; +#endif + char* tempUchime; + tempUchime= new char[uchimeCommand.length()+1]; + *tempUchime = '\0'; + strncat(tempUchime, uchimeCommand.c_str(), uchimeCommand.length()); cPara.push_back(tempUchime); - char* tempIn = new char[8]; - *tempIn = '\0'; strncat(tempIn, "--input", 7); - //strcpy(tempIn, "--input"); - cPara.push_back(tempIn); - char* temp = new char[filename.length()+1]; - *temp = '\0'; strncat(temp, filename.c_str(), filename.length()); - //strcpy(temp, filename.c_str()); - cPara.push_back(temp); - - //are you using a reference file + //are you using a reference file if (templatefile != "self") { + string outputFileName = filename.substr(1, filename.length()-2) + ".uchime_formatted"; + prepFile(filename.substr(1, filename.length()-2), outputFileName); + filename = outputFileName; + filename = "\"" + filename + "\""; //add reference file char* tempRef = new char[5]; //strcpy(tempRef, "--db"); @@ -977,6 +1196,15 @@ int ChimeraUchimeCommand::driver(string outputFName, string filename, string acc cPara.push_back(tempR); } + char* tempIn = new char[8]; + *tempIn = '\0'; strncat(tempIn, "--input", 7); + //strcpy(tempIn, "--input"); + cPara.push_back(tempIn); + char* temp = new char[filename.length()+1]; + *temp = '\0'; strncat(temp, filename.c_str(), filename.length()); + //strcpy(temp, filename.c_str()); + cPara.push_back(temp); + char* tempO = new char[12]; *tempO = '\0'; strncat(tempO, "--uchimeout", 11); //strcpy(tempO, "--uchimeout"); @@ -1174,15 +1402,28 @@ int ChimeraUchimeCommand::driver(string outputFName, string filename, string acc char** uchimeParameters; uchimeParameters = new char*[cPara.size()]; - for (int i = 0; i < cPara.size(); i++) { uchimeParameters[i] = cPara[i]; } - int numArgs = cPara.size(); - - uchime_main(numArgs, uchimeParameters); + string commandString = ""; + for (int i = 0; i < cPara.size(); i++) { uchimeParameters[i] = cPara[i]; commandString += toString(cPara[i]) + " "; } + //int numArgs = cPara.size(); + + //uchime_main(numArgs, uchimeParameters); + //cout << "commandString = " << commandString << endl; +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) +#else + commandString = "\"" + commandString + "\""; +#endif + if (m->debug) { m->mothurOut("[DEBUG]: uchime command = " + commandString + ".\n"); } + system(commandString.c_str()); //free memory - for(int i = 0; i < cPara.size(); i++) { delete[] cPara[i]; } + for(int i = 0; i < cPara.size(); i++) { delete cPara[i]; } delete[] uchimeParameters; + //remove "" from filenames + outputFName = outputFName.substr(1, outputFName.length()-2); + filename = filename.substr(1, filename.length()-2); + alns = alns.substr(1, alns.length()-2); + if (m->control_pressed) { return 0; } //create accnos file from uchime results @@ -1217,6 +1458,8 @@ int ChimeraUchimeCommand::driver(string outputFName, string filename, string acc in.close(); out.close(); + //if (templatefile != "self") { m->mothurRemove(filename); } + return num; } catch(exception& e) { @@ -1225,6 +1468,34 @@ int ChimeraUchimeCommand::driver(string outputFName, string filename, string acc } } /**************************************************************************************************/ +//uchime can't handle some of the things allowed in mothurs fasta files. This functions "cleans up" the file. +int ChimeraUchimeCommand::prepFile(string filename, string output) { + try { + + ifstream in; + m->openInputFile(filename, in); + + ofstream out; + m->openOutputFile(output, out); + + while (!in.eof()) { + if (m->control_pressed) { break; } + + Sequence seq(in); m->gobble(in); + + if (seq.getName() != "") { seq.printSequence(out); } + } + in.close(); + out.close(); + + return 0; + } + catch(exception& e) { + m->errorOut(e, "ChimeraUchimeCommand", "prepFile"); + exit(1); + } +} +/**************************************************************************************************/ int ChimeraUchimeCommand::createProcesses(string outputFileName, string filename, string accnos, string alns, int& numChimeras) { try { @@ -1232,9 +1503,10 @@ int ChimeraUchimeCommand::createProcesses(string outputFileName, string filename processIDS.clear(); int process = 1; int num = 0; -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - //break up file into multiple files vector files; + +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + //break up file into multiple files m->divideFile(filename, processors, files); if (m->control_pressed) { return 0; } @@ -1287,10 +1559,92 @@ int ChimeraUchimeCommand::createProcesses(string outputFileName, string filename } in.close(); m->mothurRemove(tempFile); } +#else + ////////////////////////////////////////////////////////////////////////////////////////////////////// + //Windows version shared memory, so be careful when passing variables through the preClusterData struct. + //Above fork() will clone, so memory is separate, but that's not the case with windows, + ////////////////////////////////////////////////////////////////////////////////////////////////////// + + //divide file + int count = 0; + int spot = 0; + map filehandles; + map::iterator it3; + + ofstream* temp; + for (int i = 0; i < processors; i++) { + temp = new ofstream; + filehandles[i] = temp; + m->openOutputFile(filename+toString(i)+".temp", *(temp)); + files.push_back(filename+toString(i)+".temp"); + } + + ifstream in; + m->openInputFile(filename, in); + + while(!in.eof()) { + + if (m->control_pressed) { in.close(); for (it3 = filehandles.begin(); it3 != filehandles.end(); it3++) { (*(it3->second)).close(); delete it3->second; } return 0; } + + Sequence tempSeq(in); m->gobble(in); + + if (tempSeq.getName() != "") { + tempSeq.printSequence(*(filehandles[spot])); + spot++; count++; + if (spot == processors) { spot = 0; } + } + } + in.close(); + + //delete memory + for (it3 = filehandles.begin(); it3 != filehandles.end(); it3++) { + (*(it3->second)).close(); + delete it3->second; + } + + //sanity check for number of processors + if (count < processors) { processors = count; } + + vector pDataArray; + DWORD dwThreadIdArray[processors-1]; + HANDLE hThreadArray[processors-1]; + vector dummy; //used so that we can use the same struct for MyUchimeSeqsThreadFunction and MyUchimeThreadFunction + + //Create processor worker threads. + for( int i=1; isetBooleans(useAbskew, chimealns, useMinH, useMindiv, useXn, useDn, useXa, useChunks, useMinchunk, useIdsmoothwindow, useMinsmoothid, useMaxp, skipgaps, skipgaps2, useMinlen, useMaxlen, ucl, useQueryfract, hasCount); + tempUchime->setVariables(abskew, minh, mindiv, xn, dn, xa, chunks, minchunk, idsmoothwindow, minsmoothid, maxp, minlen, maxlen, queryfract); + + pDataArray.push_back(tempUchime); + processIDS.push_back(i); + + //MySeqSumThreadFunction is in header. It must be global or static to work with the threads. + //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier + hThreadArray[i-1] = CreateThread(NULL, 0, MyUchimeSeqsThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]); + } + + + //using the main process as a worker saves time and memory + num = driver(outputFileName, files[0], accnos, alns, numChimeras); + + //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; + numChimeras += pDataArray[i]->numChimeras; + CloseHandle(hThreadArray[i]); + delete pDataArray[i]; + } +#endif //append output files - for(int i=0;iappendFiles((outputFileName + toString(processIDS[i]) + ".temp"), outputFileName); m->mothurRemove((outputFileName + toString(processIDS[i]) + ".temp")); @@ -1305,7 +1659,6 @@ int ChimeraUchimeCommand::createProcesses(string outputFileName, string filename //get rid of the file pieces. for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } -#endif return num; } catch(exception& e) { @@ -1315,7 +1668,7 @@ int ChimeraUchimeCommand::createProcesses(string outputFileName, string filename } /**************************************************************************************************/ -int ChimeraUchimeCommand::createProcessesGroups(SequenceParser& parser, string outputFName, string filename, string accnos, string alns, vector groups) { +int ChimeraUchimeCommand::createProcessesGroups(string outputFName, string filename, string accnos, string alns, vector groups, string nameFile, string groupFile, string fastaFile) { try { processIDS.clear(); @@ -1335,7 +1688,7 @@ int ChimeraUchimeCommand::createProcessesGroups(SequenceParser& parser, string o lines.push_back(linePair(startIndex, endIndex)); } -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) //loop through and create all the processes you want while (process != processors) { @@ -1345,7 +1698,7 @@ int ChimeraUchimeCommand::createProcessesGroups(SequenceParser& parser, string o 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 = driverGroups(parser, outputFName + toString(getpid()) + ".temp", filename + toString(getpid()) + ".temp", accnos + toString(getpid()) + ".temp", alns + toString(getpid()) + ".temp", lines[process].start, lines[process].end, groups); + num = driverGroups(outputFName + toString(getpid()) + ".temp", filename + toString(getpid()) + ".temp", accnos + toString(getpid()) + ".temp", alns + toString(getpid()) + ".temp", lines[process].start, lines[process].end, groups); //pass numSeqs to parent ofstream out; @@ -1363,14 +1716,13 @@ int ChimeraUchimeCommand::createProcessesGroups(SequenceParser& parser, string o } //do my part - num = driverGroups(parser, outputFName, filename, accnos, alns, lines[0].start, lines[0].end, groups); + num = driverGroups(outputFName, filename, accnos, alns, lines[0].start, lines[0].end, groups); //force parent to wait until all the processes are done for (int i=0;i> tempNum; num += tempNum; } in.close(); m->mothurRemove(tempFile); } + +#else + ////////////////////////////////////////////////////////////////////////////////////////////////////// + //Windows version shared memory, so be careful when passing variables through the uchimeData 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=1; isetBooleans(useAbskew, chimealns, useMinH, useMindiv, useXn, useDn, useXa, useChunks, useMinchunk, useIdsmoothwindow, useMinsmoothid, useMaxp, skipgaps, skipgaps2, useMinlen, useMaxlen, ucl, useQueryfract, hasCount); + tempUchime->setVariables(abskew, minh, mindiv, xn, dn, xa, chunks, minchunk, idsmoothwindow, minsmoothid, maxp, minlen, maxlen, queryfract); + + pDataArray.push_back(tempUchime); + processIDS.push_back(i); + + //MyUchimeThreadFunction 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-1] = CreateThread(NULL, 0, MyUchimeThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]); + } + + //using the main process as a worker saves time and memory + num = driverGroups(outputFName, filename, accnos, alns, lines[0].start, lines[0].end, groups); + //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 output files - for(int i=0;iappendFiles((outputFName + toString(processIDS[i]) + ".temp"), outputFName); m->mothurRemove((outputFName + toString(processIDS[i]) + ".temp"));