X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=removeseqscommand.cpp;h=25d45a7c578c47167fac24164cac3322f3462772;hb=250e3b11b1c9c1e1ad458ab6c7e71ac2e67e11d9;hp=1fb94462ac48fa13ac871c6b9d3cf4f3ae2ceff6;hpb=ea4f373c28543cd1002b0dd7dc6e55c526647d59;p=mothur.git diff --git a/removeseqscommand.cpp b/removeseqscommand.cpp index 1fb9446..25d45a7 100644 --- a/removeseqscommand.cpp +++ b/removeseqscommand.cpp @@ -10,22 +10,24 @@ #include "removeseqscommand.h" #include "sequence.hpp" #include "listvector.hpp" +#include "counttable.h" //********************************************************************************************************************** vector RemoveSeqsCommand::setParameters(){ try { - CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "FNGLT", "none",false,false); parameters.push_back(pfasta); - CommandParameter pname("name", "InputTypes", "", "", "NameCount", "FNGLT", "none",false,false); parameters.push_back(pname); - CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "FNGLT", "none",false,false); parameters.push_back(pcount); - CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "FNGLT", "none",false,false); parameters.push_back(pgroup); - CommandParameter plist("list", "InputTypes", "", "", "none", "FNGLT", "none",false,false); parameters.push_back(plist); - CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "none", "FNGLT", "none",false,false); parameters.push_back(ptaxonomy); - CommandParameter palignreport("alignreport", "InputTypes", "", "", "none", "FNGLT", "none",false,false); parameters.push_back(palignreport); - CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "FNGLT", "none",false,false); parameters.push_back(pqfile); - CommandParameter paccnos("accnos", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(paccnos); - CommandParameter pdups("dups", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pdups); - CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); - CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); + CommandParameter pfastq("fastq", "InputTypes", "", "", "none", "FNGLT", "none","fastq",false,false,true); parameters.push_back(pfastq); + CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "FNGLT", "none","fasta",false,false,true); parameters.push_back(pfasta); + CommandParameter pname("name", "InputTypes", "", "", "NameCount", "FNGLT", "none","name",false,false,true); parameters.push_back(pname); + CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "FNGLT", "none","count",false,false,true); parameters.push_back(pcount); + CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "FNGLT", "none","group",false,false,true); parameters.push_back(pgroup); + CommandParameter plist("list", "InputTypes", "", "", "none", "FNGLT", "none","list",false,false,true); parameters.push_back(plist); + CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "none", "FNGLT", "none","taxonomy",false,false,true); parameters.push_back(ptaxonomy); + CommandParameter palignreport("alignreport", "InputTypes", "", "", "none", "FNGLT", "none","alignreport",false,false); parameters.push_back(palignreport); + CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "FNGLT", "none","qfile",false,false); parameters.push_back(pqfile); + CommandParameter paccnos("accnos", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(paccnos); + CommandParameter pdups("dups", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pdups); + 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); } @@ -40,9 +42,9 @@ vector RemoveSeqsCommand::setParameters(){ string RemoveSeqsCommand::getHelpString(){ try { string helpString = ""; - helpString += "The remove.seqs command reads an .accnos file and at least one of the following file types: fasta, name, group, count, list, taxonomy, quality or alignreport file.\n"; + helpString += "The remove.seqs command reads an .accnos file and at least one of the following file types: fasta, name, group, count, list, taxonomy, quality, fastq or alignreport file.\n"; helpString += "It outputs a file containing the sequences NOT in the .accnos file.\n"; - helpString += "The remove.seqs command parameters are accnos, fasta, name, group, count, list, taxonomy, qfile, alignreport and dups. You must provide accnos and at least one of the file parameters.\n"; + helpString += "The remove.seqs command parameters are accnos, fasta, name, group, count, list, taxonomy, qfile, alignreport, fastq and dups. You must provide accnos and at least one of the file parameters.\n"; helpString += "The dups parameter allows you to remove the entire line from a name file if you remove any name from the line. default=true. \n"; helpString += "The remove.seqs command should be in the following format: remove.seqs(accnos=yourAccnos, fasta=yourFasta).\n"; helpString += "Example remove.seqs(accnos=amazon.accnos, fasta=amazon.fasta).\n"; @@ -55,33 +57,28 @@ string RemoveSeqsCommand::getHelpString(){ } } //********************************************************************************************************************** -string RemoveSeqsCommand::getOutputFileNameTag(string type, string inputName=""){ - try { - string outputFileName = ""; - map >::iterator it; +string RemoveSeqsCommand::getOutputPattern(string type) { + try { + string pattern = ""; - //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 == "fasta") { outputFileName = "pick" + m->getExtension(inputName); } - else if (type == "taxonomy") { outputFileName = "pick" + m->getExtension(inputName); } - else if (type == "name") { outputFileName = "pick" + m->getExtension(inputName); } - else if (type == "group") { outputFileName = "pick" + m->getExtension(inputName); } - else if (type == "list") { outputFileName = "pick" + m->getExtension(inputName); } - else if (type == "qfile") { outputFileName = "pick" + m->getExtension(inputName); } - else if (type == "alignreport") { outputFileName = "pick.align.report"; } - else if (type == "count") { outputFileName = "pick.count.table"; } - 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, "RemoveSeqsCommand", "getOutputFileNameTag"); - exit(1); - } + if (type == "fasta") { pattern = "[filename],pick,[extension]"; } + else if (type == "fastq") { pattern = "[filename],pick,[extension]"; } + else if (type == "taxonomy") { pattern = "[filename],pick,[extension]"; } + else if (type == "name") { pattern = "[filename],pick,[extension]"; } + else if (type == "group") { pattern = "[filename],pick,[extension]"; } + else if (type == "count") { pattern = "[filename],pick,[extension]"; } + else if (type == "list") { pattern = "[filename],[distance],pick,[extension]"; } + else if (type == "qfile") { pattern = "[filename],pick,[extension]"; } + else if (type == "alignreport") { pattern = "[filename],pick.align.report"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } + + return pattern; + } + catch(exception& e) { + m->errorOut(e, "GetSeqsCommand", "getOutputPattern"); + exit(1); + } } - //********************************************************************************************************************** RemoveSeqsCommand::RemoveSeqsCommand(){ try { @@ -89,6 +86,7 @@ RemoveSeqsCommand::RemoveSeqsCommand(){ setParameters(); vector tempOutNames; outputTypes["fasta"] = tempOutNames; + outputTypes["fastq"] = tempOutNames; outputTypes["taxonomy"] = tempOutNames; outputTypes["name"] = tempOutNames; outputTypes["group"] = tempOutNames; @@ -128,6 +126,7 @@ RemoveSeqsCommand::RemoveSeqsCommand(string option) { //initialize outputTypes vector tempOutNames; outputTypes["fasta"] = tempOutNames; + outputTypes["fastq"] = tempOutNames; outputTypes["taxonomy"] = tempOutNames; outputTypes["name"] = tempOutNames; outputTypes["group"] = tempOutNames; @@ -215,6 +214,14 @@ RemoveSeqsCommand::RemoveSeqsCommand(string option) { //if the user has not given a path then, add inputdir. else leave path alone. if (path == "") { parameters["count"] = inputDir + it->second; } } + + it = parameters.find("fastq"); + //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["fastq"] = inputDir + it->second; } + } } @@ -263,6 +270,10 @@ RemoveSeqsCommand::RemoveSeqsCommand(string option) { if (qualfile == "not open") { abort = true; } else if (qualfile == "not found") { qualfile = ""; } else { m->setQualFile(qualfile); } + + fastqfile = validParameter.validFile(parameters, "fastq", true); + if (fastqfile == "not open") { abort = true; } + else if (fastqfile == "not found") { fastqfile = ""; } string usedDups = "true"; string temp = validParameter.validFile(parameters, "dups", false); @@ -285,7 +296,7 @@ RemoveSeqsCommand::RemoveSeqsCommand(string option) { m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true; } - if ((countfile == "") && (fastafile == "") && (namefile == "") && (groupfile == "") && (alignfile == "") && (listfile == "") && (taxfile == "") && (qualfile == "")) { m->mothurOut("You must provide at least one of the following: fasta, name, group, taxonomy, quality, alignreport or list."); m->mothurOutEndLine(); abort = true; } + if ((fastqfile == "") && (countfile == "") && (fastafile == "") && (namefile == "") && (groupfile == "") && (alignfile == "") && (listfile == "") && (taxfile == "") && (qualfile == "")) { m->mothurOut("You must provide at least one of the following: fasta, name, group, taxonomy, quality, alignreport, fastq or list."); m->mothurOutEndLine(); abort = true; } if (countfile == "") { if ((fastafile != "") && (namefile == "")) { @@ -322,6 +333,7 @@ int RemoveSeqsCommand::execute(){ //read through the correct file and output lines you want to keep if (namefile != "") { readName(); } if (fastafile != "") { readFasta(); } + if (fastqfile != "") { readFastq(); } if (groupfile != "") { readGroup(); } if (alignfile != "") { readAlign(); } if (listfile != "") { readList(); } @@ -389,7 +401,10 @@ int RemoveSeqsCommand::readFasta(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(fastafile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)) + getOutputFileNameTag("fasta", fastafile); + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)); + variables["[extension]"] = m->getExtension(fastafile); + string outputFileName = getOutputFileName("fasta", variables); ofstream out; m->openOutputFile(outputFileName, out); @@ -405,6 +420,12 @@ int RemoveSeqsCommand::readFasta(){ if (m->control_pressed) { in.close(); out.close(); m->mothurRemove(outputFileName); return 0; } Sequence currSeq(in); + + if (!dups) {//adjust name if needed + map::iterator it = uniqueMap.find(currSeq.getName()); + if (it != uniqueMap.end()) { currSeq.setName(it->second); } + } + name = currSeq.getName(); if (name != "") { @@ -412,7 +433,7 @@ int RemoveSeqsCommand::readFasta(){ if (names.count(name) == 0) { wroteSomething = true; - currSeq.printSequence(out); + currSeq.printSequence(out); }else { removedCount++; } } m->gobble(in); @@ -434,11 +455,79 @@ int RemoveSeqsCommand::readFasta(){ } } //********************************************************************************************************************** +int RemoveSeqsCommand::readFastq(){ + try { + bool wroteSomething = false; + int removedCount = 0; + + ifstream in; + m->openInputFile(fastqfile, in); + + string thisOutputDir = outputDir; + if (outputDir == "") { thisOutputDir += m->hasPath(fastqfile); } + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(fastqfile)); + variables["[extension]"] = m->getExtension(fastqfile); + string outputFileName = getOutputFileName("fastq", variables); + ofstream out; + m->openOutputFile(outputFileName, out); + + + while(!in.eof()){ + + if (m->control_pressed) { in.close(); out.close(); m->mothurRemove(outputFileName); return 0; } + + //read sequence name + string input = m->getline(in); m->gobble(in); + + string outputString = input + "\n"; + + if (input[0] == '@') { + //get rest of lines + outputString += m->getline(in) + "\n"; m->gobble(in); + outputString += m->getline(in) + "\n"; m->gobble(in); + outputString += m->getline(in) + "\n"; m->gobble(in); + + vector splits = m->splitWhiteSpace(input); + string name = splits[0]; + name = name.substr(1); + m->checkName(name); + + if (names.count(name) == 0) { + wroteSomething = true; + out << outputString; + }else { removedCount++; } + } + + m->gobble(in); + } + in.close(); + out.close(); + + + if (wroteSomething == false) { m->mothurOut("Your file contains only sequences from the .accnos file."); m->mothurOutEndLine(); } + outputTypes["fasta"].push_back(outputFileName); outputNames.push_back(outputFileName); + + m->mothurOut("Removed " + toString(removedCount) + " sequences from your fastq file."); m->mothurOutEndLine(); + + + return 0; + + } + catch(exception& e) { + m->errorOut(e, "RemoveSeqsCommand", "readFastq"); + exit(1); + } +} +//********************************************************************************************************************** int RemoveSeqsCommand::readQual(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(qualfile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(qualfile)) + getOutputFileNameTag("qfile", qualfile); + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(qualfile)); + variables["[extension]"] = m->getExtension(qualfile); + string outputFileName = getOutputFileName("qfile", variables); ofstream out; m->openOutputFile(outputFileName, out); @@ -462,7 +551,7 @@ int RemoveSeqsCommand::readQual(){ saveName = name.substr(1); while (!in.eof()) { char c = in.get(); - if (c == 10 || c == 13){ break; } + if (c == 10 || c == 13 || c == -1){ break; } else { name += c; } } m->gobble(in); @@ -476,6 +565,11 @@ int RemoveSeqsCommand::readQual(){ m->gobble(in); + if (!dups) {//adjust name if needed + map::iterator it = uniqueMap.find(saveName); + if (it != uniqueMap.end()) { name = ">" + it->second; saveName = it->second; } + } + if (names.count(saveName) == 0) { wroteSomething = true; @@ -507,7 +601,10 @@ int RemoveSeqsCommand::readCount(){ string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(countfile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(countfile)) + getOutputFileNameTag("count", countfile); + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(countfile)); + variables["[extension]"] = m->getExtension(countfile); + string outputFileName = getOutputFileName("count", variables); ofstream out; m->openOutputFile(outputFileName, out); @@ -538,6 +635,14 @@ int RemoveSeqsCommand::readCount(){ } in.close(); out.close(); + + //check for groups that have been eliminated + CountTable ct; + if (ct.testGroups(outputFileName)) { + ct.readTable(outputFileName, true, false); + ct.printTable(outputFileName); + } + if (wroteSomething == false) { m->mothurOut("Your file contains only sequences from the .accnos file."); m->mothurOutEndLine(); } outputTypes["count"].push_back(outputFileName); outputNames.push_back(outputFileName); @@ -556,10 +661,10 @@ int RemoveSeqsCommand::readList(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(listfile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(listfile)) + getOutputFileNameTag("list", listfile); - ofstream out; - m->openOutputFile(outputFileName, out); - + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(listfile)); + variables["[extension]"] = m->getExtension(listfile); + ifstream in; m->openInputFile(listfile, in); @@ -576,48 +681,60 @@ int RemoveSeqsCommand::readList(){ //make a new list vector ListVector newList; newList.setLabel(list.getLabel()); + + variables["[distance]"] = list.getLabel(); + string outputFileName = getOutputFileName("list", variables); + ofstream out; + m->openOutputFile(outputFileName, out); + outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName); + + vector binLabels = list.getLabels(); + vector newBinLabels; + + if (m->control_pressed) { in.close(); out.close(); return 0; } + //for each bin for (int i = 0; i < list.getNumBins(); i++) { if (m->control_pressed) { in.close(); out.close(); m->mothurRemove(outputFileName); return 0; } //parse out names that are in accnos file - string binnames = list.get(i); + string bin = list.get(i); + vector bnames; + m->splitAtComma(bin, bnames); string newNames = ""; - while (binnames.find_first_of(',') != -1) { - string name = binnames.substr(0,binnames.find_first_of(',')); - binnames = binnames.substr(binnames.find_first_of(',')+1, binnames.length()); - - //if that name is in the .accnos file, add it + for (int j = 0; j < bnames.size(); j++) { + string name = bnames[j]; + //if that name is in the .accnos file, add it if (names.count(name) == 0) { newNames += name + ","; } else { removedCount++; } - } - - //get last name - if (names.count(binnames) == 0) { newNames += binnames + ","; } - else { removedCount++; } + } //if there are names in this bin add to new list if (newNames != "") { newNames = newNames.substr(0, newNames.length()-1); //rip off extra comma - newList.push_back(newNames); + newList.push_back(newNames); + newBinLabels.push_back(binLabels[i]); } } //print new listvector if (newList.getNumBins() != 0) { wroteSomething = true; + newList.setLabels(newBinLabels); + newList.printHeaders(out); newList.print(out); + } m->gobble(in); + out.close(); } in.close(); - out.close(); + if (wroteSomething == false) { m->mothurOut("Your file contains only sequences from the .accnos file."); m->mothurOutEndLine(); } - outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName); m->mothurOut("Removed " + toString(removedCount) + " sequences from your list file."); m->mothurOutEndLine(); @@ -634,7 +751,10 @@ int RemoveSeqsCommand::readName(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(namefile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(namefile)) + getOutputFileNameTag("name", namefile); + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(namefile)); + variables["[extension]"] = m->getExtension(namefile); + string outputFileName = getOutputFileName("name", variables); ofstream out; m->openOutputFile(outputFileName, out); @@ -686,6 +806,8 @@ int RemoveSeqsCommand::readName(){ wroteSomething = true; out << validSecond[0] << '\t'; + //we are changing the unique name in the fasta file + uniqueMap[firstCol] = validSecond[0]; //you know you have at least one valid second since first column is valid for (int i = 0; i < validSecond.size()-1; i++) { out << validSecond[i] << ','; } @@ -716,7 +838,10 @@ int RemoveSeqsCommand::readGroup(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(groupfile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(groupfile)) + getOutputFileNameTag("group", groupfile); + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(groupfile)); + variables["[extension]"] = m->getExtension(groupfile); + string outputFileName = getOutputFileName("group", variables); ofstream out; m->openOutputFile(outputFileName, out); @@ -762,7 +887,10 @@ int RemoveSeqsCommand::readTax(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(taxfile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(taxfile)) + getOutputFileNameTag("taxonomy", taxfile); + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(taxfile)); + variables["[extension]"] = m->getExtension(taxfile); + string outputFileName = getOutputFileName("taxonomy", variables); ofstream out; m->openOutputFile(outputFileName, out); @@ -776,12 +904,18 @@ int RemoveSeqsCommand::readTax(){ while(!in.eof()){ if (m->control_pressed) { in.close(); out.close(); m->mothurRemove(outputFileName); return 0; } - in >> name; //read from first column + in >> name; m->gobble(in); //read from first column in >> tax; //read from second column + if (!dups) {//adjust name if needed + map::iterator it = uniqueMap.find(name); + if (it != uniqueMap.end()) { name = it->second; } + } + //if this name is in the accnos file if (names.count(name) == 0) { wroteSomething = true; + out << name << '\t' << tax << endl; }else { removedCount++; } @@ -808,7 +942,9 @@ int RemoveSeqsCommand::readAlign(){ try { string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(alignfile); } - string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(alignfile)) + getOutputFileNameTag("alignreport"); + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(alignfile)); + string outputFileName = getOutputFileName("alignreport", variables); ofstream out; m->openOutputFile(outputFileName, out); @@ -831,6 +967,11 @@ int RemoveSeqsCommand::readAlign(){ if (m->control_pressed) { in.close(); out.close(); m->mothurRemove(outputFileName); return 0; } in >> name; //read from first column + + if (!dups) {//adjust name if needed + map::iterator it = uniqueMap.find(name); + if (it != uniqueMap.end()) { name = it->second; } + } //if this name is in the accnos file if (names.count(name) == 0) {