unifracunweightedcommand.cpp

   1 /*
   2  *  unifracunweightedcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 2/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "unifracunweightedcommand.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> UnifracUnweightedCommand::setParameters(){
  14         try {
  15                 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
  16                 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
  17                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
  18                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  19                 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
  20                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  21                 CommandParameter prandom("random", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prandom);
  22                 CommandParameter pdistance("distance", "Multiple", "column-lt-square", "column", "", "", "",false,false); parameters.push_back(pdistance);
  23                 CommandParameter proot("root", "Boolean", "F", "", "", "", "",false,false); parameters.push_back(proot);
  24                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  25                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  26
  27                 vector<string> myArray;
  28                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  29                 return myArray;
  30         }
  31         catch(exception& e) {
  32                 m->errorOut(e, "UnifracUnweightedCommand", "setParameters");
  33                 exit(1);
  34         }
  35 }
  36 //**********************************************************************************************************************
  37 string UnifracUnweightedCommand::getHelpString(){
  38         try {
  39                 string helpString = "";
  40                 helpString += "The unifrac.unweighted command parameters are tree, group, name, groups, iters, distance, processors, root and random.  tree parameter is required unless you have valid current tree file.\n";
  41                 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed.  You must enter at least 1 valid group.\n";
  42                 helpString += "The group names are separated by dashes.  The iters parameter allows you to specify how many random trees you would like compared to your tree.\n";
  43                 helpString += "The distance parameter allows you to create a distance file from the results. The default is false. You may set distance to lt, square or column.\n";
  44                 helpString += "The random parameter allows you to shut off the comparison to random trees. The default is false, meaning compare don't your trees with randomly generated trees.\n";
  45                 helpString += "The root parameter allows you to include the entire root in your calculations. The default is false, meaning stop at the root for this comparision instead of the root of the entire tree.\n";
  46                 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
  47                 helpString += "The unifrac.unweighted command should be in the following format: unifrac.unweighted(groups=yourGroups, iters=yourIters).\n";
  48                 helpString += "Example unifrac.unweighted(groups=A-B-C, iters=500).\n";
  49                 helpString += "The default value for groups is all the groups in your groupfile, and iters is 1000.\n";
  50                 helpString += "The unifrac.unweighted command output two files: .unweighted and .uwsummary their descriptions are in the manual.\n";
  51                 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
  52                 return helpString;
  53         }
  54         catch(exception& e) {
  55                 m->errorOut(e, "UnifracUnweightedCommand", "getHelpString");
  56                 exit(1);
  57         }
  58 }
  59 //**********************************************************************************************************************
  60 UnifracUnweightedCommand::UnifracUnweightedCommand(){
  61         try {
  62                 abort = true; calledHelp = true;
  63                 setParameters();
  64                 vector<string> tempOutNames;
  65                 outputTypes["unweighted"] = tempOutNames;
  66                 outputTypes["uwsummary"] = tempOutNames;
  67                 outputTypes["phylip"] = tempOutNames;
  68                 outputTypes["column"] = tempOutNames;
  69         }
  70         catch(exception& e) {
  71                 m->errorOut(e, "UnifracUnweightedCommand", "UnifracUnweightedCommand");
  72                 exit(1);
  73         }
  74 }
  75 /***********************************************************/
  76 UnifracUnweightedCommand::UnifracUnweightedCommand(string option)  {
  77         try {
  78                 abort = false; calledHelp = false;
  79
  80
  81                 //allow user to run help
  82                 if(option == "help") { help(); abort = true; calledHelp = true; }
  83                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  84
  85                 else {
  86                         vector<string> myArray = setParameters();
  87
  88                         OptionParser parser(option);
  89                         map<string,string> parameters = parser.getParameters();
  90                         map<string,string>::iterator it;
  91
  92                         ValidParameters validParameter;
  93
  94                         //check to make sure all parameters are valid for command
  95                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  96                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  97                         }
  98
  99                         //initialize outputTypes
 100                         vector<string> tempOutNames;
 101                         outputTypes["unweighted"] = tempOutNames;
 102                         outputTypes["uwsummary"] = tempOutNames;
 103                         outputTypes["phylip"] = tempOutNames;
 104                         outputTypes["column"] = tempOutNames;
 105
 106                         //if the user changes the input directory command factory will send this info to us in the output parameter
 107                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 108                         if (inputDir == "not found"){   inputDir = "";          }
 109                         else {
 110                                 string path;
 111                                 it = parameters.find("tree");
 112                                 //user has given a template file
 113                                 if(it != parameters.end()){
 114                                         path = m->hasPath(it->second);
 115                                         //if the user has not given a path then, add inputdir. else leave path alone.
 116                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
 117                                 }
 118
 119                                 it = parameters.find("group");
 120                                 //user has given a template file
 121                                 if(it != parameters.end()){
 122                                         path = m->hasPath(it->second);
 123                                         //if the user has not given a path then, add inputdir. else leave path alone.
 124                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 125                                 }
 126
 127                                 it = parameters.find("name");
 128                                 //user has given a template file
 129                                 if(it != parameters.end()){
 130                                         path = m->hasPath(it->second);
 131                                         //if the user has not given a path then, add inputdir. else leave path alone.
 132                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 133                                 }
 134                         }
 135
 136                         m->runParse = true;
 137                         m->Groups.clear();
 138                         m->namesOfGroups.clear();
 139                         m->Treenames.clear();
 140                         m->names.clear();
 141
 142                         //check for required parameters
 143                         treefile = validParameter.validFile(parameters, "tree", true);
 144                         if (treefile == "not open") { abort = true; }
 145                         else if (treefile == "not found") {                             //if there is a current design file, use it
 146                                 treefile = m->getTreeFile();
 147                                 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
 148                                 else {  m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
 149                         }else { m->setTreeFile(treefile); }
 150
 151                         //check for required parameters
 152                         groupfile = validParameter.validFile(parameters, "group", true);
 153                         if (groupfile == "not open") { abort = true; }
 154                         else if (groupfile == "not found") { groupfile = ""; }
 155                         else { m->setGroupFile(groupfile); }
 156
 157                         namefile = validParameter.validFile(parameters, "name", true);
 158                         if (namefile == "not open") { abort = true; }
 159                         else if (namefile == "not found") { namefile = ""; }
 160                         else { m->setNameFile(namefile); }
 161
 162                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 163
 164                         //check for optional parameter and set defaults
 165                         // ...at some point should added some additional type checking...
 166                         groups = validParameter.validFile(parameters, "groups", false);
 167                         if (groups == "not found") { groups = ""; }
 168                         else {
 169                                 m->splitAtDash(groups, Groups);
 170                                 m->Groups = Groups;
 171                         }
 172
 173                         itersString = validParameter.validFile(parameters, "iters", false);                             if (itersString == "not found") { itersString = "1000"; }
 174                         convert(itersString, iters);
 175
 176                         string temp = validParameter.validFile(parameters, "distance", false);
 177                         if (temp == "not found") { phylip = false; outputForm = ""; }
 178                         else{
 179                                 if ((temp == "lt") || (temp == "column") || (temp == "square")) {  phylip = true;  outputForm = temp; }
 180                                 else { m->mothurOut("Options for distance are: lt, square, or column. Using lt."); m->mothurOutEndLine(); phylip = true; outputForm = "lt"; }
 181                         }
 182
 183                         temp = validParameter.validFile(parameters, "random", false);                                   if (temp == "not found") { temp = "f"; }
 184                         random = m->isTrue(temp);
 185
 186                         temp = validParameter.validFile(parameters, "root", false);                                     if (temp == "not found") { temp = "F"; }
 187                         includeRoot = m->isTrue(temp);
 188
 189                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 190                         m->setProcessors(temp);
 191                         convert(temp, processors);
 192
 193                         if (!random) {  iters = 0;  } //turn off random calcs
 194
 195                         //if user selects distance = true and no groups it won't calc the pairwise
 196                         if ((phylip) && (Groups.size() == 0)) {
 197                                 groups = "all";
 198                                 m->splitAtDash(groups, Groups);
 199                                 m->Groups = Groups;
 200                         }
 201                 }
 202
 203         }
 204         catch(exception& e) {
 205                 m->errorOut(e, "UnifracUnweightedCommand", "UnifracUnweightedCommand");
 206                 exit(1);
 207         }
 208 }
 209
 210 /***********************************************************/
 211 int UnifracUnweightedCommand::execute() {
 212         try {
 213
 214                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 215
 216                 m->setTreeFile(treefile);
 217
 218                 if (groupfile != "") {
 219                         //read in group map info.
 220                         tmap = new TreeMap(groupfile);
 221                         tmap->readMap();
 222                 }else{ //fake out by putting everyone in one group
 223                         Tree* tree = new Tree(treefile); delete tree;  //extracts names from tree to make faked out groupmap
 224                         tmap = new TreeMap();
 225
 226                         for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
 227                 }
 228
 229                 if (namefile != "") { readNamesFile(); }
 230
 231                 read = new ReadNewickTree(treefile);
 232                 int readOk = read->read(tmap);
 233
 234                 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
 235
 236                 read->AssembleTrees();
 237                 T = read->getTrees();
 238                 delete read;
 239
 240                 //make sure all files match
 241                 //if you provide a namefile we will use the numNames in the namefile as long as the number of unique match the tree names size.
 242                 int numNamesInTree;
 243                 if (namefile != "")  {
 244                         if (numUniquesInName == m->Treenames.size()) {  numNamesInTree = nameMap.size();  }
 245                         else {   numNamesInTree = m->Treenames.size();  }
 246                 }else {  numNamesInTree = m->Treenames.size();  }
 247
 248
 249                 //output any names that are in group file but not in tree
 250                 if (numNamesInTree < tmap->getNumSeqs()) {
 251                         for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
 252                                 //is that name in the tree?
 253                                 int count = 0;
 254                                 for (int j = 0; j < m->Treenames.size(); j++) {
 255                                         if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
 256                                         count++;
 257                                 }
 258
 259                                 if (m->control_pressed) {
 260                                         delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 261                                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } outputTypes.clear();
 262                                         m->Groups.clear();
 263                                         return 0;
 264                                 }
 265
 266                                 //then you did not find it so report it
 267                                 if (count == m->Treenames.size()) {
 268                                         //if it is in your namefile then don't remove
 269                                         map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
 270
 271                                         if (it == nameMap.end()) {
 272                                                 m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
 273                                                 tmap->removeSeq(tmap->namesOfSeqs[i]);
 274                                                 i--; //need this because removeSeq removes name from namesOfSeqs
 275                                         }
 276                                 }
 277                         }
 278                 }
 279
 280                 sumFile = outputDir + m->getSimpleName(treefile) + ".uwsummary";
 281                 outputNames.push_back(sumFile); outputTypes["uwsummary"].push_back(sumFile);
 282                 m->openOutputFile(sumFile, outSum);
 283
 284                 util = new SharedUtil();
 285                 util->setGroups(m->Groups, tmap->namesOfGroups, allGroups, numGroups, "unweighted");    //sets the groups the user wants to analyze
 286                 util->getCombos(groupComb, m->Groups, numComp);
 287                 delete util;
 288
 289                 if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
 290
 291                 unweighted = new Unweighted(tmap, includeRoot);
 292
 293                 int start = time(NULL);
 294
 295                 userData.resize(numComp,0);  //data[0] = unweightedscore
 296                 randomData.resize(numComp,0); //data[0] = unweightedscore
 297                 //create new tree with same num nodes and leaves as users
 298
 299                 if (numComp < processors) { processors = numComp;  }
 300
 301                 outSum << "Tree#" << '\t' << "Groups" << '\t'  <<  "UWScore" <<'\t' << "UWSig" <<  endl;
 302                 m->mothurOut("Tree#\tGroups\tUWScore\tUWSig"); m->mothurOutEndLine();
 303
 304                 //get pscores for users trees
 305                 for (int i = 0; i < T.size(); i++) {
 306                         if (m->control_pressed) {
 307                                 delete tmap; delete unweighted;
 308                                 for (int i = 0; i < T.size(); i++) { delete T[i]; }
 309                                 outSum.close();
 310                                 for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  }
 311                                 return 0;
 312                         }
 313
 314                         counter = 0;
 315
 316                         if (random)  {
 317                                 output = new ColumnFile(outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted", itersString);
 318                                 outputNames.push_back(outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted");
 319                                 outputTypes["unweighted"].push_back(outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted");
 320                         }
 321
 322
 323                         //get unweighted for users tree
 324                         rscoreFreq.resize(numComp);
 325                         rCumul.resize(numComp);
 326                         utreeScores.resize(numComp);
 327                         UWScoreSig.resize(numComp);
 328
 329                         userData = unweighted->getValues(T[i], processors, outputDir);  //userData[0] = unweightedscore
 330
 331                         if (m->control_pressed) { delete tmap; delete unweighted;
 332                                 for (int i = 0; i < T.size(); i++) { delete T[i]; }if (random) { delete output;  } outSum.close();  for (int i = 0; i < outputNames.size(); i++) {      remove(outputNames[i].c_str());  }return 0; }
 333
 334                         //output scores for each combination
 335                         for(int k = 0; k < numComp; k++) {
 336                                 //saves users score
 337                                 utreeScores[k].push_back(userData[k]);
 338
 339                                 //add users score to validscores
 340                                 validScores[userData[k]] = userData[k];
 341                         }
 342
 343                         //get unweighted scores for random trees - if random is false iters = 0
 344                         for (int j = 0; j < iters; j++) {
 345
 346                                 //we need a different getValues because when we swap the labels we only want to swap those in each pairwise comparison
 347                                 randomData = unweighted->getValues(T[i], "", "", processors, outputDir);
 348
 349                                 if (m->control_pressed) { delete tmap; delete unweighted;
 350                                         for (int i = 0; i < T.size(); i++) { delete T[i]; }if (random) { delete output;  } outSum.close(); for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str());  } return 0; }
 351
 352                                 for(int k = 0; k < numComp; k++) {
 353                                         //add trees unweighted score to map of scores
 354                                         map<float,float>::iterator it = rscoreFreq[k].find(randomData[k]);
 355                                         if (it != rscoreFreq[k].end()) {//already have that score
 356                                                 rscoreFreq[k][randomData[k]]++;
 357                                         }else{//first time we have seen this score
 358                                                 rscoreFreq[k][randomData[k]] = 1;
 359                                         }
 360
 361                                         //add randoms score to validscores
 362                                         validScores[randomData[k]] = randomData[k];
 363                                 }
 364
 365                                 //report progress
 366 //                              m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();
 367                         }
 368
 369                         for(int a = 0; a < numComp; a++) {
 370                                 float rcumul = 1.0000;
 371
 372                                 if (random) {
 373                                         //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
 374                                         for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 375                                                 //make rscoreFreq map and rCumul
 376                                                 map<float,float>::iterator it2 = rscoreFreq[a].find(it->first);
 377                                                 rCumul[a][it->first] = rcumul;
 378                                                 //get percentage of random trees with that info
 379                                                 if (it2 != rscoreFreq[a].end()) {  rscoreFreq[a][it->first] /= iters; rcumul-= it2->second;  }
 380                                                 else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
 381                                         }
 382                                         UWScoreSig[a].push_back(rCumul[a][userData[a]]);
 383                                 }else           {       UWScoreSig[a].push_back(0.0);                                           }
 384
 385                         }
 386
 387                         if (m->control_pressed) { delete tmap; delete unweighted;
 388                                 for (int i = 0; i < T.size(); i++) { delete T[i]; }if (random) { delete output;  } outSum.close(); for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str());  } return 0;  }
 389
 390                         //print output files
 391                         printUWSummaryFile(i);
 392                         if (random)  {  printUnweightedFile();  delete output;  }
 393                         if (phylip) {   createPhylipFile(i);            }
 394
 395                         rscoreFreq.clear();
 396                         rCumul.clear();
 397                         validScores.clear();
 398                         utreeScores.clear();
 399                         UWScoreSig.clear();
 400                 }
 401
 402
 403                 outSum.close();
 404                 m->Groups.clear();
 405                 delete tmap; delete unweighted;
 406                 for (int i = 0; i < T.size(); i++) { delete T[i]; }
 407
 408                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  }      return 0; }
 409
 410                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.unweighted."); m->mothurOutEndLine();
 411
 412                 //set phylip file as new current phylipfile
 413                 string current = "";
 414                 itTypes = outputTypes.find("phylip");
 415                 if (itTypes != outputTypes.end()) {
 416                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
 417                 }
 418
 419                 //set column file as new current columnfile
 420                 itTypes = outputTypes.find("column");
 421                 if (itTypes != outputTypes.end()) {
 422                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
 423                 }
 424
 425                 m->mothurOutEndLine();
 426                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 427                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 428                 m->mothurOutEndLine();
 429
 430                 return 0;
 431
 432         }
 433         catch(exception& e) {
 434                 m->errorOut(e, "UnifracUnweightedCommand", "execute");
 435                 exit(1);
 436         }
 437 }
 438 /***********************************************************/
 439 void UnifracUnweightedCommand::printUnweightedFile() {
 440         try {
 441                 vector<double> data;
 442                 vector<string> tags;
 443
 444                 tags.push_back("Score");
 445                 tags.push_back("RandFreq"); tags.push_back("RandCumul");
 446
 447                 for(int a = 0; a < numComp; a++) {
 448                         output->initFile(groupComb[a], tags);
 449                         //print each line
 450                         for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 451                                 data.push_back(it->first);  data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 452                                 output->output(data);
 453                                 data.clear();
 454                         }
 455                         output->resetFile();
 456                 }
 457         }
 458         catch(exception& e) {
 459                 m->errorOut(e, "UnifracUnweightedCommand", "printUnweightedFile");
 460                 exit(1);
 461         }
 462 }
 463
 464 /***********************************************************/
 465 void UnifracUnweightedCommand::printUWSummaryFile(int i) {
 466         try {
 467
 468                 //format output
 469                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 470
 471                 //print each line
 472
 473                 for(int a = 0; a < numComp; a++) {
 474                         outSum << i+1 << '\t';
 475                         m->mothurOut(toString(i+1) + "\t");
 476
 477                         if (random) {
 478                                 if (UWScoreSig[a][0] > (1/(float)iters)) {
 479                                         outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
 480                                         cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
 481                                         m->mothurOutJustToLog(groupComb[a]  + "\t" + toString(utreeScores[a][0])  + "\t" + toString(UWScoreSig[a][0])+ "\n");
 482                                 }else {
 483                                         outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 484                                         cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 485                                         m->mothurOutJustToLog(groupComb[a]  + "\t" + toString(utreeScores[a][0])  + "\t<" + toString((1/float(iters))) + "\n");
 486                                 }
 487                         }else{
 488                                 outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
 489                                 cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
 490                                 m->mothurOutJustToLog(groupComb[a]  + "\t" + toString(utreeScores[a][0])  + "\t0.00\n");
 491                         }
 492                 }
 493
 494         }
 495         catch(exception& e) {
 496                 m->errorOut(e, "UnifracUnweightedCommand", "printUWSummaryFile");
 497                 exit(1);
 498         }
 499 }
 500 /***********************************************************/
 501 void UnifracUnweightedCommand::createPhylipFile(int i) {
 502         try {
 503                 string phylipFileName;
 504                 if ((outputForm == "lt") || (outputForm == "square")) {
 505                         phylipFileName = outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted.phylip.dist";
 506                         outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
 507                 }else { //column
 508                         phylipFileName = outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted.column.dist";
 509                         outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
 510                 }
 511
 512                 ofstream out;
 513                 m->openOutputFile(phylipFileName, out);
 514
 515                 if ((outputForm == "lt") || (outputForm == "square")) {
 516                         //output numSeqs
 517                         out << m->Groups.size() << endl;
 518                 }
 519
 520                 //make matrix with scores in it
 521                 vector< vector<float> > dists;  dists.resize(m->Groups.size());
 522                 for (int i = 0; i < m->Groups.size(); i++) {
 523                         dists[i].resize(m->Groups.size(), 0.0);
 524                 }
 525
 526                 //flip it so you can print it
 527                 int count = 0;
 528                 for (int r=0; r<m->Groups.size(); r++) {
 529                         for (int l = 0; l < r; l++) {
 530                                 dists[r][l] = utreeScores[count][0];
 531                                 dists[l][r] = utreeScores[count][0];
 532                                 count++;
 533                         }
 534                 }
 535
 536                 //output to file
 537                 for (int r=0; r<m->Groups.size(); r++) {
 538                         //output name
 539                         string name = m->Groups[r];
 540                         if (name.length() < 10) { //pad with spaces to make compatible
 541                                 while (name.length() < 10) {  name += " ";  }
 542                         }
 543
 544                         if (outputForm == "lt") {
 545                                 out << name << '\t';
 546
 547                                 //output distances
 548                                 for (int l = 0; l < r; l++) {   out  << dists[r][l] << '\t';  }
 549                                 out << endl;
 550                         }else if (outputForm == "square") {
 551                                 out << name << '\t';
 552
 553                                 //output distances
 554                                 for (int l = 0; l < m->Groups.size(); l++) {    out << dists[r][l] << '\t';  }
 555                                 out << endl;
 556                         }else{
 557                                 //output distances
 558                                 for (int l = 0; l < r; l++) {
 559                                         string otherName = m->Groups[l];
 560                                         if (otherName.length() < 10) { //pad with spaces to make compatible
 561                                                 while (otherName.length() < 10) {  otherName += " ";  }
 562                                         }
 563
 564                                         out  << name << '\t' << otherName << dists[r][l] << endl;
 565                                 }
 566                         }
 567                 }
 568                 out.close();
 569         }
 570         catch(exception& e) {
 571                 m->errorOut(e, "UnifracUnweightedCommand", "createPhylipFile");
 572                 exit(1);
 573         }
 574 }/*****************************************************************/
 575 int UnifracUnweightedCommand::readNamesFile() {
 576         try {
 577                 m->names.clear();
 578                 numUniquesInName = 0;
 579
 580                 ifstream in;
 581                 m->openInputFile(namefile, in);
 582
 583                 string first, second;
 584                 map<string, string>::iterator itNames;
 585
 586                 while(!in.eof()) {
 587                         in >> first >> second; m->gobble(in);
 588
 589                         numUniquesInName++;
 590
 591                         itNames = m->names.find(first);
 592                         if (itNames == m->names.end()) {
 593                                 m->names[first] = second;
 594
 595                                 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
 596                                 vector<string> dupNames;
 597                                 m->splitAtComma(second, dupNames);
 598
 599                                 for (int i = 0; i < dupNames.size(); i++) {
 600                                         nameMap[dupNames[i]] = dupNames[i];
 601                                         if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
 602                                 }
 603                         }else {  m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
 604                 }
 605                 in.close();
 606
 607                 return 0;
 608         }
 609         catch(exception& e) {
 610                 m->errorOut(e, "UnifracUnweightedCommand", "readNamesFile");
 611                 exit(1);
 612         }
 613 }
 614 /***********************************************************/
 615
 616
 617
 618