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