X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=preclustercommand.h;h=0bd56c8662a5b8c622b57a4628f23cc829828f90;hb=deba0af0ccdcb6005ed5b2b82649b137c63fbdf7;hp=64efc7fc1d43aca03d17497e28bee2fdb61cc74a;hpb=91a27e0483827c06c21c4fe89558923bbfe86573;p=mothur.git diff --git a/preclustercommand.h b/preclustercommand.h index 64efc7f..0bd56c8 100644 --- a/preclustercommand.h +++ b/preclustercommand.h @@ -15,6 +15,7 @@ #include "command.hpp" #include "sequence.hpp" #include "sequenceparser.h" +#include "sequencecountparser.h" /************************************************************/ struct seqPNode { @@ -28,7 +29,13 @@ struct seqPNode { ~seqPNode() {} }; /************************************************************/ -inline bool comparePriority(seqPNode first, seqPNode second) { return (first.numIdentical > second.numIdentical); } +inline bool comparePriority(seqPNode first, seqPNode second) { + if (first.numIdentical > second.numIdentical) { return true; } + else if (first.numIdentical == second.numIdentical) { + if (first.seq.getName() > second.seq.getName()) { return true; } + } + return false; +} //************************************************************/ class PreClusterCommand : public Command { @@ -41,8 +48,10 @@ public: vector setParameters(); string getCommandName() { return "pre.cluster"; } string getCommandCategory() { return "Sequence Processing"; } + string getHelpString(); - string getCitation() { return "http://www.mothur.org/wiki/Pre.cluster"; } + string getOutputPattern(string); + string getCitation() { return "Schloss PD, Gevers D, Westcott SL (2011). Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies. PLoS ONE. 6:e27310.\nhttp://www.mothur.org/wiki/Pre.cluster"; } string getDescription() { return "implements a pseudo-single linkage algorithm with the goal of removing sequences that are likely due to pyrosequencing errors"; } @@ -57,9 +66,13 @@ private: linePair(int i, int j) : start(i), end(j) {} }; + SequenceParser* parser; + SequenceCountParser* cparser; + CountTable ct; + int diffs, length, processors; bool abort, bygroup; - string fastafile, namefile, outputDir, groupfile; + string fastafile, namefile, outputDir, groupfile, countfile; vector alignSeqs; //maps the number of identical seqs to a sequence map names; //represents the names file first column maps to second column map sizes; //this map a seq name to the number of identical seqs in the names file @@ -72,11 +85,12 @@ private: void readNameFile(); //int readNamesFASTA(); int calcMisMatches(string, string); - void printData(string, string); //fasta filename, names file name + void printData(string, string, string); //fasta filename, names file name int process(string); - int loadSeqs(map&, vector&); - int driverGroups(SequenceParser*, string, string, string, int, int, vector groups); - int createProcessesGroups(SequenceParser*, string, string, string, vector); + int loadSeqs(map&, vector&, string); + int driverGroups(string, string, string, int, int, vector groups); + int createProcessesGroups(string, string, string, vector); + int mergeGroupCounts(string, string, string); }; /**************************************************************************************************/ @@ -86,7 +100,7 @@ private: struct preClusterData { string fastafile; string namefile; - string groupfile; + string groupfile, countfile; string newFName, newNName, newMName; MothurOut* m; int start; @@ -96,7 +110,7 @@ struct preClusterData { vector mapFileNames; preClusterData(){} - preClusterData(string f, string n, string g, string nff, string nnf, string nmf, vector gr, MothurOut* mout, int st, int en, int d, int tid) { + preClusterData(string f, string n, string g, string c, string nff, string nnf, string nmf, vector gr, MothurOut* mout, int st, int en, int d, int tid) { fastafile = f; namefile = n; groupfile = g; @@ -109,6 +123,7 @@ struct preClusterData { diffs = d; threadID = tid; groups = gr; + countfile = c; } }; @@ -123,10 +138,15 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ //parse fasta and name file by group SequenceParser* parser; - if (pDataArray->namefile != "") { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile, pDataArray->namefile); } - else { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile); } - - int numSeqs = 0; + SequenceCountParser* cparser; + if (pDataArray->countfile != "") { + cparser = new SequenceCountParser(pDataArray->countfile, pDataArray->fastafile); + }else { + if (pDataArray->namefile != "") { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile, pDataArray->namefile); } + else { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile); } + } + + int numSeqs = 0; vector alignSeqs; //clear out old files ofstream outF; pDataArray->m->openOutputFile(pDataArray->newFName, outF); outF.close(); @@ -142,8 +162,13 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ pDataArray->m->mothurOutEndLine(); pDataArray->m->mothurOut("Processing group " + pDataArray->groups[k] + ":"); pDataArray->m->mothurOutEndLine(); map thisNameMap; - if (pDataArray->namefile != "") { thisNameMap = parser->getNameMap(pDataArray->groups[k]); } - vector thisSeqs = parser->getSeqs(pDataArray->groups[k]); + vector thisSeqs; + if (pDataArray->groupfile != "") { + thisSeqs = parser->getSeqs(pDataArray->groups[k]); + }else if (pDataArray->countfile != "") { + thisSeqs = cparser->getSeqs(pDataArray->groups[k]); + } + if (pDataArray->namefile != "") { thisNameMap = parser->getNameMap(pDataArray->groups[k]); } //fill alignSeqs with this groups info. //////////////////////////////////////////////////// @@ -153,6 +178,9 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ alignSeqs.clear(); map::iterator it; bool error = false; + map thisCount; + if (pDataArray->countfile != "") { thisCount = cparser->getCountTable(pDataArray->groups[k]); } + for (int i = 0; i < thisSeqs.size(); i++) { @@ -175,8 +203,16 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ if (thisSeqs[i].getAligned().length() > length) { length = thisSeqs[i].getAligned().length(); } } }else { //no names file, you are identical to yourself - seqPNode tempNode(1, thisSeqs[i], thisSeqs[i].getName()); - alignSeqs.push_back(tempNode); + int numRep = 1; + if (pDataArray->countfile != "") { + map::iterator it2 = thisCount.find(thisSeqs[i].getName()); + + //should never be true since parser checks for this + if (it2 == thisCount.end()) { pDataArray->m->mothurOut(thisSeqs[i].getName() + " is not in your count file, please correct."); pDataArray->m->mothurOutEndLine(); error = true; } + else { numRep = it2->second; } + } + seqPNode tempNode(numRep, thisSeqs[i], thisSeqs[i].getName()); + alignSeqs.push_back(tempNode); if (thisSeqs[i].getAligned().length() > length) { length = thisSeqs[i].getAligned().length(); } } } @@ -273,7 +309,9 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ for (int i = 0; i < alignSeqs.size(); i++) { if (alignSeqs[i].numIdentical != 0) { alignSeqs[i].seq.printSequence(outFasta); - outNames << alignSeqs[i].seq.getName() << '\t' << alignSeqs[i].names << endl; + if (pDataArray->countfile != "") { outNames << pDataArray->groups[k] << '\t' << alignSeqs[i].seq.getName() << '\t' << alignSeqs[i].names << endl; + }else { outNames << alignSeqs[i].seq.getName() << '\t' << alignSeqs[i].names << endl; } + } }