X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=preclustercommand.cpp;h=05c8b7d42a1ba7de160449f02707097800357bea;hp=dadc9186f243952445db24766cf3a08c9dc4d5d0;hb=b206f634aae1b4ce13978d203247fb64757d5482;hpb=96dbe925073caefaed6e6db85659c144a806aeb1 diff --git a/preclustercommand.cpp b/preclustercommand.cpp index dadc918..05c8b7d 100644 --- a/preclustercommand.cpp +++ b/preclustercommand.cpp @@ -13,14 +13,16 @@ //********************************************************************************************************************** vector PreClusterCommand::setParameters(){ try { - CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta); - CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none",false,false); parameters.push_back(pname); - CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none",false,false); parameters.push_back(pcount); - CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none",false,false); parameters.push_back(pgroup); - CommandParameter pdiffs("diffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(pdiffs); - CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); - CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); - CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); + CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","fasta-name",false,true,true); parameters.push_back(pfasta); + 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", "", "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); + + CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir); + CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir); vector myArray; for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } @@ -37,11 +39,12 @@ string PreClusterCommand::getHelpString(){ string helpString = ""; helpString += "The pre.cluster command groups sequences that are within a given number of base mismatches.\n"; helpString += "The pre.cluster command outputs a new fasta and name file.\n"; - helpString += "The pre.cluster command parameters are fasta, name, group, count, processors and diffs. The fasta parameter is required. \n"; + helpString += "The pre.cluster command parameters are fasta, name, group, count, topdown, processors and diffs. The fasta parameter is required. \n"; helpString += "The name parameter allows you to give a list of seqs that are identical. This file is 2 columns, first column is name or representative sequence, second column is a list of its identical sequences separated by commas.\n"; 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, 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"; @@ -54,27 +57,22 @@ string PreClusterCommand::getHelpString(){ } } //********************************************************************************************************************** -string PreClusterCommand::getOutputFileNameTag(string type, string inputName=""){ - try { - string outputFileName = ""; - map >::iterator it; +string PreClusterCommand::getOutputPattern(string type) { + try { + string pattern = ""; - //is this a type this command creates - it = outputTypes.find(type); - if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); } - else { - if (type == "fasta") { outputFileName = "precluster" + m->getExtension(inputName); } - else if (type == "name") { outputFileName = "precluster.names"; } - else if (type == "count") { outputFileName = "precluster.count_table"; } - else if (type == "map") { outputFileName = "precluster.map"; } - else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; } - } - return outputFileName; - } - catch(exception& e) { - m->errorOut(e, "PreClusterCommand", "getOutputFileNameTag"); - exit(1); - } + if (type == "fasta") { pattern = "[filename],precluster,[extension]"; } + else if (type == "name") { pattern = "[filename],precluster.names"; } + else if (type == "count") { pattern = "[filename],precluster.count_table"; } + else if (type == "map") { pattern = "[filename],precluster.map"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } + + return pattern; + } + catch(exception& e) { + m->errorOut(e, "PreClusterCommand", "getOutputPattern"); + exit(1); + } } //********************************************************************************************************************** PreClusterCommand::PreClusterCommand(){ @@ -194,7 +192,7 @@ PreClusterCommand::PreClusterCommand(string option) { else if (countfile == "not open") { abort = true; countfile = ""; } else { m->setCountTableFile(countfile); - ct.readTable(countfile); + ct.readTable(countfile, true, false); if (ct.hasGroupInfo()) { bygroup = true; } else { bygroup = false; } } @@ -215,6 +213,9 @@ PreClusterCommand::PreClusterCommand(string option) { m->setProcessors(temp); m->mothurConvert(temp, processors); + temp = validParameter.validFile(parameters, "topdown", false); if(temp == "not found"){ temp = "T"; } + topdown = m->isTrue(temp); + if (countfile == "") { if (namefile == "") { vector files; files.push_back(fastafile); @@ -239,10 +240,13 @@ int PreClusterCommand::execute(){ int start = time(NULL); string fileroot = outputDir + m->getRootName(m->getSimpleName(fastafile)); - string newFastaFile = fileroot + getOutputFileNameTag("fasta", fastafile); - string newNamesFile = fileroot + getOutputFileNameTag("name"); - string newCountFile = fileroot + getOutputFileNameTag("count"); - string newMapFile = fileroot + getOutputFileNameTag("map"); //add group name if by group + map variables; + variables["[filename]"] = fileroot; + string newNamesFile = getOutputFileName("name",variables); + string newCountFile = getOutputFileName("count",variables); + string newMapFile = getOutputFileName("map",variables); //add group name if by group + variables["[extension]"] = m->getExtension(fastafile); + string newFastaFile = getOutputFileName("fasta", variables); outputNames.push_back(newFastaFile); outputTypes["fasta"].push_back(newFastaFile); if (countfile == "") { outputNames.push_back(newNamesFile); outputTypes["name"].push_back(newNamesFile); } else { outputNames.push_back(newCountFile); outputTypes["count"].push_back(newCountFile); } @@ -296,6 +300,7 @@ int PreClusterCommand::execute(){ m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run pre.cluster."); m->mothurOutEndLine(); }else { + if (processors != 1) { m->mothurOut("When using running without group information mothur can only use 1 processor, continuing."); m->mothurOutEndLine(); processors = 1; } if (namefile != "") { readNameFile(); } //reads fasta file and return number of seqs @@ -364,13 +369,15 @@ int PreClusterCommand::createProcessesGroups(string newFName, string newNName, s //divide the groups between the processors vector lines; - int numGroupsPerProcessor = groups.size() / processors; - for (int i = 0; i < processors; i++) { - int startIndex = i * numGroupsPerProcessor; - int endIndex = (i+1) * numGroupsPerProcessor; - if(i == (processors - 1)){ endIndex = groups.size(); } - lines.push_back(linePair(startIndex, endIndex)); - } + int remainingPairs = groups.size(); + int startIndex = 0; + for (int remainingProcessors = processors; remainingProcessors > 0; remainingProcessors--) { + int numPairs = remainingPairs; //case for last processor + if (remainingProcessors != 1) { numPairs = ceil(remainingPairs / remainingProcessors); } + lines.push_back(linePair(startIndex, (startIndex+numPairs))); //startIndex, endIndex + startIndex = startIndex + numPairs; + remainingPairs = remainingPairs - numPairs; + } #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) @@ -442,7 +449,7 @@ int PreClusterCommand::createProcessesGroups(string newFName, string newNName, s // Allocate memory for thread data. string extension = toString(i) + ".temp"; - preClusterData* tempPreCluster = new preClusterData(fastafile, namefile, groupfile, countfile, (newFName+extension), (newNName+extension), newMFile, groups, m, lines[i].start, lines[i].end, diffs, i); + preClusterData* tempPreCluster = new preClusterData(fastafile, namefile, groupfile, countfile, (newFName+extension), (newNName+extension), newMFile, groups, m, lines[i].start, lines[i].end, diffs, topdown, i); pDataArray.push_back(tempPreCluster); processIDS.push_back(i); @@ -460,6 +467,9 @@ int PreClusterCommand::createProcessesGroups(string newFName, string newNName, s //Close all thread handles and free memory allocations. for(int i=0; i < pDataArray.size(); i++){ + if (pDataArray[i]->count != (pDataArray[i]->end-pDataArray[i]->start)) { + m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->end-pDataArray[i]->start) + " groups assigned to it, quitting. \n"); m->control_pressed = true; + } for (int j = 0; j < pDataArray[i]->mapFileNames.size(); j++) { outputNames.push_back(pDataArray[i]->mapFileNames[j]); outputTypes["map"].push_back(pDataArray[i]->mapFileNames[j]); } @@ -547,52 +557,94 @@ int PreClusterCommand::process(string newMapFile){ m->openOutputFile(newMapFile, out); //sort seqs by number of identical seqs - sort(alignSeqs.begin(), alignSeqs.end(), comparePriority); + if (topdown) { sort(alignSeqs.begin(), alignSeqs.end(), comparePriorityTopDown); } + else { sort(alignSeqs.begin(), alignSeqs.end(), comparePriorityDownTop); } 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(); } @@ -853,13 +905,13 @@ void PreClusterCommand::readNameFile(){ while (!in.eof()) { in >> firstCol >> secondCol; m->gobble(in); + + m->checkName(firstCol); + m->checkName(secondCol); + int size = m->getNumNames(secondCol); + names[firstCol] = secondCol; - int size = 1; - - for(int i=0;i