X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=preclustercommand.h;h=2b93328e64a038695a3abbd10cfe9be1165b783d;hp=9844b8ffc1a951c47d606ec23e7caf45b4bf4ce1;hb=b206f634aae1b4ce13978d203247fb64757d5482;hpb=96dbe925073caefaed6e6db85659c144a806aeb1 diff --git a/preclustercommand.h b/preclustercommand.h index 9844b8f..2b93328 100644 --- a/preclustercommand.h +++ b/preclustercommand.h @@ -29,13 +29,21 @@ struct seqPNode { ~seqPNode() {} }; /************************************************************/ -inline bool comparePriority(seqPNode first, seqPNode second) { +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 { @@ -48,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"; } @@ -70,7 +79,7 @@ private: CountTable ct; int diffs, length, processors; - bool abort, bygroup; + 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 @@ -78,7 +87,6 @@ private: 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(); @@ -103,13 +111,14 @@ struct preClusterData { 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 c, 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; @@ -123,6 +132,8 @@ struct preClusterData { threadID = tid; groups = gr; countfile = c; + topdown = td; + count=0; } }; @@ -154,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; } @@ -235,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; }