2 // getmetacommunitycommand.cpp
5 // Created by SarahsWork on 4/9/13.
6 // Copyright (c) 2013 Schloss Lab. All rights reserved.
9 #include "getmetacommunitycommand.h"
10 #include "communitytype.h"
12 #include "validcalculator.h"
13 #include "subsample.h"
15 //**********************************************************************************************************************
16 vector<string> GetMetaCommunityCommand::setParameters(){
18 CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","outputType",false,true); parameters.push_back(pshared);
19 CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
20 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
21 CommandParameter pcalc("calc", "Multiple", "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan-kstest-sharednseqs-ochiai-anderberg-kulczynski-kulczynskicody-lennon-morisitahorn-braycurtis-whittaker-odum-canberra-structeuclidean-structchord-hellinger-manhattan-structpearson-soergel-spearman-structkulczynski-speciesprofile-hamming-structchi2-gower-memchi2-memchord-memeuclidean-mempearson-jsd-rjsd", "rjsd", "", "", "","",false,false,true); parameters.push_back(pcalc);
22 CommandParameter psubsample("subsample", "String", "", "", "", "", "","",false,false); parameters.push_back(psubsample);
23 CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
24 CommandParameter pminpartitions("minpartitions", "Number", "", "5", "", "", "","",false,false,true); parameters.push_back(pminpartitions);
25 CommandParameter pmaxpartitions("maxpartitions", "Number", "", "100", "", "", "","",false,false,true); parameters.push_back(pmaxpartitions);
26 CommandParameter poptimizegap("optimizegap", "Number", "", "3", "", "", "","",false,false,true); parameters.push_back(poptimizegap);
27 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
28 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
29 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
30 CommandParameter pmethod("method", "Multiple", "dmm-kmeans-pam", "dmm", "", "", "","",false,false,true); parameters.push_back(pmethod);
32 vector<string> myArray;
33 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
37 m->errorOut(e, "NewCommand", "setParameters");
41 //**********************************************************************************************************************
42 string GetMetaCommunityCommand::getHelpString(){
44 string helpString = "";
45 helpString += "The get.communitytype command parameters are shared, method, label, groups, minpartitions, maxpartitions, optimizegap and processors. The shared file is required. \n";
46 helpString += "The label parameter is used to analyze specific labels in your input. labels are separated by dashes.\n";
47 helpString += "The groups parameter allows you to specify which of the groups in your shared file you would like analyzed. Group names are separated by dashes.\n";
48 helpString += "The method parameter allows to select the method you would like to use. Options are dmm, kmeans and pam. Default=dmm.\n";
49 helpString += "The calc parameter allows to select the calculator you would like to use to calculate the distance matrix used by the pam method. By default the jsd calculator is used.\n";
50 helpString += "The iters parameter allows you to choose the number of times you would like to run the subsample while calculating the distance matirx for the pam method.\n";
51 helpString += "The subsample parameter allows you to enter the size pergroup of the sample or you can set subsample=T and mothur will use the size of your smallest group while calculating the distance matrix for the pam method.\n";
52 helpString += "The minpartitions parameter is used to .... Default=5.\n";
53 helpString += "The maxpartitions parameter is used to .... Default=10.\n";
54 helpString += "The optimizegap parameter is used to .... Default=3.\n";
55 helpString += "The processors parameter allows you to specify number of processors to use. The default is 1.\n";
56 helpString += "The get.communitytype command should be in the following format: get.communitytype(shared=yourSharedFile).\n";
60 m->errorOut(e, "GetMetaCommunityCommand", "getHelpString");
64 //**********************************************************************************************************************
65 string GetMetaCommunityCommand::getOutputPattern(string type) {
69 if (type == "fit") { pattern = "[filename],[distance],mix.fit"; }
70 else if (type == "relabund") { pattern = "[filename],[distance],[tag],mix.relabund"; }
71 else if (type == "design") { pattern = "[filename],[distance],mix.design"; }
72 else if (type == "matrix") { pattern = "[filename],[distance],[tag],mix.posterior"; }
73 else if (type == "parameters") { pattern = "[filename],[distance],mix.parameters"; }
74 else if (type == "summary") { pattern = "[filename],[distance],mix.summary"; }
75 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
80 m->errorOut(e, "GetMetaCommunityCommand", "getOutputPattern");
84 //**********************************************************************************************************************
85 GetMetaCommunityCommand::GetMetaCommunityCommand(){
87 abort = true; calledHelp = true;
89 vector<string> tempOutNames;
90 outputTypes["fit"] = tempOutNames;
91 outputTypes["relabund"] = tempOutNames;
92 outputTypes["matrix"] = tempOutNames;
93 outputTypes["design"] = tempOutNames;
94 outputTypes["parameters"] = tempOutNames;
95 outputTypes["summary"] = tempOutNames;
98 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
102 //**********************************************************************************************************************
103 GetMetaCommunityCommand::GetMetaCommunityCommand(string option) {
105 abort = false; calledHelp = false;
108 //allow user to run help
109 if(option == "help") { help(); abort = true; calledHelp = true; }
110 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
113 //valid paramters for this command
114 vector<string> myArray = setParameters();
116 OptionParser parser(option);
117 map<string,string> parameters = parser.getParameters();
119 ValidParameters validParameter;
120 map<string,string>::iterator it;
121 //check to make sure all parameters are valid for command
122 for (it = parameters.begin(); it != parameters.end(); it++) {
123 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
126 vector<string> tempOutNames;
127 outputTypes["fit"] = tempOutNames;
128 outputTypes["relabund"] = tempOutNames;
129 outputTypes["matrix"] = tempOutNames;
130 outputTypes["design"] = tempOutNames;
131 outputTypes["parameters"] = tempOutNames;
132 outputTypes["summary"] = tempOutNames;
134 //if the user changes the input directory command factory will send this info to us in the output parameter
135 string inputDir = validParameter.validFile(parameters, "inputdir", false);
136 if (inputDir == "not found"){ inputDir = ""; }
139 it = parameters.find("shared");
140 if(it != parameters.end()){
141 path = m->hasPath(it->second);
142 if (path == "") { parameters["shared"] = inputDir + it->second; }
146 //get shared file, it is required
147 sharedfile = validParameter.validFile(parameters, "shared", true);
148 if (sharedfile == "not open") { sharedfile = ""; abort = true; }
149 else if (sharedfile == "not found") {
150 //if there is a current shared file, use it
151 sharedfile = m->getSharedFile();
152 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
153 else { m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
154 }else { m->setSharedFile(sharedfile); }
156 //if the user changes the output directory command factory will send this info to us in the output parameter
157 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
158 outputDir = m->hasPath(sharedfile); //if user entered a file with a path then preserve it
161 string temp = validParameter.validFile(parameters, "minpartitions", false); if (temp == "not found"){ temp = "5"; }
162 m->mothurConvert(temp, minpartitions);
164 temp = validParameter.validFile(parameters, "maxpartitions", false); if (temp == "not found"){ temp = "10"; }
165 m->mothurConvert(temp, maxpartitions);
167 temp = validParameter.validFile(parameters, "optimizegap", false); if (temp == "not found"){ temp = "3"; }
168 m->mothurConvert(temp, optimizegap);
170 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
171 m->setProcessors(temp);
172 m->mothurConvert(temp, processors);
174 string groups = validParameter.validFile(parameters, "groups", false);
175 if (groups == "not found") { groups = ""; }
176 else { m->splitAtDash(groups, Groups); }
177 m->setGroups(Groups);
179 string label = validParameter.validFile(parameters, "label", false);
180 if (label == "not found") { label = ""; }
182 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
183 else { allLines = 1; }
186 method = validParameter.validFile(parameters, "method", false);
187 if (method == "not found") { method = "dmm"; }
189 if ((method == "dmm") || (method == "kmeans") || (method == "pam")) { }
190 else { m->mothurOut("[ERROR]: " + method + " is not a valid method. Valid algorithms are dmm, kmeans and pam."); m->mothurOutEndLine(); abort = true; }
192 calc = validParameter.validFile(parameters, "calc", false);
193 if (calc == "not found") { calc = "rjsd"; }
195 if (calc == "default") { calc = "rjsd"; }
197 m->splitAtDash(calc, Estimators);
198 if (m->inUsersGroups("citation", Estimators)) {
199 ValidCalculators validCalc; validCalc.printCitations(Estimators);
200 //remove citation from list of calcs
201 for (int i = 0; i < Estimators.size(); i++) { if (Estimators[i] == "citation") { Estimators.erase(Estimators.begin()+i); break; } }
203 if (Estimators.size() != 1) { abort = true; m->mothurOut("[ERROR]: only one calculator is allowed.\n"); }
205 temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
206 m->mothurConvert(temp, iters);
208 temp = validParameter.validFile(parameters, "subsample", false); if (temp == "not found") { temp = "F"; }
209 if (m->isNumeric1(temp)) { m->mothurConvert(temp, subsampleSize); subsample = true; }
211 if (m->isTrue(temp)) { subsample = true; subsampleSize = -1; } //we will set it to smallest group later
212 else { subsample = false; }
215 if (subsample == false) { iters = 0; }
219 catch(exception& e) {
220 m->errorOut(e, "GetMetaCommunityCommand", "GetMetaCommunityCommand");
224 //**********************************************************************************************************************
226 int GetMetaCommunityCommand::execute(){
229 if (abort == true) { if (calledHelp) { return 0; } return 2; }
231 InputData input(sharedfile, "sharedfile");
232 vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
233 string lastLabel = lookup[0]->getLabel();
235 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
236 set<string> processedLabels;
237 set<string> userLabels = labels;
240 if (subsampleSize == -1) { //user has not set size, set size = smallest samples size
241 subsampleSize = lookup[0]->getNumSeqs();
242 for (int i = 1; i < lookup.size(); i++) {
243 int thisSize = lookup[i]->getNumSeqs();
245 if (thisSize < subsampleSize) { subsampleSize = thisSize; }
250 vector<SharedRAbundVector*> temp;
251 for (int i = 0; i < lookup.size(); i++) {
252 if (lookup[i]->getNumSeqs() < subsampleSize) {
253 m->mothurOut(lookup[i]->getGroup() + " contains " + toString(lookup[i]->getNumSeqs()) + ". Eliminating."); m->mothurOutEndLine();
256 Groups.push_back(lookup[i]->getGroup());
257 temp.push_back(lookup[i]);
261 m->setGroups(Groups);
264 if (lookup.size() < 2) { m->mothurOut("You have not provided enough valid groups. I cannot run the command."); m->mothurOutEndLine(); m->control_pressed = true; return 0; }
268 //as long as you are not at the end of the file or done wih the lines you want
269 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
271 if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
273 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
275 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
277 createProcesses(lookup);
279 processedLabels.insert(lookup[0]->getLabel());
280 userLabels.erase(lookup[0]->getLabel());
283 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
284 string saveLabel = lookup[0]->getLabel();
286 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
287 lookup = input.getSharedRAbundVectors(lastLabel);
288 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
290 createProcesses(lookup);
292 processedLabels.insert(lookup[0]->getLabel());
293 userLabels.erase(lookup[0]->getLabel());
295 //restore real lastlabel to save below
296 lookup[0]->setLabel(saveLabel);
299 lastLabel = lookup[0]->getLabel();
300 //prevent memory leak
301 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; }
303 if (m->control_pressed) { return 0; }
305 //get next line to process
306 lookup = input.getSharedRAbundVectors();
309 if (m->control_pressed) { return 0; }
311 //output error messages about any remaining user labels
312 set<string>::iterator it;
313 bool needToRun = false;
314 for (it = userLabels.begin(); it != userLabels.end(); it++) {
315 m->mothurOut("Your file does not include the label " + *it);
316 if (processedLabels.count(lastLabel) != 1) {
317 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
320 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
324 //run last label if you need to
325 if (needToRun == true) {
326 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
327 lookup = input.getSharedRAbundVectors(lastLabel);
329 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
331 createProcesses(lookup);
333 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
336 //output files created by command
337 m->mothurOutEndLine();
338 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
339 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
340 m->mothurOutEndLine();
344 catch(exception& e) {
345 m->errorOut(e, "GetMetaCommunityCommand", "execute");
349 //**********************************************************************************************************************
350 int GetMetaCommunityCommand::createProcesses(vector<SharedRAbundVector*>& thislookup){
353 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
355 processors=1; //qFinderDMM not thread safe
358 vector<int> processIDS;
361 int minPartition = 0;
364 if (maxpartitions < processors) { processors = maxpartitions; }
366 map<string, string> variables;
367 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
368 variables["[distance]"] = thislookup[0]->getLabel();
369 string outputFileName = getOutputFileName("fit", variables);
370 outputNames.push_back(outputFileName); outputTypes["fit"].push_back(outputFileName);
372 //divide the partitions between the processors
373 vector< vector<int> > dividedPartitions;
374 vector< vector<string> > rels, matrix;
375 vector<string> doneFlags;
376 dividedPartitions.resize(processors);
377 rels.resize(processors);
378 matrix.resize(processors);
380 //for each file group figure out which process will complete it
381 //want to divide the load intelligently so the big files are spread between processes
382 for (int i=1; i<=maxpartitions; i++) {
384 int processToAssign = (i+1) % processors;
385 if (processToAssign == 0) { processToAssign = processors; }
387 if (m->debug) { m->mothurOut("[DEBUG]: assigning " + toString(i) + " to process " + toString(processToAssign-1) + "\n"); }
388 dividedPartitions[(processToAssign-1)].push_back(i);
390 variables["[tag]"] = toString(i);
391 string relName = getOutputFileName("relabund", variables);
392 string mName = getOutputFileName("matrix", variables);
393 rels[(processToAssign-1)].push_back(relName);
394 matrix[(processToAssign-1)].push_back(mName);
397 for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
398 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
399 doneFlags.push_back(tempDoneFile);
401 m->openOutputFile(tempDoneFile, out); //clear out
406 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
408 //loop through and create all the processes you want
409 while (process != processors) {
413 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
417 num = processDriver(thislookup, dividedPartitions[process], (outputFileName + toString(getpid())), rels[process], matrix[process], doneFlags, process);
419 //pass numSeqs to parent
421 string tempFile = toString(getpid()) + ".outputNames.temp";
422 m->openOutputFile(tempFile, out);
424 out << outputNames.size() << endl;
425 for (int i = 0; i < outputNames.size(); i++) { out << outputNames[i] << endl; }
430 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
431 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
437 if (method == "dmm") { m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n"); }
439 m->mothurOut("K\tCH\t");
440 for (int i = 0; i < thislookup.size(); i++) { m->mothurOut(thislookup[i]->getGroup() + '\t'); }
443 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
445 //force parent to wait until all the processes are done
446 for (int i=0;i<processIDS.size();i++) {
447 int temp = processIDS[i];
451 vector<string> tempOutputNames = outputNames;
452 for (int i=0;i<processIDS.size();i++) {
454 string tempFile = toString(processIDS[i]) + ".outputNames.temp";
455 m->openInputFile(tempFile, in);
458 in >> tempNum; m->gobble(in);
459 if (tempNum < minPartition) { minPartition = tempNum; }
460 in >> tempNum; m->gobble(in);
461 for (int i = 0; i < tempNum; i++) {
462 string tempName = "";
463 in >> tempName; m->gobble(in);
464 tempOutputNames.push_back(tempName);
467 in.close(); m->mothurRemove(tempFile);
469 m->appendFilesWithoutHeaders(outputFileName + toString(processIDS[i]), outputFileName);
470 m->mothurRemove(outputFileName + toString(processIDS[i]));
473 if (processors > 1) {
475 for (int i = 0; i < tempOutputNames.size(); i++) { //remove files if needed
476 string name = tempOutputNames[i];
477 vector<string> parts;
478 m->splitAtChar(name, parts, '.');
480 if (((parts[parts.size()-1] == "relabund") || (parts[parts.size()-1] == "posterior")) && (parts[parts.size()-2] == "mix")) {
481 string tempNum = parts[parts.size()-3];
482 int num; m->mothurConvert(tempNum, num);
483 //if (num > minPartition) {
484 // m->mothurRemove(tempOutputNames[i]);
485 // keep = false; if (m->debug) { m->mothurOut("[DEBUG]: removing " + tempOutputNames[i] + ".\n"); }
488 if (keep) { outputNames.push_back(tempOutputNames[i]); }
493 m->openInputFile(outputFileName, in);
494 string headers = m->getline(in); m->gobble(in);
496 map<int, string> file;
498 string numString, line;
500 in >> numString; line = m->getline(in); m->gobble(in);
501 m->mothurConvert(numString, num);
506 m->openOutputFile(outputFileName, out);
507 out << headers << endl;
508 for (map<int, string>::iterator it = file.begin(); it != file.end(); it++) {
509 out << it->first << '\t' << it->second << endl;
510 if (m->debug) { m->mothurOut("[DEBUG]: printing: " + toString(it->first) + '\t' + it->second + ".\n"); }
516 m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");
517 minPartition = processDriver(thislookup, dividedPartitions[0], outputFileName, rels[0], matrix[0], doneFlags, 0);
519 for (int i = 0; i < processors; i++) { //read from everyone elses, just write to yours
520 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(i) + ".done.temp";
521 m->mothurRemove(tempDoneFile);
524 if (m->control_pressed) { return 0; }
526 if (m->debug) { m->mothurOut("[DEBUG]: minPartition = " + toString(minPartition) + "\n"); }
528 //run generate Summary function for smallest minPartition
529 variables["[tag]"] = toString(minPartition);
530 vector<double> piValues = generateDesignFile(minPartition, variables);
531 if (method == "dmm") { generateSummaryFile(minPartition, variables, piValues); } //pam doesn't make a relabund file
536 catch(exception& e) {
537 m->errorOut(e, "GetMetaCommunityCommand", "createProcesses");
541 //**********************************************************************************************************************
542 int GetMetaCommunityCommand::processDriver(vector<SharedRAbundVector*>& thislookup, vector<int>& parts, string outputFileName, vector<string> relabunds, vector<string> matrix, vector<string> doneFlags, int processID){
545 double minLaplace = 1e10;
546 int minPartition = 1;
547 vector<double> minSilhouettes; minSilhouettes.resize(thislookup.size(), 0);
549 ofstream fitData, silData;
550 if (method == "dmm") {
551 m->openOutputFile(outputFileName, fitData);
552 fitData.setf(ios::fixed, ios::floatfield);
553 fitData.setf(ios::showpoint);
554 fitData << "K\tNLE\tlogDet\tBIC\tAIC\tLaplace" << endl;
555 }else if((method == "pam") || (method == "kmeans")) { //because ch is looking of maximal value
557 m->openOutputFile(outputFileName, silData);
558 silData.setf(ios::fixed, ios::floatfield);
559 silData.setf(ios::showpoint);
560 silData << "K\tCH\t";
561 for (int i = 0; i < thislookup.size(); i++) { silData << thislookup[i]->getGroup() << '\t'; }
565 cout.setf(ios::fixed, ios::floatfield);
566 cout.setf(ios::showpoint);
568 vector< vector<int> > sharedMatrix;
569 vector<string> thisGroups;
570 for (int i = 0; i < thislookup.size(); i++) { sharedMatrix.push_back(thislookup[i]->getAbundances()); thisGroups.push_back(thislookup[i]->getGroup()); }
572 vector< vector<double> > dists; //do we want to output this matrix??
573 if ((method == "pam") || (method == "kmeans")) { dists = generateDistanceMatrix(thislookup); }
576 m->mothurOut("[DEBUG]: dists = \n");
577 for (int i = 0; i < dists.size(); i++) {
578 if (m->control_pressed) { break; }
579 m->mothurOut("[DEBUG]: i = " + toString(i) + '\t');
580 for (int j = 0; j < i; j++) { m->mothurOut(toString(dists[i][j]) +"\t"); }
585 for(int i=0;i<parts.size();i++){
587 int numPartitions = parts[i];
589 if (m->debug) { m->mothurOut("[DEBUG]: running partition " + toString(numPartitions) + "\n"); }
591 if (m->control_pressed) { break; }
593 //check to see if anyone else is done
594 for (int j = 0; j < doneFlags.size(); j++) {
595 if (!m->isBlank(doneFlags[j])) { //another process has finished
596 //are they done at a lower partition?
598 m->openInputFile(doneFlags[j], in);
600 in >> tempNum; in.close();
601 if (tempNum < numPartitions) { break; } //quit, because someone else has finished
605 CommunityTypeFinder* finder = NULL;
606 if (method == "dmm") { finder = new qFinderDMM(sharedMatrix, numPartitions); }
607 else if (method == "kmeans") { finder = new KMeans(sharedMatrix, numPartitions); }
608 else if (method == "pam") { finder = new Pam(sharedMatrix, dists, numPartitions); }
610 if (i == 0) { m->mothurOut(method + " is not a valid method option. I will run the command using dmm.\n"); }
611 finder = new qFinderDMM(sharedMatrix, numPartitions);
614 string relabund = relabunds[i];
615 string matrixName = matrix[i];
616 outputNames.push_back(matrixName); outputTypes["matrix"].push_back(matrixName);
618 finder->printZMatrix(matrixName, thisGroups);
620 double chi; vector<double> silhouettes;
621 if (method == "dmm") {
622 double laplace = finder->getLaplace();
623 if(laplace < minLaplace){
624 minPartition = numPartitions;
625 minLaplace = laplace;
628 chi = finder->calcCHIndex(dists);
629 silhouettes = finder->calcSilhouettes(dists);
630 if (chi > minLaplace) { //save partition with maximum ch index score
631 minPartition = numPartitions;
633 minSilhouettes = silhouettes;
637 if (method == "dmm") {
638 finder->printFitData(cout, minLaplace);
639 finder->printFitData(fitData);
640 finder->printRelAbund(relabund, m->currentSharedBinLabels);
641 outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
642 }else if ((method == "pam") || (method == "kmeans")) { //print silouettes and ch values
643 finder->printSilData(cout, chi, silhouettes);
644 finder->printSilData(silData, chi, silhouettes);
645 if (method == "kmeans") {
646 finder->printRelAbund(relabund, m->currentSharedBinLabels);
647 outputNames.push_back(relabund); outputTypes["relabund"].push_back(relabund);
652 if(optimizegap != -1 && (numPartitions - minPartition) >= optimizegap && numPartitions >= minpartitions){
653 string tempDoneFile = m->getRootName(m->getSimpleName(sharedfile)) + toString(processID) + ".done.temp";
655 m->openOutputFile(tempDoneFile, outDone);
656 outDone << minPartition << endl;
661 if (method == "dmm") { fitData.close(); }
663 if (m->control_pressed) { return 0; }
667 catch(exception& e) {
668 m->errorOut(e, "GetMetaCommunityCommand", "processDriver");
672 /**************************************************************************************************/
674 vector<double> GetMetaCommunityCommand::generateDesignFile(int numPartitions, map<string,string> variables){
676 vector<double> piValues(numPartitions, 0);
679 variables["[tag]"] = toString(numPartitions);
680 string input = getOutputFileName("matrix", variables);
681 m->openInputFile(input, postFile);//((fileRoot + toString(numPartitions) + "mix.posterior").c_str()); //matrix file
682 variables.erase("[tag]");
683 string outputFileName = getOutputFileName("design", variables);
685 m->openOutputFile(outputFileName, designFile);
686 outputNames.push_back(outputFileName); outputTypes["design"].push_back(outputFileName);
689 vector<string> titles(numPartitions);
691 for(int i=0;i<numPartitions;i++){ postFile >> titles[i]; }
699 if (m->control_pressed) { break; }
701 double maxPosterior = 0.0000;
702 int maxPartition = -1;
704 postFile >> sampleName;
706 for(int i=0;i<numPartitions;i++){
708 postFile >> posterior;
709 if(posterior > maxPosterior){
710 maxPosterior = posterior;
713 piValues[i] += posterior;
717 designFile << sampleName << '\t' << titles[maxPartition] << endl;
722 for(int i=0;i<numPartitions;i++){
723 piValues[i] /= (double)numSamples;
732 catch(exception& e) {
733 m->errorOut(e, "GetMetaCommunityCommand", "generateDesignFile");
738 /**************************************************************************************************/
740 inline bool summaryFunction(summaryData i, summaryData j){ return i.difference > j.difference; }
742 /**************************************************************************************************/
743 int GetMetaCommunityCommand::generateSummaryFile(int numPartitions, map<string,string> v, vector<double> piValues){
745 vector<summaryData> summary;
747 vector<double> pMean(numPartitions, 0);
748 vector<double> pLCI(numPartitions, 0);
749 vector<double> pUCI(numPartitions, 0);
752 double mean, lci, uci;
754 ifstream referenceFile;
755 map<string, string> variables;
756 variables["[filename]"] = v["[filename]"];
757 variables["[distance]"] = v["[distance]"];
758 variables["[tag]"] = "1";
759 string reference = getOutputFileName("relabund", variables);
760 m->openInputFile(reference, referenceFile); //((fileRoot + label + ".1mix.relabund").c_str());
761 variables["[tag]"] = toString(numPartitions);
762 string partFile = getOutputFileName("relabund", variables);
763 ifstream partitionFile;
764 m->openInputFile(partFile, partitionFile); //((fileRoot + toString(numPartitions) + "mix.relabund").c_str());
766 header = m->getline(referenceFile);
767 header = m->getline(partitionFile);
768 stringstream head(header);
771 vector<string> thetaValues(numPartitions, "");
772 for(int i=0;i<numPartitions;i++){
773 head >> label >> dummy >> dummy;
774 thetaValues[i] = label.substr(label.find_last_of('_')+1);
778 vector<double> partitionDiff(numPartitions, 0.0000);
780 while(referenceFile){
782 if (m->control_pressed) { break; }
784 referenceFile >> name >> mean >> lci >> uci;
786 summaryData tempData;
787 tempData.name = name;
788 tempData.refMean = mean;
790 double difference = 0.0000;
792 partitionFile >> name;
793 for(int j=0;j<numPartitions;j++){
794 partitionFile >> pMean[j] >> pLCI[j] >> pUCI[j];
795 difference += abs(mean - pMean[j]);
796 partitionDiff[j] += abs(mean - pMean[j]);;
799 tempData.partMean = pMean;
800 tempData.partLCI = pLCI;
801 tempData.partUCI = pUCI;
802 tempData.difference = difference;
803 summary.push_back(tempData);
805 m->gobble(referenceFile);
806 m->gobble(partitionFile);
808 referenceFile.close();
809 partitionFile.close();
811 if (m->control_pressed) { return 0; }
813 int numOTUs = (int)summary.size();
815 sort(summary.begin(), summary.end(), summaryFunction);
817 variables.erase("[tag]");
818 string outputFileName = getOutputFileName("parameters", variables);
819 outputNames.push_back(outputFileName); outputTypes["parameters"].push_back(outputFileName);
821 ofstream parameterFile;
822 m->openOutputFile(outputFileName, parameterFile); //((fileRoot + "mix.parameters").c_str());
823 parameterFile.setf(ios::fixed, ios::floatfield);
824 parameterFile.setf(ios::showpoint);
826 double totalDifference = 0.0000;
827 parameterFile << "Part\tDif2Ref_i\ttheta_i\tpi_i\n";
828 for(int i=0;i<numPartitions;i++){
829 if (m->control_pressed) { break; }
830 parameterFile << i+1 << '\t' << setprecision(2) << partitionDiff[i] << '\t' << thetaValues[i] << '\t' << piValues[i] << endl;
831 totalDifference += partitionDiff[i];
833 parameterFile.close();
835 if (m->control_pressed) { return 0; }
837 string summaryFileName = getOutputFileName("summary", variables);
838 outputNames.push_back(summaryFileName); outputTypes["summary"].push_back(summaryFileName);
840 ofstream summaryFile;
841 m->openOutputFile(summaryFileName, summaryFile); //((fileRoot + "mix.summary").c_str());
842 summaryFile.setf(ios::fixed, ios::floatfield);
843 summaryFile.setf(ios::showpoint);
846 summaryFile << "OTU\tP0.mean";
847 for(int i=0;i<numPartitions;i++){
848 summaryFile << "\tP" << i+1 << ".mean\tP" << i+1 << ".lci\tP" << i+1 << ".uci";
850 summaryFile << "\tDifference\tCumFraction" << endl;
852 double cumDiff = 0.0000;
854 for(int i=0;i<numOTUs;i++){
855 if (m->control_pressed) { break; }
856 summaryFile << summary[i].name << setprecision(2) << '\t' << summary[i].refMean;
857 for(int j=0;j<numPartitions;j++){
858 summaryFile << '\t' << summary[i].partMean[j] << '\t' << summary[i].partLCI[j] << '\t' << summary[i].partUCI[j];
861 cumDiff += summary[i].difference/totalDifference;
862 summaryFile << '\t' << summary[i].difference << '\t' << cumDiff << endl;
869 catch(exception& e) {
870 m->errorOut(e, "GetMetaCommunityCommand", "generateSummaryFile");
875 //**********************************************************************************************************************
876 vector<vector<double> > GetMetaCommunityCommand::generateDistanceMatrix(vector<SharedRAbundVector*>& thisLookup){
878 vector<vector<double> > results;
880 Calculator* matrixCalculator;
881 ValidCalculators validCalculator;
884 if (validCalculator.isValidCalculator("matrix", Estimators[i]) == true) {
885 if (Estimators[i] == "sharedsobs") {
886 matrixCalculator = new SharedSobsCS();
887 }else if (Estimators[i] == "sharedchao") {
888 matrixCalculator = new SharedChao1();
889 }else if (Estimators[i] == "sharedace") {
890 matrixCalculator = new SharedAce();
891 }else if (Estimators[i] == "jabund") {
892 matrixCalculator = new JAbund();
893 }else if (Estimators[i] == "sorabund") {
894 matrixCalculator = new SorAbund();
895 }else if (Estimators[i] == "jclass") {
896 matrixCalculator = new Jclass();
897 }else if (Estimators[i] == "sorclass") {
898 matrixCalculator = new SorClass();
899 }else if (Estimators[i] == "jest") {
900 matrixCalculator = new Jest();
901 }else if (Estimators[i] == "sorest") {
902 matrixCalculator = new SorEst();
903 }else if (Estimators[i] == "thetayc") {
904 matrixCalculator = new ThetaYC();
905 }else if (Estimators[i] == "thetan") {
906 matrixCalculator = new ThetaN();
907 }else if (Estimators[i] == "kstest") {
908 matrixCalculator = new KSTest();
909 }else if (Estimators[i] == "sharednseqs") {
910 matrixCalculator = new SharedNSeqs();
911 }else if (Estimators[i] == "ochiai") {
912 matrixCalculator = new Ochiai();
913 }else if (Estimators[i] == "anderberg") {
914 matrixCalculator = new Anderberg();
915 }else if (Estimators[i] == "kulczynski") {
916 matrixCalculator = new Kulczynski();
917 }else if (Estimators[i] == "kulczynskicody") {
918 matrixCalculator = new KulczynskiCody();
919 }else if (Estimators[i] == "lennon") {
920 matrixCalculator = new Lennon();
921 }else if (Estimators[i] == "morisitahorn") {
922 matrixCalculator = new MorHorn();
923 }else if (Estimators[i] == "braycurtis") {
924 matrixCalculator = new BrayCurtis();
925 }else if (Estimators[i] == "whittaker") {
926 matrixCalculator = new Whittaker();
927 }else if (Estimators[i] == "odum") {
928 matrixCalculator = new Odum();
929 }else if (Estimators[i] == "canberra") {
930 matrixCalculator = new Canberra();
931 }else if (Estimators[i] == "structeuclidean") {
932 matrixCalculator = new StructEuclidean();
933 }else if (Estimators[i] == "structchord") {
934 matrixCalculator = new StructChord();
935 }else if (Estimators[i] == "hellinger") {
936 matrixCalculator = new Hellinger();
937 }else if (Estimators[i] == "manhattan") {
938 matrixCalculator = new Manhattan();
939 }else if (Estimators[i] == "structpearson") {
940 matrixCalculator = new StructPearson();
941 }else if (Estimators[i] == "soergel") {
942 matrixCalculator = new Soergel();
943 }else if (Estimators[i] == "spearman") {
944 matrixCalculator = new Spearman();
945 }else if (Estimators[i] == "structkulczynski") {
946 matrixCalculator = new StructKulczynski();
947 }else if (Estimators[i] == "speciesprofile") {
948 matrixCalculator = new SpeciesProfile();
949 }else if (Estimators[i] == "hamming") {
950 matrixCalculator = new Hamming();
951 }else if (Estimators[i] == "structchi2") {
952 matrixCalculator = new StructChi2();
953 }else if (Estimators[i] == "gower") {
954 matrixCalculator = new Gower();
955 }else if (Estimators[i] == "memchi2") {
956 matrixCalculator = new MemChi2();
957 }else if (Estimators[i] == "memchord") {
958 matrixCalculator = new MemChord();
959 }else if (Estimators[i] == "memeuclidean") {
960 matrixCalculator = new MemEuclidean();
961 }else if (Estimators[i] == "mempearson") {
962 matrixCalculator = new MemPearson();
963 }else if (Estimators[i] == "jsd") {
964 matrixCalculator = new JSD();
965 }else if (Estimators[i] == "rjsd") {
966 matrixCalculator = new RJSD();
968 m->mothurOut("[ERROR]: " + Estimators[i] + " is not a valid calculator, please correct.\n"); m->control_pressed = true; return results;
973 vector< vector< vector<seqDist> > > calcDistsTotals; //each iter, then each groupCombos dists. this will be used to make .dist files
974 vector< vector<seqDist> > calcDists; calcDists.resize(1);
976 for (int thisIter = 0; thisIter < iters+1; thisIter++) {
978 vector<SharedRAbundVector*> thisItersLookup = thisLookup;
980 if (subsample && (thisIter != 0)) {
982 vector<string> tempLabels; //dont need since we arent printing the sampled sharedRabunds
984 //make copy of lookup so we don't get access violations
985 vector<SharedRAbundVector*> newLookup;
986 for (int k = 0; k < thisItersLookup.size(); k++) {
987 SharedRAbundVector* temp = new SharedRAbundVector();
988 temp->setLabel(thisItersLookup[k]->getLabel());
989 temp->setGroup(thisItersLookup[k]->getGroup());
990 newLookup.push_back(temp);
994 for (int k = 0; k < thisItersLookup[0]->getNumBins(); k++) {
995 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return results; }
996 for (int j = 0; j < thisItersLookup.size(); j++) { newLookup[j]->push_back(thisItersLookup[j]->getAbundance(k), thisItersLookup[j]->getGroup()); }
999 tempLabels = sample.getSample(newLookup, subsampleSize);
1000 thisItersLookup = newLookup;
1004 driver(thisItersLookup, calcDists, matrixCalculator);
1006 if (subsample && (thisIter != 0)) {
1007 if((thisIter) % 100 == 0){ m->mothurOutJustToScreen(toString(thisIter)+"\n"); }
1008 calcDistsTotals.push_back(calcDists);
1009 for (int i = 0; i < calcDists.size(); i++) {
1010 for (int j = 0; j < calcDists[i].size(); j++) {
1011 if (m->debug) { m->mothurOut("[DEBUG]: Results: iter = " + toString(thisIter) + ", " + thisLookup[calcDists[i][j].seq1]->getGroup() + " - " + thisLookup[calcDists[i][j].seq2]->getGroup() + " distance = " + toString(calcDists[i][j].dist) + ".\n"); }
1015 for (int i = 0; i < thisItersLookup.size(); i++) { delete thisItersLookup[i]; }
1016 thisItersLookup.clear();
1017 }else { //print results for whole dataset
1018 for (int i = 0; i < calcDists.size(); i++) {
1019 if (m->control_pressed) { break; }
1022 results.resize(thisLookup.size());
1023 for (int k = 0; k < thisLookup.size(); k++) { results[k].resize(thisLookup.size(), 0.0); }
1025 for (int j = 0; j < calcDists[i].size(); j++) {
1026 int row = calcDists[i][j].seq1;
1027 int column = calcDists[i][j].seq2;
1028 double dist = calcDists[i][j].dist;
1030 results[row][column] = dist;
1031 results[column][row] = dist;
1035 for (int i = 0; i < calcDists.size(); i++) { calcDists[i].clear(); }
1039 //we need to find the average distance and standard deviation for each groups distance
1040 vector< vector<seqDist> > calcAverages = m->getAverages(calcDistsTotals, "average");
1043 for (int i = 0; i < calcDists.size(); i++) {
1044 results.resize(thisLookup.size());
1045 for (int k = 0; k < thisLookup.size(); k++) { results[k].resize(thisLookup.size(), 0.0); }
1047 for (int j = 0; j < calcAverages[i].size(); j++) {
1048 int row = calcAverages[i][j].seq1;
1049 int column = calcAverages[i][j].seq2;
1050 float dist = calcAverages[i][j].dist;
1052 results[row][column] = dist;
1053 results[column][row] = dist;
1061 catch(exception& e) {
1062 m->errorOut(e, "GetMetaCommunityCommand", "generateDistanceMatrix");
1066 /**************************************************************************************************/
1067 int GetMetaCommunityCommand::driver(vector<SharedRAbundVector*> thisLookup, vector< vector<seqDist> >& calcDists, Calculator* matrixCalculator) {
1069 vector<SharedRAbundVector*> subset;
1071 for (int k = 0; k < thisLookup.size(); k++) { // pass cdd each set of groups to compare
1073 for (int l = 0; l < k; l++) {
1075 if (k != l) { //we dont need to similiarity of a groups to itself
1076 subset.clear(); //clear out old pair of sharedrabunds
1077 //add new pair of sharedrabunds
1078 subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
1082 //if this calc needs all groups to calculate the pair load all groups
1083 if (matrixCalculator->getNeedsAll()) {
1084 //load subset with rest of lookup for those calcs that need everyone to calc for a pair
1085 for (int w = 0; w < thisLookup.size(); w++) {
1086 if ((w != k) && (w != l)) { subset.push_back(thisLookup[w]); }
1090 vector<double> tempdata = matrixCalculator->getValues(subset); //saves the calculator outputs
1092 if (m->control_pressed) { return 1; }
1094 seqDist temp(l, k, tempdata[0]);
1095 //cout << l << '\t' << k << '\t' << tempdata[0] << endl;
1096 calcDists[0].push_back(temp);
1104 catch(exception& e) {
1105 m->errorOut(e, "MatrixOutputCommand", "driver");
1109 //**********************************************************************************************************************