X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=trimseqscommand.cpp;h=89ea3f9c8b44c782e2a9ebd8588c548c61f9e81a;hb=7f0cae4f4853cc3f12bc751ee06ea31c7c97496e;hp=d1b9a5c43549991bd3876164a318b8d9accb54ce;hpb=e840ba0e90a2b2be39e0910e4be23dde6a14cbac;p=mothur.git diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp index d1b9a5c..89ea3f9 100644 --- a/trimseqscommand.cpp +++ b/trimseqscommand.cpp @@ -9,67 +9,107 @@ #include "trimseqscommand.h" #include "needlemanoverlap.hpp" +#include "trimoligos.h" //********************************************************************************************************************** - -vector TrimSeqsCommand::getValidParameters(){ +vector TrimSeqsCommand::setParameters(){ try { - string Array[] = {"fasta", "flip", "oligos", "maxambig", "maxhomop","minlength", "maxlength", "qfile", - "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage", - "keepfirst", "removelast", - "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta); + CommandParameter poligos("oligos", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(poligos); + CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pqfile); + CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname); + CommandParameter pflip("flip", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pflip); + CommandParameter pmaxambig("maxambig", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pmaxambig); + CommandParameter pmaxhomop("maxhomop", "Number", "", "0", "", "", "",false,false); parameters.push_back(pmaxhomop); + CommandParameter pminlength("minlength", "Number", "", "0", "", "", "",false,false); parameters.push_back(pminlength); + CommandParameter pmaxlength("maxlength", "Number", "", "0", "", "", "",false,false); parameters.push_back(pmaxlength); + CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(ppdiffs); + CommandParameter pbdiffs("bdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(pbdiffs); + CommandParameter ptdiffs("tdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(ptdiffs); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); + CommandParameter pallfiles("allfiles", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pallfiles); + CommandParameter pqtrim("qtrim", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pqtrim); + CommandParameter pqthreshold("qthreshold", "Number", "", "0", "", "", "",false,false); parameters.push_back(pqthreshold); + CommandParameter pqaverage("qaverage", "Number", "", "0", "", "", "",false,false); parameters.push_back(pqaverage); + CommandParameter prollaverage("rollaverage", "Number", "", "0", "", "", "",false,false); parameters.push_back(prollaverage); + CommandParameter pqwindowaverage("qwindowaverage", "Number", "", "0", "", "", "",false,false); parameters.push_back(pqwindowaverage); + CommandParameter pqstepsize("qstepsize", "Number", "", "1", "", "", "",false,false); parameters.push_back(pqstepsize); + CommandParameter pqwindowsize("qwindowsize", "Number", "", "50", "", "", "",false,false); parameters.push_back(pqwindowsize); + CommandParameter pkeepfirst("keepfirst", "Number", "", "0", "", "", "",false,false); parameters.push_back(pkeepfirst); + CommandParameter premovelast("removelast", "Number", "", "0", "", "", "",false,false); parameters.push_back(premovelast); + CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); + CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); + + vector myArray; + for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } return myArray; } catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "getValidParameters"); + m->errorOut(e, "TrimSeqsCommand", "setParameters"); exit(1); } } - //********************************************************************************************************************** - -TrimSeqsCommand::TrimSeqsCommand(){ +string TrimSeqsCommand::getHelpString(){ try { - abort = true; calledHelp = true; - vector tempOutNames; - outputTypes["fasta"] = tempOutNames; - outputTypes["qfile"] = tempOutNames; - outputTypes["group"] = tempOutNames; + string helpString = ""; + helpString += "The trim.seqs command reads a fastaFile and creates 2 new fasta files, .trim.fasta and scrap.fasta, as well as group files if you provide and oligos file.\n"; + helpString += "The .trim.fasta contains sequences that meet your requirements, and the .scrap.fasta contains those which don't.\n"; + helpString += "The trim.seqs command parameters are fasta, name, flip, oligos, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, diffs, qtrim, keepfirst, removelast and allfiles.\n"; + helpString += "The fasta parameter is required.\n"; + helpString += "The flip parameter will output the reverse compliment of your trimmed sequence. The default is false.\n"; + helpString += "The oligos parameter allows you to provide an oligos file.\n"; + helpString += "The name parameter allows you to provide a names file with your fasta file.\n"; + helpString += "The maxambig parameter allows you to set the maximum number of ambigious bases allowed. The default is -1.\n"; + helpString += "The maxhomop parameter allows you to set a maximum homopolymer length. \n"; + helpString += "The minlength parameter allows you to set and minimum sequence length. \n"; + helpString += "The maxlength parameter allows you to set and maximum sequence length. \n"; + helpString += "The tdiffs parameter is used to specify the total number of differences allowed in the sequence. The default is pdiffs + bdiffs.\n"; + helpString += "The bdiffs parameter is used to specify the number of differences allowed in the barcode. The default is 0.\n"; + helpString += "The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n"; + helpString += "The qfile parameter allows you to provide a quality file.\n"; + helpString += "The qthreshold parameter allows you to set a minimum quality score allowed. \n"; + helpString += "The qaverage parameter allows you to set a minimum average quality score allowed. \n"; + helpString += "The qwindowsize parameter allows you to set a number of bases in a window. Default=50.\n"; + helpString += "The qwindowaverage parameter allows you to set a minimum average quality score allowed over a window. \n"; + helpString += "The rollaverage parameter allows you to set a minimum rolling average quality score allowed over a window. \n"; + helpString += "The qstepsize parameter allows you to set a number of bases to move the window over. Default=1.\n"; + helpString += "The allfiles parameter will create separate group and fasta file for each grouping. The default is F.\n"; + helpString += "The qtrim parameter will trim sequence from the point that they fall below the qthreshold and put it in the .trim file if set to true. The default is T.\n"; + helpString += "The keepfirst parameter trims the sequence to the first keepfirst number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements. \n"; + helpString += "The removelast removes the last removelast number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements.\n"; + helpString += "The trim.seqs command should be in the following format: \n"; + helpString += "trim.seqs(fasta=yourFastaFile, flip=yourFlip, oligos=yourOligos, maxambig=yourMaxambig, \n"; + helpString += "maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n"; + helpString += "Example trim.seqs(fasta=abrecovery.fasta, flip=..., oligos=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n"; + helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n"; + helpString += "For more details please check out the wiki http://www.mothur.org/wiki/Trim.seqs .\n"; + return helpString; } catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "TrimSeqsCommand"); + m->errorOut(e, "TrimSeqsCommand", "getHelpString"); exit(1); } } -//********************************************************************************************************************** - -vector TrimSeqsCommand::getRequiredParameters(){ - try { - string Array[] = {"fasta"}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); - return myArray; - } - catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "getRequiredParameters"); - exit(1); - } -} //********************************************************************************************************************** -vector TrimSeqsCommand::getRequiredFiles(){ +TrimSeqsCommand::TrimSeqsCommand(){ try { - vector myArray; - return myArray; + abort = true; calledHelp = true; + setParameters(); + vector tempOutNames; + outputTypes["fasta"] = tempOutNames; + outputTypes["qfile"] = tempOutNames; + outputTypes["group"] = tempOutNames; + outputTypes["name"] = tempOutNames; } catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "getRequiredFiles"); + m->errorOut(e, "TrimSeqsCommand", "TrimSeqsCommand"); exit(1); } } - //*************************************************************************************************************** TrimSeqsCommand::TrimSeqsCommand(string option) { @@ -80,15 +120,10 @@ TrimSeqsCommand::TrimSeqsCommand(string option) { //allow user to run help if(option == "help") { help(); abort = true; calledHelp = true; } + else if(option == "citation") { citation(); abort = true; calledHelp = true;} else { - //valid paramters for this command - string AlignArray[] = { "fasta", "flip", "oligos", "maxambig", "maxhomop", "minlength", "maxlength", "qfile", - "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage", - "keepfirst", "removelast", - "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"}; - - vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); + vector myArray = setParameters(); OptionParser parser(option); map parameters = parser.getParameters(); @@ -106,6 +141,7 @@ TrimSeqsCommand::TrimSeqsCommand(string option) { outputTypes["fasta"] = tempOutNames; outputTypes["qfile"] = tempOutNames; outputTypes["group"] = tempOutNames; + outputTypes["name"] = 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); @@ -136,13 +172,25 @@ TrimSeqsCommand::TrimSeqsCommand(string option) { if (path == "") { parameters["qfile"] = inputDir + it->second; } } + it = parameters.find("name"); + //user has given a template file + if(it != parameters.end()){ + path = m->hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["name"] = inputDir + it->second; } + } + } //check for required parameters fastaFile = validParameter.validFile(parameters, "fasta", true); - if (fastaFile == "not found") { m->mothurOut("fasta is a required parameter for the trim.seqs command."); m->mothurOutEndLine(); abort = true; } - else if (fastaFile == "not open") { abort = true; } + if (fastaFile == "not found") { + fastaFile = m->getFastaFile(); + if (fastaFile != "") { m->mothurOut("Using " + fastaFile + " as input file for the fasta parameter."); m->mothurOutEndLine(); } + else { m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; } + }else if (fastaFile == "not open") { abort = true; } + else { m->setFastaFile(fastaFile); } //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"){ @@ -155,13 +203,13 @@ TrimSeqsCommand::TrimSeqsCommand(string option) { // ...at some point should added some additional type checking... string temp; temp = validParameter.validFile(parameters, "flip", false); - if (temp == "not found"){ flip = 0; } - else if(m->isTrue(temp)) { flip = 1; } + if (temp == "not found") { flip = 0; } + else { flip = m->isTrue(temp); } temp = validParameter.validFile(parameters, "oligos", true); if (temp == "not found"){ oligoFile = ""; } else if(temp == "not open"){ abort = true; } - else { oligoFile = temp; } + else { oligoFile = temp; m->setOligosFile(oligoFile); } temp = validParameter.validFile(parameters, "maxambig", false); if (temp == "not found") { temp = "-1"; } @@ -190,7 +238,12 @@ TrimSeqsCommand::TrimSeqsCommand(string option) { temp = validParameter.validFile(parameters, "qfile", true); if (temp == "not found") { qFileName = ""; } else if(temp == "not open") { abort = true; } - else { qFileName = temp; } + else { qFileName = temp; m->setQualFile(qFileName); } + + temp = validParameter.validFile(parameters, "name", true); + if (temp == "not found") { nameFile = ""; } + else if(temp == "not open") { nameFile = ""; abort = true; } + else { nameFile = temp; m->setNameFile(nameFile); } temp = validParameter.validFile(parameters, "qthreshold", false); if (temp == "not found") { temp = "0"; } convert(temp, qThreshold); @@ -222,7 +275,8 @@ TrimSeqsCommand::TrimSeqsCommand(string option) { temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found") { temp = "F"; } allFiles = m->isTrue(temp); - temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found") { temp = "1"; } + temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } + m->setProcessors(temp); convert(temp, processors); @@ -246,54 +300,6 @@ TrimSeqsCommand::TrimSeqsCommand(string option) { exit(1); } } - -//********************************************************************************************************************** - -void TrimSeqsCommand::help(){ - try { - m->mothurOut("The trim.seqs command reads a fastaFile and creates 2 new fasta files, .trim.fasta and scrap.fasta, as well as group files if you provide and oligos file.\n"); - m->mothurOut("The .trim.fasta contains sequences that meet your requirements, and the .scrap.fasta contains those which don't.\n"); - m->mothurOut("The trim.seqs command parameters are fasta, flip, oligos, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, diffs, qtrim, keepfirst, removelast and allfiles.\n"); - m->mothurOut("The fasta parameter is required.\n"); - m->mothurOut("The flip parameter will output the reverse compliment of your trimmed sequence. The default is false.\n"); - m->mothurOut("The oligos parameter allows you to provide an oligos file.\n"); - m->mothurOut("The maxambig parameter allows you to set the maximum number of ambigious bases allowed. The default is -1.\n"); - m->mothurOut("The maxhomop parameter allows you to set a maximum homopolymer length. \n"); - m->mothurOut("The minlength parameter allows you to set and minimum sequence length. \n"); - m->mothurOut("The maxlength parameter allows you to set and maximum sequence length. \n"); - m->mothurOut("The tdiffs parameter is used to specify the total number of differences allowed in the sequence. The default is pdiffs + bdiffs.\n"); - m->mothurOut("The bdiffs parameter is used to specify the number of differences allowed in the barcode. The default is 0.\n"); - m->mothurOut("The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n"); - m->mothurOut("The qfile parameter allows you to provide a quality file.\n"); - m->mothurOut("The qthreshold parameter allows you to set a minimum quality score allowed. \n"); - m->mothurOut("The qaverage parameter allows you to set a minimum average quality score allowed. \n"); - m->mothurOut("The qwindowsize parameter allows you to set a number of bases in a window. Default=50.\n"); - m->mothurOut("The qwindowaverage parameter allows you to set a minimum average quality score allowed over a window. \n"); - m->mothurOut("The rollaverage parameter allows you to set a minimum rolling average quality score allowed over a window. \n"); - m->mothurOut("The qstepsize parameter allows you to set a number of bases to move the window over. Default=1.\n"); - m->mothurOut("The allfiles parameter will create separate group and fasta file for each grouping. The default is F.\n"); - m->mothurOut("The qtrim parameter will trim sequence from the point that they fall below the qthreshold and put it in the .trim file if set to true. The default is T.\n"); - m->mothurOut("The keepfirst parameter trims the sequence to the first keepfirst number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements. \n"); - m->mothurOut("The removelast removes the last removelast number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements.\n"); - m->mothurOut("The trim.seqs command should be in the following format: \n"); - m->mothurOut("trim.seqs(fasta=yourFastaFile, flip=yourFlip, oligos=yourOligos, maxambig=yourMaxambig, \n"); - m->mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n"); - m->mothurOut("Example trim.seqs(fasta=abrecovery.fasta, flip=..., oligos=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n"); - m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n"); - m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/Trim.seqs .\n\n"); - - } - catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "help"); - exit(1); - } -} - - -//*************************************************************************************************************** - -TrimSeqsCommand::~TrimSeqsCommand(){ /* do nothing */ } - //*************************************************************************************************************** int TrimSeqsCommand::execute(){ @@ -303,8 +309,10 @@ int TrimSeqsCommand::execute(){ numFPrimers = 0; //this needs to be initialized numRPrimers = 0; + createGroup = false; vector > fastaFileNames; vector > qualFileNames; + vector > nameFileNames; string trimSeqFile = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "trim.fasta"; outputNames.push_back(trimSeqFile); outputTypes["fasta"].push_back(trimSeqFile); @@ -314,6 +322,7 @@ int TrimSeqsCommand::execute(){ string trimQualFile = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "trim.qual"; string scrapQualFile = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "scrap.qual"; + if (qFileName != "") { outputNames.push_back(trimQualFile); outputNames.push_back(scrapQualFile); @@ -321,16 +330,30 @@ int TrimSeqsCommand::execute(){ outputTypes["qfile"].push_back(scrapQualFile); } + string trimNameFile = outputDir + m->getRootName(m->getSimpleName(nameFile)) + "trim.names"; + string scrapNameFile = outputDir + m->getRootName(m->getSimpleName(nameFile)) + "scrap.names"; + + if (nameFile != "") { + m->readNames(nameFile, nameMap); + outputNames.push_back(trimNameFile); + outputNames.push_back(scrapNameFile); + outputTypes["name"].push_back(trimNameFile); + outputTypes["name"].push_back(scrapNameFile); + } + + if (m->control_pressed) { return 0; } + string outputGroupFileName; - if(oligoFile != ""){ - outputGroupFileName = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "groups"; - outputNames.push_back(outputGroupFileName); outputTypes["group"].push_back(outputGroupFileName); - getOligos(fastaFileNames, qualFileNames); + createGroup = getOligos(fastaFileNames, qualFileNames, nameFileNames); + if (createGroup) { + outputGroupFileName = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "groups"; + outputNames.push_back(outputGroupFileName); outputTypes["group"].push_back(outputGroupFileName); + } } - - vector fastaFilePos; - vector qFilePos; + + vector fastaFilePos; + vector qFilePos; setLines(fastaFile, qFileName, fastaFilePos, qFilePos); @@ -342,37 +365,106 @@ int TrimSeqsCommand::execute(){ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) if(processors == 1){ - driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, outputGroupFileName, fastaFileNames, qualFileNames, lines[0], qLines[0]); + driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, trimNameFile, scrapNameFile, outputGroupFileName, fastaFileNames, qualFileNames, nameFileNames, lines[0], qLines[0]); }else{ - createProcessesCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, outputGroupFileName, fastaFileNames, qualFileNames); + createProcessesCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, trimNameFile, scrapNameFile, outputGroupFileName, fastaFileNames, qualFileNames, nameFileNames); } #else - driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, outputGroupFileName, fastaFileNames, qualFileNames, lines[0], qLines[0]); + driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, trimNameFile, scrapNameFile, outputGroupFileName, fastaFileNames, qualFileNames, nameFileNames, lines[0], qLines[0]); #endif if (m->control_pressed) { return 0; } - - + if(allFiles){ + map uniqueFastaNames;// so we don't add the same groupfile multiple times + map::iterator it; + set namesToRemove; for(int i=0;iisBlank(fastaFileNames[i][j])){ - remove(fastaFileNames[i][j].c_str()); - - if(qFileName != ""){ - remove(fastaFileNames[i][j].c_str()); + if (fastaFileNames[i][j] != "") { + if (namesToRemove.count(fastaFileNames[i][j]) == 0) { + if(m->isBlank(fastaFileNames[i][j])){ + m->mothurRemove(fastaFileNames[i][j]); + namesToRemove.insert(fastaFileNames[i][j]); + + if(qFileName != ""){ + m->mothurRemove(qualFileNames[i][j]); + namesToRemove.insert(qualFileNames[i][j]); + } + + if(nameFile != ""){ + m->mothurRemove(nameFileNames[i][j]); + namesToRemove.insert(nameFileNames[i][j]); + } + }else{ + it = uniqueFastaNames.find(fastaFileNames[i][j]); + if (it == uniqueFastaNames.end()) { + uniqueFastaNames[fastaFileNames[i][j]] = barcodeNameVector[i]; + } + } } - } } } + + //remove names for outputFileNames, just cleans up the output + vector outputNames2; + for(int i = 0; i < outputNames.size(); i++) { if (namesToRemove.count(outputNames[i]) == 0) { outputNames2.push_back(outputNames[i]); } } + outputNames = outputNames2; + + for (it = uniqueFastaNames.begin(); it != uniqueFastaNames.end(); it++) { + ifstream in; + m->openInputFile(it->first, in); + + ofstream out; + string thisGroupName = outputDir + m->getRootName(m->getSimpleName(it->first)) + "groups"; + outputNames.push_back(thisGroupName); outputTypes["group"].push_back(thisGroupName); + m->openOutputFile(thisGroupName, out); + + while (!in.eof()){ + if (m->control_pressed) { break; } + + Sequence currSeq(in); m->gobble(in); + out << currSeq.getName() << '\t' << it->second << endl; + } + in.close(); + out.close(); + } + } + + if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } + + //output group counts + m->mothurOutEndLine(); + int total = 0; + if (groupCounts.size() != 0) { m->mothurOut("Group count: \n"); } + for (map::iterator it = groupCounts.begin(); it != groupCounts.end(); it++) { + total += it->second; m->mothurOut(it->first + "\t" + toString(it->second)); m->mothurOutEndLine(); } + if (total != 0) { m->mothurOut("Total of all groups is " + toString(total)); m->mothurOutEndLine(); } + if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } + + //set fasta file as new current fastafile + string current = ""; + itTypes = outputTypes.find("fasta"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); } + } + + itTypes = outputTypes.find("name"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); } + } + itTypes = outputTypes.find("qfile"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setQualFile(current); } + } - if (m->control_pressed) { - for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } - return 0; + itTypes = outputTypes.find("group"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); } } m->mothurOutEndLine(); @@ -391,7 +483,7 @@ int TrimSeqsCommand::execute(){ /**************************************************************************************/ -int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string trimFileName, string scrapFileName, string trimQFileName, string scrapQFileName, string groupFileName, vector > fastaFileNames, vector > qualFileNames, linePair* line, linePair* qline) { +int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string trimFileName, string scrapFileName, string trimQFileName, string scrapQFileName, string trimNFileName, string scrapNFileName, string groupFileName, vector > fastaFileNames, vector > qualFileNames, vector > nameFileNames, linePair* line, linePair* qline) { try { @@ -408,16 +500,29 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string m->openOutputFile(scrapQFileName, scrapQualFile); } - ofstream outGroupsFile; - if (oligoFile != ""){ m->openOutputFile(groupFileName, outGroupsFile); } + ofstream trimNameFile; + ofstream scrapNameFile; + if(nameFile != ""){ + m->openOutputFile(trimNFileName, trimNameFile); + m->openOutputFile(scrapNFileName, scrapNameFile); + } + + ofstream outGroupsFile; + if (createGroup){ m->openOutputFile(groupFileName, outGroupsFile); } if(allFiles){ for (int i = 0; i < fastaFileNames.size(); i++) { //clears old file for (int j = 0; j < fastaFileNames[i].size(); j++) { //clears old file - ofstream temp; - m->openOutputFile(fastaFileNames[i][j], temp); temp.close(); - if(qFileName != ""){ - m->openOutputFile(qualFileNames[i][j], temp); temp.close(); + if (fastaFileNames[i][j] != "") { + ofstream temp; + m->openOutputFile(fastaFileNames[i][j], temp); temp.close(); + if(qFileName != ""){ + m->openOutputFile(qualFileNames[i][j], temp); temp.close(); + } + + if(nameFile != ""){ + m->openOutputFile(nameFileNames[i][j], temp); temp.close(); + } } } } @@ -435,17 +540,18 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string int count = 0; bool moreSeqs = 1; + TrimOligos trimOligos(pdiffs, bdiffs, primers, barcodes, revPrimer); while (moreSeqs) { if (m->control_pressed) { inFASTA.close(); trimFASTAFile.close(); scrapFASTAFile.close(); - if (oligoFile != "") { outGroupsFile.close(); } + if (createGroup) { outGroupsFile.close(); } if(qFileName != ""){ qFile.close(); } - for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } + for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } @@ -455,12 +561,12 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string int currentSeqsDiffs = 0; Sequence currSeq(inFASTA); m->gobble(inFASTA); - + //cout << currSeq.getName() << '\t' << currSeq.getUnaligned().length() << endl; QualityScores currQual; if(qFileName != ""){ currQual = QualityScores(qFile); m->gobble(qFile); } - + string origSeq = currSeq.getUnaligned(); if (origSeq != "") { @@ -468,13 +574,13 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string int primerIndex = 0; if(barcodes.size() != 0){ - success = stripBarcode(currSeq, currQual, barcodeIndex); + success = trimOligos.stripBarcode(currSeq, currQual, barcodeIndex); if(success > bdiffs) { trashCode += 'b'; } else{ currentSeqsDiffs += success; } } if(numFPrimers != 0){ - success = stripForward(currSeq, currQual, primerIndex); + success = trimOligos.stripForward(currSeq, currQual, primerIndex); if(success > pdiffs) { trashCode += 'f'; } else{ currentSeqsDiffs += success; } } @@ -482,7 +588,7 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string if (currentSeqsDiffs > tdiffs) { trashCode += 't'; } if(numRPrimers != 0){ - success = stripReverse(currSeq, currQual); + success = trimOligos.stripReverse(currSeq, currQual); if(!success) { trashCode += 'r'; } } @@ -539,10 +645,44 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string currQual.printQScores(trimQualFile); } - if(barcodes.size() != 0){ - outGroupsFile << currSeq.getName() << '\t' << barcodeNameVector[barcodeIndex] << endl; + if(nameFile != ""){ + map::iterator itName = nameMap.find(currSeq.getName()); + if (itName != nameMap.end()) { trimNameFile << itName->first << '\t' << itName->second << endl; } + else { m->mothurOut("[ERROR]: " + currSeq.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); } } + if (createGroup) { + if(barcodes.size() != 0){ + string thisGroup = barcodeNameVector[barcodeIndex]; + if (primers.size() != 0) { + if (primerNameVector[primerIndex] != "") { + if(thisGroup != "") { + thisGroup += "." + primerNameVector[primerIndex]; + }else { + thisGroup = primerNameVector[primerIndex]; + } + } + } + + outGroupsFile << currSeq.getName() << '\t' << thisGroup << endl; + + if (nameFile != "") { + map::iterator itName = nameMap.find(currSeq.getName()); + if (itName != nameMap.end()) { + vector thisSeqsNames; + m->splitAtChar(itName->second, thisSeqsNames, ','); + for (int k = 1; k < thisSeqsNames.size(); k++) { //start at 1 to skip self + outGroupsFile << thisSeqsNames[k] << '\t' << thisGroup << endl; + } + }else { m->mothurOut("[ERROR]: " + currSeq.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); } + } + + map::iterator it = groupCounts.find(thisGroup); + if (it == groupCounts.end()) { groupCounts[thisGroup] = 1; } + else { groupCounts[it->first]++; } + + } + } if(allFiles){ ofstream output; @@ -555,9 +695,23 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string currQual.printQScores(output); output.close(); } + + if(nameFile != ""){ + map::iterator itName = nameMap.find(currSeq.getName()); + if (itName != nameMap.end()) { + m->openOutputFileAppend(nameFileNames[barcodeIndex][primerIndex], output); + output << itName->first << '\t' << itName->second << endl; + output.close(); + }else { m->mothurOut("[ERROR]: " + currSeq.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); } + } } } else{ + if(nameFile != ""){ //needs to be before the currSeq name is changed + map::iterator itName = nameMap.find(currSeq.getName()); + if (itName != nameMap.end()) { scrapNameFile << itName->first << '\t' << itName->second << endl; } + else { m->mothurOut("[ERROR]: " + currSeq.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); } + } currSeq.setName(currSeq.getName() + '|' + trashCode); currSeq.setUnaligned(origSeq); currSeq.setAligned(origSeq); @@ -570,25 +724,27 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string } #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - unsigned long int pos = inFASTA.tellg(); + unsigned long long pos = inFASTA.tellg(); if ((pos == -1) || (pos >= line->end)) { break; } + #else if (inFASTA.eof()) { break; } #endif - + //report progress if((count) % 1000 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } } //report progress if((count) % 1000 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } - + inFASTA.close(); trimFASTAFile.close(); scrapFASTAFile.close(); - if (oligoFile != "") { outGroupsFile.close(); } + if (createGroup) { outGroupsFile.close(); } if(qFileName != "") { qFile.close(); scrapQualFile.close(); trimQualFile.close(); } + if(nameFile != "") { scrapNameFile.close(); trimNameFile.close(); } return count; } @@ -600,7 +756,7 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string /**************************************************************************************************/ -int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName, string trimFASTAFileName, string scrapFASTAFileName, string trimQualFileName, string scrapQualFileName, string groupFile, vector > fastaFileNames, vector > qualFileNames) { +int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName, string trimFASTAFileName, string scrapFASTAFileName, string trimQualFileName, string scrapQualFileName, string trimNameFileName, string scrapNameFileName, string groupFile, vector > fastaFileNames, vector > qualFileNames, vector > nameFileNames) { try { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) int process = 1; @@ -618,18 +774,25 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName vector > tempFASTAFileNames = fastaFileNames; vector > tempPrimerQualFileNames = qualFileNames; + vector > tempNameFileNames = nameFileNames; if(allFiles){ ofstream temp; for(int i=0;iopenOutputFile(tempFASTAFileNames[i][j], temp); temp.close(); - - if(qFileName != ""){ - tempPrimerQualFileNames[i][j] += toString(getpid()) + ".temp"; - m->openOutputFile(tempPrimerQualFileNames[i][j], temp); temp.close(); + if (tempFASTAFileNames[i][j] != "") { + tempFASTAFileNames[i][j] += toString(getpid()) + ".temp"; + m->openOutputFile(tempFASTAFileNames[i][j], temp); temp.close(); + + if(qFileName != ""){ + tempPrimerQualFileNames[i][j] += toString(getpid()) + ".temp"; + m->openOutputFile(tempPrimerQualFileNames[i][j], temp); temp.close(); + } + if(nameFile != ""){ + tempNameFileNames[i][j] += toString(getpid()) + ".temp"; + m->openOutputFile(tempNameFileNames[i][j], temp); temp.close(); + } } } } @@ -641,12 +804,28 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName (scrapFASTAFileName + toString(getpid()) + ".temp"), (trimQualFileName + toString(getpid()) + ".temp"), (scrapQualFileName + toString(getpid()) + ".temp"), + (trimNameFileName + toString(getpid()) + ".temp"), + (scrapNameFileName + toString(getpid()) + ".temp"), (groupFile + toString(getpid()) + ".temp"), tempFASTAFileNames, tempPrimerQualFileNames, + tempNameFileNames, lines[process], qLines[process]); + //pass groupCounts to parent + if(createGroup){ + ofstream out; + string tempFile = filename + toString(getpid()) + ".num.temp"; + m->openOutputFile(tempFile, out); + + out << groupCounts.size() << endl; + + for (map::iterator it = groupCounts.begin(); it != groupCounts.end(); it++) { + out << it->first << '\t' << it->second << endl; + } + out.close(); + } exit(0); }else { m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); @@ -659,13 +838,16 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName ofstream temp; m->openOutputFile(trimFASTAFileName, temp); temp.close(); m->openOutputFile(scrapFASTAFileName, temp); temp.close(); - m->openOutputFile(trimQualFileName, temp); temp.close(); - m->openOutputFile(scrapQualFileName, temp); temp.close(); + if(qFileName != ""){ + m->openOutputFile(trimQualFileName, temp); temp.close(); + m->openOutputFile(scrapQualFileName, temp); temp.close(); + } + if (nameFile != "") { + m->openOutputFile(trimNameFileName, temp); temp.close(); + m->openOutputFile(scrapNameFileName, temp); temp.close(); + } - - - driverCreateTrim(filename, qFileName, trimFASTAFileName, scrapFASTAFileName, trimQualFileName, scrapQualFileName, groupFile, fastaFileNames, qualFileNames, lines[0], qLines[0]); - + driverCreateTrim(filename, qFileName, trimFASTAFileName, scrapFASTAFileName, trimQualFileName, scrapQualFileName, trimNameFileName, scrapNameFileName, groupFile, fastaFileNames, qualFileNames, nameFileNames, lines[0], qLines[0]); //force parent to wait until all the processes are done for (int i=0;imothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine(); m->appendFiles((trimFASTAFileName + toString(processIDS[i]) + ".temp"), trimFASTAFileName); - remove((trimFASTAFileName + toString(processIDS[i]) + ".temp").c_str()); + m->mothurRemove((trimFASTAFileName + toString(processIDS[i]) + ".temp")); m->appendFiles((scrapFASTAFileName + toString(processIDS[i]) + ".temp"), scrapFASTAFileName); - remove((scrapFASTAFileName + toString(processIDS[i]) + ".temp").c_str()); + m->mothurRemove((scrapFASTAFileName + toString(processIDS[i]) + ".temp")); if(qFileName != ""){ m->appendFiles((trimQualFileName + toString(processIDS[i]) + ".temp"), trimQualFileName); - remove((trimQualFileName + toString(processIDS[i]) + ".temp").c_str()); + m->mothurRemove((trimQualFileName + toString(processIDS[i]) + ".temp")); m->appendFiles((scrapQualFileName + toString(processIDS[i]) + ".temp"), scrapQualFileName); - remove((scrapQualFileName + toString(processIDS[i]) + ".temp").c_str()); + m->mothurRemove((scrapQualFileName + toString(processIDS[i]) + ".temp")); } - m->appendFiles((groupFile + toString(processIDS[i]) + ".temp"), groupFile); - remove((groupFile + toString(processIDS[i]) + ".temp").c_str()); + if(nameFile != ""){ + m->appendFiles((trimNameFileName + toString(processIDS[i]) + ".temp"), trimNameFileName); + m->mothurRemove((trimNameFileName + toString(processIDS[i]) + ".temp")); + m->appendFiles((scrapNameFileName + toString(processIDS[i]) + ".temp"), scrapNameFileName); + m->mothurRemove((scrapNameFileName + toString(processIDS[i]) + ".temp")); + } + + if(createGroup){ + m->appendFiles((groupFile + toString(processIDS[i]) + ".temp"), groupFile); + m->mothurRemove((groupFile + toString(processIDS[i]) + ".temp")); + } if(allFiles){ for(int j=0;jappendFiles((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp"), fastaFileNames[j][k]); - remove((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp").c_str()); - - if(qFileName != ""){ - m->appendFiles((qualFileNames[j][k] + toString(processIDS[i]) + ".temp"), qualFileNames[j][k]); - remove((qualFileNames[j][k] + toString(processIDS[i]) + ".temp").c_str()); + if (fastaFileNames[j][k] != "") { + m->appendFiles((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp"), fastaFileNames[j][k]); + m->mothurRemove((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp")); + + if(qFileName != ""){ + m->appendFiles((qualFileNames[j][k] + toString(processIDS[i]) + ".temp"), qualFileNames[j][k]); + m->mothurRemove((qualFileNames[j][k] + toString(processIDS[i]) + ".temp")); + } + + if(nameFile != ""){ + m->appendFiles((nameFileNames[j][k] + toString(processIDS[i]) + ".temp"), nameFileNames[j][k]); + m->mothurRemove((nameFileNames[j][k] + toString(processIDS[i]) + ".temp")); + } } } } } + if(createGroup){ + ifstream in; + string tempFile = filename + toString(processIDS[i]) + ".num.temp"; + m->openInputFile(tempFile, in); + int tempNum; + string group; + + in >> tempNum; m->gobble(in); + + if (tempNum != 0) { + while (!in.eof()) { + in >> group >> tempNum; m->gobble(in); + + map::iterator it = groupCounts.find(group); + if (it == groupCounts.end()) { groupCounts[group] = tempNum; } + else { groupCounts[it->first] += tempNum; } + } + } + in.close(); m->mothurRemove(tempFile); + } + } return exitCommand; @@ -721,9 +940,9 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName /**************************************************************************************************/ -int TrimSeqsCommand::setLines(string filename, string qfilename, vector& fastaFilePos, vector& qfileFilePos) { +int TrimSeqsCommand::setLines(string filename, string qfilename, vector& fastaFilePos, vector& qfileFilePos) { try { - + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) //set file positions for fasta file fastaFilePos = m->divideFile(filename, processors); @@ -760,7 +979,7 @@ int TrimSeqsCommand::setLines(string filename, string qfilename, vector::iterator it = firstSeqNames.find(sname); if(it != firstSeqNames.end()) { //this is the start of a new chunk - unsigned long int pos = inQual.tellg(); + unsigned long long pos = inQual.tellg(); qfileFilePos.push_back(pos - input.length() - 1); firstSeqNames.erase(it); } @@ -782,7 +1001,7 @@ int TrimSeqsCommand::setLines(string filename, string qfilename, vectorerrorOut(e, "TrimSeqsCommand", "setLines"); @@ -805,7 +1058,7 @@ int TrimSeqsCommand::setLines(string filename, string qfilename, vector >& fastaFileNames, vector >& qualFileNames){ +bool TrimSeqsCommand::getOligos(vector >& fastaFileNames, vector >& qualFileNames, vector >& nameFileNames){ try { ifstream inOligos; m->openInputFile(oligoFile, inOligos); @@ -819,12 +1072,14 @@ void TrimSeqsCommand::getOligos(vector >& fastaFileNames, vector< while(!inOligos.eof()){ - inOligos >> type; m->gobble(inOligos); + inOligos >> type; if(type[0] == '#'){ - while (!inOligos.eof()) { char c = inOligos.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there + while (!inOligos.eof()) { char c = inOligos.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there + m->gobble(inOligos); } else{ + m->gobble(inOligos); //make type case insensitive for(int i=0;i >& fastaFileNames, vector< for(int i=0;i uniqueNames; //used to cleanup outputFileNames for(map::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){ for(map::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){ @@ -903,6 +1160,7 @@ void TrimSeqsCommand::getOligos(vector >& fastaFileNames, vector< string comboGroupName = ""; string fastaFileName = ""; string qualFileName = ""; + string nameFileName = ""; if(primerName == ""){ comboGroupName = barcodeNameVector[itBar->second]; @@ -915,306 +1173,75 @@ void TrimSeqsCommand::getOligos(vector >& fastaFileNames, vector< comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second]; } } - + + ofstream temp; fastaFileName = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + comboGroupName + ".fasta"; - outputNames.push_back(fastaFileName); - outputTypes["fasta"].push_back(fastaFileName); + if (uniqueNames.count(fastaFileName) == 0) { + outputNames.push_back(fastaFileName); + outputTypes["fasta"].push_back(fastaFileName); + uniqueNames.insert(fastaFileName); + } + fastaFileNames[itBar->second][itPrimer->second] = fastaFileName; m->openOutputFile(fastaFileName, temp); temp.close(); - + if(qFileName != ""){ qualFileName = outputDir + m->getRootName(m->getSimpleName(qFileName)) + comboGroupName + ".qual"; - outputNames.push_back(qualFileName); - outputTypes["qfile"].push_back(qualFileName); + if (uniqueNames.count(qualFileName) == 0) { + outputNames.push_back(qualFileName); + outputTypes["qfile"].push_back(qualFileName); + } + qualFileNames[itBar->second][itPrimer->second] = qualFileName; m->openOutputFile(qualFileName, temp); temp.close(); } + + if(nameFile != ""){ + nameFileName = outputDir + m->getRootName(m->getSimpleName(nameFile)) + comboGroupName + ".names"; + if (uniqueNames.count(nameFileName) == 0) { + outputNames.push_back(nameFileName); + outputTypes["name"].push_back(nameFileName); + } + + nameFileNames[itBar->second][itPrimer->second] = nameFileName; + m->openOutputFile(nameFileName, temp); temp.close(); + } + } } } numFPrimers = primers.size(); numRPrimers = revPrimer.size(); - - } - catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "getOligos"); - exit(1); - } -} - -//*************************************************************************************************************** - -int TrimSeqsCommand::stripBarcode(Sequence& seq, QualityScores& qual, int& group){ - try { - string rawSequence = seq.getUnaligned(); - int success = bdiffs + 1; //guilty until proven innocent - - //can you find the barcode - for(map::iterator it=barcodes.begin();it!=barcodes.end();it++){ - string oligo = it->first; - if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length - success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out - break; - } - - if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){ - group = it->second; - seq.setUnaligned(rawSequence.substr(oligo.length())); - - if(qual.getName() != ""){ - qual.trimQScores(oligo.length(), -1); - } - - success = 0; + bool allBlank = true; + for (int i = 0; i < barcodeNameVector.size(); i++) { + if (barcodeNameVector[i] != "") { + allBlank = false; break; } } - - //if you found the barcode or if you don't want to allow for diffs - if ((bdiffs == 0) || (success == 0)) { return success; } - - else { //try aligning and see if you can find it - - int maxLength = 0; - - Alignment* alignment; - if (barcodes.size() > 0) { - map::iterator it=barcodes.begin(); - - for(it;it!=barcodes.end();it++){ - if(it->first.length() > maxLength){ - maxLength = it->first.length(); - } - } - alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1)); - - }else{ alignment = NULL; } - - //can you find the barcode - int minDiff = 1e6; - int minCount = 1; - int minGroup = -1; - int minPos = 0; - - for(map::iterator it=barcodes.begin();it!=barcodes.end();it++){ - string oligo = it->first; -// int length = oligo.length(); - - if(rawSequence.length() < maxLength){ //let's just assume that the barcodes are the same length - success = bdiffs + 10; - break; - } - - //use needleman to align first barcode.length()+numdiffs of sequence to each barcode - alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs)); - oligo = alignment->getSeqAAln(); - string temp = alignment->getSeqBAln(); - - int alnLength = oligo.length(); - - for(int i=oligo.length()-1;i>=0;i--){ - if(oligo[i] != '-'){ alnLength = i+1; break; } - } - oligo = oligo.substr(0,alnLength); - temp = temp.substr(0,alnLength); - - int numDiff = countDiffs(oligo, temp); - - if(numDiff < minDiff){ - minDiff = numDiff; - minCount = 1; - minGroup = it->second; - minPos = 0; - for(int i=0;i bdiffs) { success = minDiff; } //no good matches - else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes - else{ //use the best match - group = minGroup; - seq.setUnaligned(rawSequence.substr(minPos)); - - if(qual.getName() != ""){ - qual.trimQScores(minPos, -1); - } - success = minDiff; - } - - if (alignment != NULL) { delete alignment; } - - } - - return success; - - } - catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "stripBarcode"); - exit(1); - } - -} - -//*************************************************************************************************************** - -int TrimSeqsCommand::stripForward(Sequence& seq, QualityScores& qual, int& group){ - try { - string rawSequence = seq.getUnaligned(); - int success = pdiffs + 1; //guilty until proven innocent - - //can you find the primer - for(map::iterator it=primers.begin();it!=primers.end();it++){ - string oligo = it->first; - if(rawSequence.length() < oligo.length()){ //let's just assume that the primers are the same length - success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out - break; - } - - if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){ - group = it->second; - seq.setUnaligned(rawSequence.substr(oligo.length())); - if(qual.getName() != ""){ - qual.trimQScores(oligo.length(), -1); - } - success = 0; + for (int i = 0; i < primerNameVector.size(); i++) { + if (primerNameVector[i] != "") { + allBlank = false; break; } } - - //if you found the barcode or if you don't want to allow for diffs - if ((pdiffs == 0) || (success == 0)) { return success; } - - else { //try aligning and see if you can find it - - int maxLength = 0; - - Alignment* alignment; - if (primers.size() > 0) { - map::iterator it=primers.begin(); - - for(it;it!=primers.end();it++){ - if(it->first.length() > maxLength){ - maxLength = it->first.length(); - } - } - alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1)); - - }else{ alignment = NULL; } - - //can you find the barcode - int minDiff = 1e6; - int minCount = 1; - int minGroup = -1; - int minPos = 0; - - for(map::iterator it=primers.begin();it!=primers.end();it++){ - string oligo = it->first; -// int length = oligo.length(); - - if(rawSequence.length() < maxLength){ - success = pdiffs + 100; - break; - } - - //use needleman to align first barcode.length()+numdiffs of sequence to each barcode - alignment->align(oligo, rawSequence.substr(0,oligo.length()+pdiffs)); - oligo = alignment->getSeqAAln(); - string temp = alignment->getSeqBAln(); - int alnLength = oligo.length(); - - for(int i=oligo.length()-1;i>=0;i--){ - if(oligo[i] != '-'){ alnLength = i+1; break; } - } - oligo = oligo.substr(0,alnLength); - temp = temp.substr(0,alnLength); - - int numDiff = countDiffs(oligo, temp); - - if(numDiff < minDiff){ - minDiff = numDiff; - minCount = 1; - minGroup = it->second; - minPos = 0; - for(int i=0;i pdiffs) { success = minDiff; } //no good matches - else if(minCount > 1) { success = pdiffs + 10; } //can't tell the difference between multiple primers - else{ //use the best match - group = minGroup; - seq.setUnaligned(rawSequence.substr(minPos)); - if(qual.getName() != ""){ - qual.trimQScores(minPos, -1); - } - success = minDiff; - } - - if (alignment != NULL) { delete alignment; } - + if (allBlank) { + m->mothurOut("[WARNING]: your oligos file does not contain any group names. mothur will not create a groupfile."); m->mothurOutEndLine(); + allFiles = false; + return false; } - return success; - - } - catch(exception& e) { - m->errorOut(e, "TrimSeqsCommand", "stripForward"); - exit(1); - } -} - -//*************************************************************************************************************** - -bool TrimSeqsCommand::stripReverse(Sequence& seq, QualityScores& qual){ - try { - string rawSequence = seq.getUnaligned(); - bool success = 0; //guilty until proven innocent - - for(int i=0;ierrorOut(e, "TrimSeqsCommand", "stripReverse"); + m->errorOut(e, "TrimSeqsCommand", "getOligos"); exit(1); } } - //*************************************************************************************************************** bool TrimSeqsCommand::keepFirstTrim(Sequence& sequence, QualityScores& qscores){ @@ -1320,80 +1347,4 @@ bool TrimSeqsCommand::cullAmbigs(Sequence& seq){ } } - -//*************************************************************************************************************** - -bool TrimSeqsCommand::compareDNASeq(string oligo, string seq){ - try { - bool success = 1; - int length = oligo.length(); - - for(int i=0;ierrorOut(e, "TrimSeqsCommand", "compareDNASeq"); - exit(1); - } - -} - -//*************************************************************************************************************** - -int TrimSeqsCommand::countDiffs(string oligo, string seq){ - try { - - int length = oligo.length(); - int countDiffs = 0; - - for(int i=0;ierrorOut(e, "TrimSeqsCommand", "countDiffs"); - exit(1); - } - -} - //***************************************************************************************************************