X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=preclustercommand.cpp;h=bbf8df27ac5ee691f091f150638b781f85486a11;hp=a1d5f48bdd77c81bb12d0095628922202a4bcacd;hb=c48d91112209b841444923670dca5454da0e2a4d;hpb=c85db0a4be3a1f8037a71a23ca73f9762184e28a diff --git a/preclustercommand.cpp b/preclustercommand.cpp index a1d5f48..bbf8df2 100644 --- a/preclustercommand.cpp +++ b/preclustercommand.cpp @@ -17,7 +17,7 @@ vector PreClusterCommand::setParameters(){ CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","name",false,false,true); parameters.push_back(pname); CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","count",false,false,true); parameters.push_back(pcount); CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","",false,false,true); parameters.push_back(pgroup); - CommandParameter pdiffs("diffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(pdiffs); + CommandParameter pdiffs("diffs", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pdiffs); CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors); CommandParameter ptopdown("topdown", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(ptopdown); @@ -44,7 +44,7 @@ string PreClusterCommand::getHelpString(){ helpString += "The group parameter allows you to provide a group file so you can cluster by group. \n"; helpString += "The count parameter allows you to provide a count file so you can cluster by group. \n"; helpString += "The diffs parameter allows you to specify maximum number of mismatched bases allowed between sequences in a grouping. The default is 1.\n"; - helpString += "The topdown parameter allows you to specify whether to cluster from largest abundance to smallest or smallest to largest. Default=T, meanging largest to smallest.\n"; + helpString += "The topdown parameter allows you to specify whether to cluster from largest abundance to smallest or smallest to largest. Default=T, meaning largest to smallest.\n"; helpString += "The pre.cluster command should be in the following format: \n"; helpString += "pre.cluster(fasta=yourFastaFile, names=yourNamesFile, diffs=yourMaxDiffs) \n"; helpString += "Example pre.cluster(fasta=amazon.fasta, diffs=2).\n"; @@ -560,47 +560,88 @@ int PreClusterCommand::process(string newMapFile){ int count = 0; int numSeqs = alignSeqs.size(); - //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 (m->control_pressed) { out.close(); 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()); - - if (mismatch <= diffs) { - //merge - alignSeqs[i].names += ',' + alignSeqs[j].names; - alignSeqs[i].numIdentical += alignSeqs[j].numIdentical; - - chunk += alignSeqs[j].seq.getName() + "\t" + toString(alignSeqs[j].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[j].seq.getAligned() + "\n"; - - alignSeqs[j].active = 0; - alignSeqs[j].numIdentical = 0; - count++; - } - }//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) { m->mothurOut(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); m->mothurOutEndLine(); } - } + if (topdown) { + //think about running through twice... + for (int i = 0; i < numSeqs; i++) { + + 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 (m->control_pressed) { out.close(); 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()); + + if (mismatch <= diffs) { + //merge + alignSeqs[i].names += ',' + alignSeqs[j].names; + alignSeqs[i].numIdentical += alignSeqs[j].numIdentical; + + chunk += alignSeqs[j].seq.getName() + "\t" + toString(alignSeqs[j].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[j].seq.getAligned() + "\n"; + + alignSeqs[j].active = 0; + alignSeqs[j].numIdentical = 0; + count++; + } + }//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) { 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 (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()); + + if (mismatch <= 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) { 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) { m->mothurOut(toString(numSeqs) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); m->mothurOutEndLine(); } @@ -861,13 +902,19 @@ void PreClusterCommand::readNameFile(){ while (!in.eof()) { in >> firstCol >> secondCol; m->gobble(in); + + for (int i = 0; i < firstCol.length(); i++) { + if (firstCol[i] == ':') { firstCol[i] = '_'; m->changedSeqNames = true; } + } + + int size = 1; + for (int i = 0; i < secondCol.length(); i++) { + if (secondCol[i] == ':') { secondCol[i] = '_'; m->changedSeqNames = true; } + else if(secondCol[i] == ','){ size++; } + } + names[firstCol] = secondCol; - int size = 1; - - for(int i=0;i