From bc4c4d63e983a3e91e74b5038bfec6705e1a3a2e Mon Sep 17 00:00:00 2001 From: pschloss Date: Tue, 6 Jul 2010 20:30:15 +0000 Subject: [PATCH] sens.spec changes --- Mothur.xcodeproj/project.pbxproj | 4 + commandfactory.cpp | 3 + makefile | 1 + sensspeccommand.cpp | 316 +++++++++++++++++++++++++++++++ sensspeccommand.h | 39 ++++ trimseqscommand.cpp | 1 - 6 files changed, 363 insertions(+), 1 deletion(-) create mode 100644 sensspeccommand.cpp create mode 100644 sensspeccommand.h diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 6e28507..6770cec 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -7,6 +7,8 @@ objects = { /* Begin PBXFileReference section */ + 7EA299BA11E384940022D8D3 /* sensspeccommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sensspeccommand.h; sourceTree = ""; }; + 7EA299BB11E384940022D8D3 /* sensspeccommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sensspeccommand.cpp; sourceTree = ""; }; 7EC61A0911BEA6AF00F668D9 /* weightedlinkage.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = weightedlinkage.cpp; sourceTree = ""; }; A703FE931194645F002C397E /* makegroupcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = makegroupcommand.h; sourceTree = SOURCE_ROOT; }; A703FE941194645F002C397E /* makegroupcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = makegroupcommand.cpp; sourceTree = SOURCE_ROOT; }; @@ -655,6 +657,8 @@ A7CB593E11402F110010EB83 /* commands */ = { isa = PBXGroup; children = ( + 7EA299BA11E384940022D8D3 /* sensspeccommand.h */, + 7EA299BB11E384940022D8D3 /* sensspeccommand.cpp */, A7DA202B113FECD400BF472F /* command.hpp */, A7DA1FEF113FECD400BF472F /* aligncommand.h */, A7DA1FEE113FECD400BF472F /* aligncommand.cpp */, diff --git a/commandfactory.cpp b/commandfactory.cpp index 2b19f72..725ce8c 100644 --- a/commandfactory.cpp +++ b/commandfactory.cpp @@ -82,6 +82,7 @@ #include "classifyotucommand.h" #include "degapseqscommand.h" #include "getrelabundcommand.h" +#include "sensspeccommand.h" /*******************************************************/ @@ -183,6 +184,7 @@ CommandFactory::CommandFactory(){ commands["screen.seqs"] = "MPIEnabled"; commands["summary.seqs"] = "MPIEnabled"; commands["cluster.split"] = "MPIEnabled"; + commands["sens.spec"] = "sens.spec"; commands["quit"] = "MPIEnabled"; } @@ -297,6 +299,7 @@ Command* CommandFactory::getCommand(string commandName, string optionString){ else if(commandName == "classify.otu") { command = new ClassifyOtuCommand(optionString); } else if(commandName == "degap.seqs") { command = new DegapSeqsCommand(optionString); } else if(commandName == "get.relabund") { command = new GetRelAbundCommand(optionString); } + else if(commandName == "sens.spec") { command = new SensSpecCommand(optionString); } else { command = new NoCommand(optionString); } return command; diff --git a/makefile b/makefile index 052519e..af46fd9 100644 --- a/makefile +++ b/makefile @@ -66,3 +66,4 @@ install : mothur clean : @rm -f $(OBJECTS) + diff --git a/sensspeccommand.cpp b/sensspeccommand.cpp new file mode 100644 index 0000000..0416dcb --- /dev/null +++ b/sensspeccommand.cpp @@ -0,0 +1,316 @@ +/* + * sensspeccommand.cpp + * Mothur + * + * Created by Pat Schloss on 7/6/10. + * Copyright 2010 Schloss Lab. All rights reserved. + * + */ + +#include "sensspeccommand.h" + +//*************************************************************************************************************** + +SensSpecCommand::SensSpecCommand(string option) { + try { + + abort = false; + + //allow user to run help + if(option == "help") { help(); abort = true; } + + else { + string temp; + + //valid paramters for this command + string AlignArray[] = {"list", "phylip", "column", "name", "hard", "label", "cutoff", "outputdir", "inputdir"}; + + vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); + + OptionParser parser(option); + map parameters = parser.getParameters(); + + ValidParameters validParameter; + map::iterator it; + + //check to make sure all parameters are valid for command + for (it = parameters.begin(); it != parameters.end(); it++) { + if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; } + } + + //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); + if (inputDir == "not found"){ inputDir = ""; } + else { + string path; + it = parameters.find("list"); + //user has given a template file + if(it != parameters.end()){ + path = hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["list"] = inputDir + it->second; } + } + + it = parameters.find("phylip"); + //user has given a template file + if(it != parameters.end()){ + path = hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["phylip"] = inputDir + it->second; } + } + + it = parameters.find("column"); + //user has given a template file + if(it != parameters.end()){ + path = hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["column"] = inputDir + it->second; } + } + + it = parameters.find("name"); + //user has given a template file + if(it != parameters.end()){ + path = 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 + listFile = validParameter.validFile(parameters, "list", true); + if (listFile == "not found") { m->mothurOut("list is a required parameter for the sens.spec command."); m->mothurOutEndLine(); abort = true; } + else if (listFile == "not open") { abort = true; } + + distFile = validParameter.validFile(parameters, "column", true); + format = "column"; + if(distFile == "not found") { + distFile = validParameter.validFile(parameters, "phylip", true); + format = "phylip"; + } + if(distFile == "not found") { m->mothurOut("either column or phylip are required for the sens.spec command."); m->mothurOutEndLine(); abort = true; } + else if (distFile == "not open") { 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 = ""; + outputDir += hasPath(listFile); //if user entered a file with a path then preserve it + } + + //check for optional parameter and set defaults + // ...at some point should added some additional type checking... + temp = validParameter.validFile(parameters, "hard", false); + if (temp == "not found"){ hard = 0; } + else if(!isTrue(temp)) { hard = 0; } + else if(isTrue(temp)) { hard = 1; } + +// temp = validParameter.validFile(parameters, "name", true); +// if (temp == "not found") { nameFile = ""; } +// else if(temp == "not open") { abort = true; } +// else { nameFile = temp; } +// cout << "name:\t" << nameFile << endl; + + temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found") { temp = ""; } + convert(temp, cutoff); + cout << "cutoff:\t" << cutoff << endl; + +// cutoff = 0.0349; + + lineLabel = validParameter.validFile(parameters, "label", false); if (lineLabel == "not found") { lineLabel = ""; } + + } + + } + catch(exception& e) { + m->errorOut(e, "SensSpecCommand", "SensSpecCommand"); + exit(1); + } +} +//********************************************************************************************************************** + +void SensSpecCommand::help(){ + try { + m->mothurOut("The sens.spec command reads a fastaFile and creates .....\n"); + + + + m->mothurOut("Example sens.spec(...).\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, "SensSpecCommand", "help"); + exit(1); + } +} + + +//*************************************************************************************************************** + +SensSpecCommand::~SensSpecCommand(){ /* do nothing */ } + +//*************************************************************************************************************** + +int SensSpecCommand::execute(){ + try{ + if (abort == true) { return 0; } + + if(format == "phylip") { processPhylip(); } +// else if(format == "column") { processColumn(seqMap); } + + +// string seqList; +// map seqMap; + + + + return 0; + } + catch(exception& e) { + m->errorOut(e, "SensSpecCommand", "execute"); + exit(1); + } +} + +//*************************************************************************************************************** + +void SensSpecCommand::processPhylip(){ + try{ + ifstream inputListFile; + openInputFile(listFile, inputListFile); + + string label; + int numOTUs; + + map seqMap; + string seqList; + + while(inputListFile){ + inputListFile >> label >> numOTUs; + for(int i=0;i> seqList; + int seqListLength = seqList.length(); + string seqName = ""; + for(int j=0;j> pNumSeqs; + if(pNumSeqs != lNumSeqs){ cout << "numSeq mismatch!" << endl; } + + string seqName; + double distance; + vector otuIndices(lNumSeqs, -1); + + truePositives = 0; + falsePositives = 0; + trueNegatives = 0; + falseNegatives = 0; + + + for(int i=0;i> seqName; + otuIndices[i] = seqMap[seqName]; + + for(int j=0;j> distance; + if(distance <= cutoff){ + if(otuIndices[i] == otuIndices[j]){ + truePositives++; + } + else{ + falseNegatives++; + } + } + else{ + if(otuIndices[i] == otuIndices[j]){ + falsePositives++; + } + else{ + trueNegatives++; + } + } + } + } + phylipFile.close(); + + cout << "truePositives:\t" << truePositives << endl; + cout << "trueNegatives:\t" << trueNegatives << endl; + cout << "falsePositives:\t" << falsePositives << endl; + cout << "falseNegatives:\t" << falseNegatives << endl; + } + catch(exception& e) { + m->errorOut(e, "SensSpecCommand", "processPhylip"); + exit(1); + } +} + +//*************************************************************************************************************** + +//void SensSpecCommand::processColumn(map seqMap){ +// +// truePositives = 0; +// falsePositives = 0; +// trueNegatives = numDists; +// falseNegatives = 0; +// +// ifstream columnFile; +// openInputFile(distFile, columnFile); +// +// string seqNameA, seqNameB, oldSeqNameA; +// int otuA, otuB, oldOTUA; +// double distance; +// +// while(columnFile){ +// columnFile >> seqNameA >> seqNameB >> distance; +// +// if(seqNameA == oldSeqNameA) { otuA = oldOTUA; } +// else { otuA = seqMap[seqNameA]; oldOTUA = otuA; } +// +// otuB = seqMap[seqNameB]; +// +// if(distance <= cutoff){ +// if(otuA == otuB){ +// truePositives++; +// } +// else{ +// falseNegatives++; +// } +// trueNegatives--; +// } +// else{ +// if(otuA == otuB){ +// falsePositives++; +// trueNegatives--; +// } +// } +// +// gobble(columnFile); +// } +// columnFile.close(); +// +// cout << "truePositives:\t" << truePositives << endl; +// cout << "trueNegatives:\t" << trueNegatives << endl; +// cout << "falsePositives:\t" << falsePositives << endl; +// cout << "falseNegatives:\t" << falseNegatives << endl; +//} + +//*************************************************************************************************************** diff --git a/sensspeccommand.h b/sensspeccommand.h new file mode 100644 index 0000000..52ab31e --- /dev/null +++ b/sensspeccommand.h @@ -0,0 +1,39 @@ +#ifndef SENSSPECCOMMAND_H +#define SENSSPECCOMMAND_H + + +/* + * sensspeccommand.h + * Mothur + * + * Created by Pat Schloss on 7/6/10. + * Copyright 2010 Schloss Lab. All rights reserved. + * + */ + +#include "mothur.h" +#include "command.hpp" + +class SensSpecCommand : public Command { + +public: + SensSpecCommand(string); + ~SensSpecCommand(); + int execute(); + void help(); + +private: + void processPhylip(); +/// void processColumn(map); + + string listFile, distFile, nameFile, outputDir; + string format; +// int numSeqs, numDists; + int truePositives, falsePositives, trueNegatives, falseNegatives; + bool abort; + bool hard; + string lineLabel; + double cutoff; +}; + +#endif \ No newline at end of file diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp index 1975f9a..cbf5e13 100644 --- a/trimseqscommand.cpp +++ b/trimseqscommand.cpp @@ -378,7 +378,6 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string int success = 1; Sequence currSeq(inFASTA); - string origSeq = currSeq.getUnaligned(); if (origSeq != "") { -- 2.39.2