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