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