X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=preclustercommand.h;h=0bd56c8662a5b8c622b57a4628f23cc829828f90;hb=55bbd10379db27def51cec72a8819d775f73e45b;hp=2868004825b8336ab1c4e4feb25dc415103b1881;hpb=192b4636345c51d962d4711206535e34cb2fd97c;p=mothur.git diff --git a/preclustercommand.h b/preclustercommand.h index 2868004..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 { @@ -22,12 +23,19 @@ struct seqPNode { Sequence seq; string names; bool active; + int diffs; seqPNode() {} - seqPNode(int n, Sequence s, string nm) : numIdentical(n), seq(s), names(nm), active(1) {} + seqPNode(int n, Sequence s, string nm) : numIdentical(n), seq(s), names(nm), active(1) { diffs = 0; } ~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 { @@ -40,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"; } @@ -56,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 @@ -71,46 +85,50 @@ private: void readNameFile(); //int readNamesFASTA(); int calcMisMatches(string, string); - void printData(string, string); //fasta filename, names file name - int process(); - int loadSeqs(map&, vector&); - int driverGroups(SequenceParser*, string, string, int, int, vector groups); - int createProcessesGroups(SequenceParser*, string, string, vector); + void printData(string, string, string); //fasta filename, names file name + int process(string); + 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); }; /**************************************************************************************************/ //custom data structure for threads to use. // This is passed by void pointer so it can be any data type // that can be passed using a single void pointer (LPVOID). -typedef struct preClusterData { +struct preClusterData { string fastafile; string namefile; - string groupfile; - string newFName, newNName; + string groupfile, countfile; + string newFName, newNName, newMName; MothurOut* m; int start; int end; int diffs, threadID; vector groups; + vector mapFileNames; preClusterData(){} - preClusterData(string f, string n, string g, string nff, string nnf, 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; newFName = nff; newNName = nnf; + newMName = nmf; m = mout; start = st; end = en; diffs = d; threadID = tid; groups = gr; + countfile = c; } }; /**************************************************************************************************/ -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) #else static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ preClusterData* pDataArray; @@ -120,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(); @@ -139,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. //////////////////////////////////////////////////// @@ -150,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++) { @@ -172,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(); } } } @@ -193,6 +232,10 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ //////////////////////////////////////////////////// //int count = process(); - same function below + ofstream out; + pDataArray->m->openOutputFile(pDataArray->newMName+pDataArray->groups[k]+".map", out); + pDataArray->mapFileNames.push_back(pDataArray->newMName+pDataArray->groups[k]+".map"); + //sort seqs by number of identical seqs sort(alignSeqs.begin(), alignSeqs.end(), comparePriority); @@ -206,6 +249,8 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ if (alignSeqs[i].active) { //this sequence has not been merged yet + string chunk = alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(0) + "\t" + alignSeqs[i].seq.getAligned() + "\n"; + //try to merge it with all smaller seqs for (int j = i+1; j < numSeqs; j++) { @@ -229,18 +274,22 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ alignSeqs[j].active = 0; alignSeqs[j].numIdentical = 0; + alignSeqs[j].diffs = mismatch; count++; + chunk += alignSeqs[j].seq.getName() + "\t" + toString(alignSeqs[j].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[j].seq.getAligned() + "\n"; } }//end if j active - }//end if i != j + }//end for loop j //remove from active list alignSeqs[i].active = 0; + out << "ideal_seq_" << (i+1) << '\t' << alignSeqs[i].numIdentical << endl << chunk << endl; + }//end if active i if(i % 100 == 0) { pDataArray->m->mothurOut(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); pDataArray->m->mothurOutEndLine(); } } - + out.close(); if(numSeqs % 100 != 0) { pDataArray->m->mothurOut(toString(numSeqs) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); pDataArray->m->mothurOutEndLine(); } //////////////////////////////////////////////////// @@ -260,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; } + } } @@ -277,7 +328,7 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ } catch(exception& e) { - pDataArray->m->errorOut(e, "AlignCommand", "MyPreclusterThreadFunction"); + pDataArray->m->errorOut(e, "PreClusterCommand", "MyPreclusterThreadFunction"); exit(1); } }