X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=preclustercommand.h;h=2b93328e64a038695a3abbd10cfe9be1165b783d;hp=084bdc61ee464c665e2f37e200b38c09b85e8ef6;hb=b206f634aae1b4ce13978d203247fb64757d5482;hpb=49d2b7459c5027557564b21e9487dadafbbbdc96 diff --git a/preclustercommand.h b/preclustercommand.h index 084bdc6..2b93328 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,21 @@ struct seqPNode { ~seqPNode() {} }; /************************************************************/ -inline bool comparePriority(seqPNode first, seqPNode second) { return (first.numIdentical > second.numIdentical); } +inline bool comparePriorityTopDown(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; +} +/************************************************************/ +inline bool comparePriorityDownTop(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,9 +56,10 @@ public: vector setParameters(); string getCommandName() { return "pre.cluster"; } string getCommandCategory() { return "Sequence Processing"; } - string getOutputFileNameTag(string, string); + 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"; } @@ -58,26 +74,30 @@ 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; + bool abort, bygroup, topdown; + 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 map::iterator itSize; // map active; //maps sequence name to whether it has already been merged or not. vector outputNames; - map > outputTypes; int readFASTA(); 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); }; /**************************************************************************************************/ @@ -87,17 +107,18 @@ private: struct preClusterData { string fastafile; string namefile; - string groupfile; + string groupfile, countfile; string newFName, newNName, newMName; MothurOut* m; int start; - int end; + int end, count; int diffs, threadID; vector groups; vector mapFileNames; + bool topdown; 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, bool td, int tid) { fastafile = f; namefile = n; groupfile = g; @@ -110,6 +131,9 @@ struct preClusterData { diffs = d; threadID = tid; groups = gr; + countfile = c; + topdown = td; + count=0; } }; @@ -124,10 +148,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(); @@ -136,6 +165,8 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ //precluster each group for (int k = pDataArray->start; k < pDataArray->end; k++) { + pDataArray->count++; + int start = time(NULL); if (pDataArray->m->control_pressed) { delete parser; return 0; } @@ -143,8 +174,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. //////////////////////////////////////////////////// @@ -154,6 +190,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++) { @@ -176,8 +215,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(); } } } @@ -201,61 +248,114 @@ static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){ 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); - + //sort seqs by number of identical seqs + if (pDataArray->topdown) { sort(alignSeqs.begin(), alignSeqs.end(), comparePriorityTopDown); } + else { sort(alignSeqs.begin(), alignSeqs.end(), comparePriorityDownTop); } + int count = 0; - //think about running through twice... - for (int i = 0; i < numSeqs; i++) { - - //are you active - // itActive = active.find(alignSeqs[i].seq.getName()); - - 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"; + if (pDataArray->topdown) { + //think about running through twice... + for (int i = 0; i < numSeqs; i++) { + + //are you active + // itActive = active.find(alignSeqs[i].seq.getName()); + + 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++) { - - if (pDataArray->m->control_pressed) { delete parser; return 0; } - - if (alignSeqs[j].active) { //this sequence has not been merged yet - //are you within "diff" bases - //int mismatch = calcMisMatches(alignSeqs[i].seq.getAligned(), alignSeqs[j].seq.getAligned()); - int mismatch = 0; - - for (int k = 0; k < alignSeqs[i].seq.getAligned().length(); k++) { - //do they match - if (alignSeqs[i].seq.getAligned()[k] != alignSeqs[j].seq.getAligned()[k]) { mismatch++; } - if (mismatch > pDataArray->diffs) { mismatch = length; break; } //to far to cluster - } - - if (mismatch <= pDataArray->diffs) { - //merge - alignSeqs[i].names += ',' + alignSeqs[j].names; - alignSeqs[i].numIdentical += alignSeqs[j].numIdentical; - - 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 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(); } + //try to merge it with all smaller seqs + for (int j = i+1; j < numSeqs; j++) { + + if (pDataArray->m->control_pressed) { delete parser; return 0; } + + if (alignSeqs[j].active) { //this sequence has not been merged yet + //are you within "diff" bases + //int mismatch = calcMisMatches(alignSeqs[i].seq.getAligned(), alignSeqs[j].seq.getAligned()); + int mismatch = 0; + + for (int k = 0; k < alignSeqs[i].seq.getAligned().length(); k++) { + //do they match + if (alignSeqs[i].seq.getAligned()[k] != alignSeqs[j].seq.getAligned()[k]) { mismatch++; } + if (mismatch > pDataArray->diffs) { mismatch = length; break; } //to far to cluster + } + + if (mismatch <= pDataArray->diffs) { + //merge + alignSeqs[i].names += ',' + alignSeqs[j].names; + alignSeqs[i].numIdentical += alignSeqs[j].numIdentical; + + 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 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->mothurOutJustToScreen(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)+"\n"); } + } + + }else { + map mapFile; + map originalCount; + map::iterator itCount; + for (int i = 0; i < numSeqs; i++) { mapFile[i] = ""; originalCount[i] = alignSeqs[i].numIdentical; } + + //think about running through twice... + for (int i = 0; i < numSeqs; i++) { + + //try to merge it into larger seqs + for (int j = i+1; j < numSeqs; j++) { + + if (pDataArray->m->control_pressed) { out.close(); return 0; } + + if (originalCount[j] > originalCount[i]) { //this sequence is more abundant than I am + //are you within "diff" bases + //int mismatch = calcMisMatches(alignSeqs[i].seq.getAligned(), alignSeqs[j].seq.getAligned()); + int mismatch = 0; + + for (int k = 0; k < alignSeqs[i].seq.getAligned().length(); k++) { + //do they match + if (alignSeqs[i].seq.getAligned()[k] != alignSeqs[j].seq.getAligned()[k]) { mismatch++; } + if (mismatch > pDataArray->diffs) { mismatch = length; break; } //to far to cluster + } + + if (mismatch <= pDataArray->diffs) { + //merge + alignSeqs[j].names += ',' + alignSeqs[i].names; + alignSeqs[j].numIdentical += alignSeqs[i].numIdentical; + + mapFile[j] = alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[i].seq.getAligned() + "\n" + mapFile[i]; + alignSeqs[i].numIdentical = 0; + originalCount.erase(i); + mapFile[i] = ""; + count++; + j+=numSeqs; //exit search, we merged this one in. + } + }//end abundance check + }//end for loop j + + if(i % 100 == 0) { pDataArray->m->mothurOutJustToScreen(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)+"\n"); } + } + + for (int i = 0; i < numSeqs; i++) { + if (alignSeqs[i].numIdentical != 0) { + out << "ideal_seq_" << (i+1) << '\t' << alignSeqs[i].numIdentical << endl << alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(0) + "\t" + alignSeqs[i].seq.getAligned() + "\n" << mapFile[i] << endl; + } + } + + } + out.close(); + if(numSeqs % 100 != 0) { pDataArray->m->mothurOut(toString(numSeqs) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); pDataArray->m->mothurOutEndLine(); } //////////////////////////////////////////////////// if (pDataArray->m->control_pressed) { delete parser; return 0; } @@ -274,7 +374,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; } + } }