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->clearGroups();
 138                         m->clearAllGroups();
 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->setGroups(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->setGroups(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++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 262                                         m->clearGroups();
 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                 Groups = m->getGroups();
 286                 vector<string> namesGroups = tmap->getNamesOfGroups();
 287                 util->setGroups(Groups, namesGroups, allGroups, numGroups, "unweighted");       //sets the groups the user wants to analyze
 288                 util->getCombos(groupComb, Groups, numComp);
 289                 m->setGroups(Groups);
 290                 delete util;
 291
 292                 if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
 293
 294                 unweighted = new Unweighted(tmap, includeRoot);
 295
 296                 int start = time(NULL);
 297
 298                 userData.resize(numComp,0);  //data[0] = unweightedscore
 299                 randomData.resize(numComp,0); //data[0] = unweightedscore
 300                 //create new tree with same num nodes and leaves as users
 301
 302                 if (numComp < processors) { processors = numComp;  }
 303
 304                 outSum << "Tree#" << '\t' << "Groups" << '\t'  <<  "UWScore" <<'\t';
 305                 m->mothurOut("Tree#\tGroups\tUWScore\t");
 306                 if (random) { outSum << "UWSig"; m->mothurOut("UWSig"); }
 307                 outSum << endl; m->mothurOutEndLine();
 308
 309                 //get pscores for users trees
 310                 for (int i = 0; i < T.size(); i++) {
 311                         if (m->control_pressed) {
 312                                 delete tmap; delete unweighted;
 313                                 for (int i = 0; i < T.size(); i++) { delete T[i]; }
 314                                 outSum.close();
 315                                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  }
 316                                 return 0;
 317                         }
 318
 319                         counter = 0;
 320
 321                         if (random)  {
 322                                 output = new ColumnFile(outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted", itersString);
 323                                 outputNames.push_back(outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted");
 324                                 outputTypes["unweighted"].push_back(outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted");
 325                         }
 326
 327
 328                         //get unweighted for users tree
 329                         rscoreFreq.resize(numComp);
 330                         rCumul.resize(numComp);
 331                         utreeScores.resize(numComp);
 332                         UWScoreSig.resize(numComp);
 333
 334                         userData = unweighted->getValues(T[i], processors, outputDir);  //userData[0] = unweightedscore
 335
 336                         if (m->control_pressed) { delete tmap; delete unweighted;
 337                                 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++) {      m->mothurRemove(outputNames[i]);  }return 0; }
 338
 339                         //output scores for each combination
 340                         for(int k = 0; k < numComp; k++) {
 341                                 //saves users score
 342                                 utreeScores[k].push_back(userData[k]);
 343
 344                                 //add users score to validscores
 345                                 validScores[userData[k]] = userData[k];
 346                         }
 347
 348                         //get unweighted scores for random trees - if random is false iters = 0
 349                         for (int j = 0; j < iters; j++) {
 350
 351                                 //we need a different getValues because when we swap the labels we only want to swap those in each pairwise comparison
 352                                 randomData = unweighted->getValues(T[i], "", "", processors, outputDir);
 353
 354                                 if (m->control_pressed) { delete tmap; delete unweighted;
 355                                         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++) {       m->mothurRemove(outputNames[i]);  } return 0; }
 356
 357                                 for(int k = 0; k < numComp; k++) {
 358                                         //add trees unweighted score to map of scores
 359                                         map<float,float>::iterator it = rscoreFreq[k].find(randomData[k]);
 360                                         if (it != rscoreFreq[k].end()) {//already have that score
 361                                                 rscoreFreq[k][randomData[k]]++;
 362                                         }else{//first time we have seen this score
 363                                                 rscoreFreq[k][randomData[k]] = 1;
 364                                         }
 365
 366                                         //add randoms score to validscores
 367                                         validScores[randomData[k]] = randomData[k];
 368                                 }
 369
 370                                 //report progress
 371 //                              m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();
 372                         }
 373
 374                         for(int a = 0; a < numComp; a++) {
 375                                 float rcumul = 1.0000;
 376
 377                                 if (random) {
 378                                         //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
 379                                         for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 380                                                 //make rscoreFreq map and rCumul
 381                                                 map<float,float>::iterator it2 = rscoreFreq[a].find(it->first);
 382                                                 rCumul[a][it->first] = rcumul;
 383                                                 //get percentage of random trees with that info
 384                                                 if (it2 != rscoreFreq[a].end()) {  rscoreFreq[a][it->first] /= iters; rcumul-= it2->second;  }
 385                                                 else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
 386                                         }
 387                                         UWScoreSig[a].push_back(rCumul[a][userData[a]]);
 388                                 }else           {       UWScoreSig[a].push_back(0.0);                                           }
 389
 390                         }
 391
 392                         if (m->control_pressed) { delete tmap; delete unweighted;
 393                                 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++) {       m->mothurRemove(outputNames[i]);  } return 0;  }
 394
 395                         //print output files
 396                         printUWSummaryFile(i);
 397                         if (random)  {  printUnweightedFile();  delete output;  }
 398                         if (phylip) {   createPhylipFile(i);            }
 399
 400                         rscoreFreq.clear();
 401                         rCumul.clear();
 402                         validScores.clear();
 403                         utreeScores.clear();
 404                         UWScoreSig.clear();
 405                 }
 406
 407
 408                 outSum.close();
 409                 m->clearGroups();
 410                 delete tmap; delete unweighted;
 411                 for (int i = 0; i < T.size(); i++) { delete T[i]; }
 412
 413                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]);  }     return 0; }
 414
 415                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.unweighted."); m->mothurOutEndLine();
 416
 417                 //set phylip file as new current phylipfile
 418                 string current = "";
 419                 itTypes = outputTypes.find("phylip");
 420                 if (itTypes != outputTypes.end()) {
 421                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
 422                 }
 423
 424                 //set column file as new current columnfile
 425                 itTypes = outputTypes.find("column");
 426                 if (itTypes != outputTypes.end()) {
 427                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
 428                 }
 429
 430                 m->mothurOutEndLine();
 431                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 432                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 433                 m->mothurOutEndLine();
 434
 435                 return 0;
 436
 437         }
 438         catch(exception& e) {
 439                 m->errorOut(e, "UnifracUnweightedCommand", "execute");
 440                 exit(1);
 441         }
 442 }
 443 /***********************************************************/
 444 void UnifracUnweightedCommand::printUnweightedFile() {
 445         try {
 446                 vector<double> data;
 447                 vector<string> tags;
 448
 449                 tags.push_back("Score");
 450                 tags.push_back("RandFreq"); tags.push_back("RandCumul");
 451
 452                 for(int a = 0; a < numComp; a++) {
 453                         output->initFile(groupComb[a], tags);
 454                         //print each line
 455                         for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 456                                 data.push_back(it->first);  data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 457                                 output->output(data);
 458                                 data.clear();
 459                         }
 460                         output->resetFile();
 461                 }
 462         }
 463         catch(exception& e) {
 464                 m->errorOut(e, "UnifracUnweightedCommand", "printUnweightedFile");
 465                 exit(1);
 466         }
 467 }
 468
 469 /***********************************************************/
 470 void UnifracUnweightedCommand::printUWSummaryFile(int i) {
 471         try {
 472
 473                 //format output
 474                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 475
 476                 //print each line
 477
 478                 for(int a = 0; a < numComp; a++) {
 479                         outSum << i+1 << '\t';
 480                         m->mothurOut(toString(i+1) + "\t");
 481
 482                         if (random) {
 483                                 if (UWScoreSig[a][0] > (1/(float)iters)) {
 484                                         outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
 485                                         cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
 486                                         m->mothurOutJustToLog(groupComb[a]  + "\t" + toString(utreeScores[a][0])  + "\t" + toString(UWScoreSig[a][0])+ "\n");
 487                                 }else {
 488                                         outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 489                                         cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
 490                                         m->mothurOutJustToLog(groupComb[a]  + "\t" + toString(utreeScores[a][0])  + "\t<" + toString((1/float(iters))) + "\n");
 491                                 }
 492                         }else{
 493                                 outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0]  << endl;
 494                                 cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0]  << endl;
 495                                 m->mothurOutJustToLog(groupComb[a]  + "\t" + toString(utreeScores[a][0]) + "\n");
 496                         }
 497                 }
 498
 499         }
 500         catch(exception& e) {
 501                 m->errorOut(e, "UnifracUnweightedCommand", "printUWSummaryFile");
 502                 exit(1);
 503         }
 504 }
 505 /***********************************************************/
 506 void UnifracUnweightedCommand::createPhylipFile(int i) {
 507         try {
 508                 string phylipFileName;
 509                 if ((outputForm == "lt") || (outputForm == "square")) {
 510                         phylipFileName = outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted.phylip.dist";
 511                         outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
 512                 }else { //column
 513                         phylipFileName = outputDir + m->getSimpleName(treefile)  + toString(i+1) + ".unweighted.column.dist";
 514                         outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
 515                 }
 516
 517                 ofstream out;
 518                 m->openOutputFile(phylipFileName, out);
 519
 520                 if ((outputForm == "lt") || (outputForm == "square")) {
 521                         //output numSeqs
 522                         out << m->getNumGroups() << endl;
 523                 }
 524
 525                 //make matrix with scores in it
 526                 vector< vector<float> > dists;  dists.resize(m->getNumGroups());
 527                 for (int i = 0; i < m->getNumGroups(); i++) {
 528                         dists[i].resize(m->getNumGroups(), 0.0);
 529                 }
 530
 531                 //flip it so you can print it
 532                 int count = 0;
 533                 for (int r=0; r<m->getNumGroups(); r++) {
 534                         for (int l = 0; l < r; l++) {
 535                                 dists[r][l] = utreeScores[count][0];
 536                                 dists[l][r] = utreeScores[count][0];
 537                                 count++;
 538                         }
 539                 }
 540
 541                 //output to file
 542                 for (int r=0; r<m->getNumGroups(); r++) {
 543                         //output name
 544                         string name = (m->getGroups())[r];
 545                         if (name.length() < 10) { //pad with spaces to make compatible
 546                                 while (name.length() < 10) {  name += " ";  }
 547                         }
 548
 549                         if (outputForm == "lt") {
 550                                 out << name << '\t';
 551
 552                                 //output distances
 553                                 for (int l = 0; l < r; l++) {   out  << dists[r][l] << '\t';  }
 554                                 out << endl;
 555                         }else if (outputForm == "square") {
 556                                 out << name << '\t';
 557
 558                                 //output distances
 559                                 for (int l = 0; l < m->getNumGroups(); l++) {   out << dists[r][l] << '\t';  }
 560                                 out << endl;
 561                         }else{
 562                                 //output distances
 563                                 for (int l = 0; l < r; l++) {
 564                                         string otherName = (m->getGroups())[l];
 565                                         if (otherName.length() < 10) { //pad with spaces to make compatible
 566                                                 while (otherName.length() < 10) {  otherName += " ";  }
 567                                         }
 568
 569                                         out  << name << '\t' << otherName << dists[r][l] << endl;
 570                                 }
 571                         }
 572                 }
 573                 out.close();
 574         }
 575         catch(exception& e) {
 576                 m->errorOut(e, "UnifracUnweightedCommand", "createPhylipFile");
 577                 exit(1);
 578         }
 579 }/*****************************************************************/
 580 int UnifracUnweightedCommand::readNamesFile() {
 581         try {
 582                 m->names.clear();
 583                 numUniquesInName = 0;
 584
 585                 ifstream in;
 586                 m->openInputFile(namefile, in);
 587
 588                 string first, second;
 589                 map<string, string>::iterator itNames;
 590
 591                 while(!in.eof()) {
 592                         in >> first >> second; m->gobble(in);
 593
 594                         numUniquesInName++;
 595
 596                         itNames = m->names.find(first);
 597                         if (itNames == m->names.end()) {
 598                                 m->names[first] = second;
 599
 600                                 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
 601                                 vector<string> dupNames;
 602                                 m->splitAtComma(second, dupNames);
 603
 604                                 for (int i = 0; i < dupNames.size(); i++) {
 605                                         nameMap[dupNames[i]] = dupNames[i];
 606                                         if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
 607                                 }
 608                         }else {  m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
 609                 }
 610                 in.close();
 611
 612                 return 0;
 613         }
 614         catch(exception& e) {
 615                 m->errorOut(e, "UnifracUnweightedCommand", "readNamesFile");
 616                 exit(1);
 617         }
 618 }
 619 /***********************************************************/
 620
 621
 622
 623