unifracunweightedcommand.cpp

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