shhhseqscommand.cpp

   1 /*
   2  *  shhhseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 11/8/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "shhhseqscommand.h"
  11
  12
  13
  14 //**********************************************************************************************************************
  15 vector<string> ShhhSeqsCommand::setParameters(){
  16         try {
  17                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
  18                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pname);
  19                 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
  20                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  21                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  22                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  23                 CommandParameter psigma("sigma", "Number", "", "0.01", "", "", "",false,false); parameters.push_back(psigma);
  24
  25                 vector<string> myArray;
  26                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  27                 return myArray;
  28         }
  29         catch(exception& e) {
  30                 m->errorOut(e, "ShhhSeqsCommand", "setParameters");
  31                 exit(1);
  32         }
  33 }
  34 //**********************************************************************************************************************
  35 string ShhhSeqsCommand::getHelpString(){
  36         try {
  37                 string helpString = "";
  38                 helpString += "The shhh.seqs command reads a fasta and name file and ....\n";
  39                 helpString += "The shhh.seqs command parameters are fasta, name, group, sigma and processors.\n";
  40                 helpString += "The fasta parameter allows you to enter the fasta file containing your potentially sequences, and is required, unless you have a valid current fasta file. \n";
  41                 helpString += "The name parameter allows you to provide a name file associated with your fasta file. It is required. \n";
  42                 helpString += "The group parameter allows you to provide a group file.  When checking sequences, only sequences from the same group as the query sequence will be used as the reference. \n";
  43                 helpString += "The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n";
  44                 helpString += "The sigma parameter ....  The default is 0.01. \n";
  45                 helpString += "The shhh.seqs command should be in the following format: \n";
  46                 helpString += "shhh.seqs(fasta=yourFastaFile, name=yourNameFile) \n";
  47                 helpString += "Example: shhh.seqs(fasta=AD.align, name=AD.names) \n";
  48                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n";
  49                 return helpString;
  50
  51         }
  52         catch(exception& e) {
  53                 m->errorOut(e, "ShhhSeqsCommand", "getHelpString");
  54                 exit(1);
  55         }
  56 }
  57 //**********************************************************************************************************************
  58
  59 ShhhSeqsCommand::ShhhSeqsCommand(){
  60         try {
  61                 abort = true; calledHelp = true;
  62                 setParameters();
  63                 vector<string> tempOutNames;
  64                 outputTypes["fasta"] = tempOutNames;
  65                 outputTypes["name"] = tempOutNames;
  66                 outputTypes["map"] = tempOutNames;
  67         }
  68         catch(exception& e) {
  69                 m->errorOut(e, "ShhhSeqsCommand", "ShhhSeqsCommand");
  70                 exit(1);
  71         }
  72 }
  73
  74 //**********************************************************************************************************************
  75 ShhhSeqsCommand::ShhhSeqsCommand(string option) {
  76         try {
  77                 abort = false; calledHelp = false;
  78
  79                 //allow user to run help
  80                 if(option == "help") { help(); abort = true; calledHelp = true; }
  81                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  82
  83                 else {
  84                         vector<string> myArray = setParameters();
  85
  86                         OptionParser parser(option);
  87                         map<string, string> parameters = parser.getParameters();
  88
  89                         ValidParameters validParameter;
  90                         map<string, string>::iterator it;
  91
  92                         //check to make sure all parameters are valid for command
  93                         for (map<string, string>::iterator it2 = parameters.begin(); it2 != parameters.end(); it2++) {
  94                                 if (validParameter.isValidParameter(it2->first, myArray, it2->second) != true) {  abort = true;  }
  95                         }
  96
  97                         //initialize outputTypes
  98                         vector<string> tempOutNames;
  99                         outputTypes["fasta"] = tempOutNames;
 100                         outputTypes["name"] = tempOutNames;
 101                         outputTypes["map"] = tempOutNames;
 102
 103                         //if the user changes the input directory command factory will send this info to us in the output parameter
 104                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 105                         if (inputDir == "not found"){   inputDir = "";          }
 106                         else {
 107                                 string path;
 108                                 it = parameters.find("fasta");
 109                                 //user has given a template file
 110                                 if(it != parameters.end()){
 111                                         path = m->hasPath(it->second);
 112                                         //if the user has not given a path then, add inputdir. else leave path alone.
 113                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
 114                                 }
 115
 116                                 it = parameters.find("name");
 117                                 //user has given a template file
 118                                 if(it != parameters.end()){
 119                                         path = m->hasPath(it->second);
 120                                         //if the user has not given a path then, add inputdir. else leave path alone.
 121                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 122                                 }
 123
 124                                 it = parameters.find("group");
 125                                 //user has given a template file
 126                                 if(it != parameters.end()){
 127                                         path = m->hasPath(it->second);
 128                                         //if the user has not given a path then, add inputdir. else leave path alone.
 129                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 130                                 }
 131                         }
 132
 133                         //check for required parameters
 134                         fastafile = validParameter.validFile(parameters, "fasta", true);
 135                         if (fastafile == "not found") {
 136                                 fastafile = m->getFastaFile();
 137                                 if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
 138                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
 139                         }
 140                         else if (fastafile == "not open") { abort = true; }
 141                         else { m->setFastaFile(fastafile); }
 142
 143                         //if the user changes the output directory command factory will send this info to us in the output parameter
 144                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 145
 146                         //check for optional parameter and set defaults
 147                         // ...at some point should added some additional type checking...
 148                         namefile = validParameter.validFile(parameters, "name", true);
 149                         if (namefile == "not found") {
 150                                 namefile = m->getNameFile();
 151                                 if (namefile != "") { m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
 152                                 else {  m->mothurOut("You have no current namefile and the name parameter is required."); m->mothurOutEndLine(); abort = true; }
 153                         }
 154                         else if (namefile == "not open") { namefile =  ""; abort = true; }
 155                         else {  m->setNameFile(namefile); }
 156
 157                         groupfile = validParameter.validFile(parameters, "group", true);
 158                         if (groupfile == "not found") { groupfile =  "";   }
 159                         else if (groupfile == "not open") { abort = true; groupfile =  ""; }
 160                         else {   m->setGroupFile(groupfile);  }
 161
 162                         string temp     = validParameter.validFile(parameters, "sigma", false);         if(temp == "not found"){        temp = "0.01"; }
 163                         convert(temp, sigma);
 164
 165                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 166                         m->setProcessors(temp);
 167                         convert(temp, processors);
 168                 }
 169         }
 170         catch(exception& e) {
 171                 m->errorOut(e, "ShhhSeqsCommand", "ShhhSeqsCommand");
 172                 exit(1);
 173         }
 174 }
 175 //**********************************************************************************************************************
 176 int ShhhSeqsCommand::execute() {
 177         try {
 178
 179                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 180
 181                 if (outputDir == "") { outputDir = m->hasPath(fastafile);  }//if user entered a file with a path then preserve it
 182                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "shhh.fasta";
 183                 string nameFileName = outputDir + m->getRootName(m->getSimpleName(fastafile))  + "shhh.names";
 184                 string mapFileName = outputDir + m->getRootName(m->getSimpleName(fastafile))  + "shhh.map";
 185
 186                 if (groupfile != "") {
 187                         //Parse sequences by group
 188                         SequenceParser parser(groupfile, fastafile, namefile);
 189                         vector<string> groups = parser.getNamesOfGroups();
 190
 191                         if (m->control_pressed) {  return 0; }
 192
 193                         //clears files
 194                         ofstream out, out1, out2;
 195                         m->openOutputFile(outputFileName, out); out.close();
 196                         m->openOutputFile(nameFileName, out1); out1.close();
 197                         mapFileName = outputDir + m->getRootName(m->getSimpleName(fastafile))  + "shhh.";
 198
 199                         if(processors == 1)     {       driverGroups(parser, outputFileName, nameFileName, mapFileName, 0, groups.size(), groups);      }
 200                         else                            {       createProcessesGroups(parser, outputFileName, nameFileName, mapFileName, groups);                       }
 201
 202                         if (m->control_pressed) {    return 0;  }
 203
 204                         //deconvolute results by running unique.seqs
 205
 206
 207                         if (m->control_pressed) {   return 0;   }
 208
 209                 }else{
 210                         vector<string> sequences;
 211                         vector<string> uniqueNames;
 212                         vector<string> redundantNames;
 213                         vector<int> seqFreq;
 214
 215                         seqNoise noise;
 216                         correctDist* correct = new correctDist(processors);
 217
 218                         //reads fasta and name file and loads them in order
 219                         readData(correct, noise, sequences, uniqueNames, redundantNames, seqFreq);
 220                         if (m->control_pressed) { return 0; }
 221
 222                         //calc distances for cluster
 223                         string distFileName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "shhh.dist";
 224                         correct->execute(distFileName);
 225                         delete correct;
 226
 227                         if (m->control_pressed) { m->mothurRemove(distFileName); return 0; }
 228
 229                         driver(noise, sequences, uniqueNames, redundantNames, seqFreq, distFileName, outputFileName, nameFileName, mapFileName);
 230                 }
 231
 232                 if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) {        m->mothurRemove(outputNames[j]);        } return 0; }
 233
 234                 outputNames.push_back(outputFileName); outputTypes["fasta"].push_back(outputFileName);
 235                 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
 236                 outputNames.push_back(mapFileName); outputTypes["map"].push_back(mapFileName);
 237
 238                 m->mothurOutEndLine();
 239                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 240                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 241                 m->mothurOutEndLine();
 242
 243                 //set accnos file as new current accnosfile
 244                 string current = "";
 245                 itTypes = outputTypes.find("fasta");
 246                 if (itTypes != outputTypes.end()) {
 247                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
 248                 }
 249
 250                 itTypes = outputTypes.find("name");
 251                 if (itTypes != outputTypes.end()) {
 252                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
 253                 }
 254
 255
 256                 return 0;
 257         }
 258         catch(exception& e) {
 259                 m->errorOut(e, "ShhhSeqsCommand", "execute");
 260                 exit(1);
 261         }
 262 }
 263 //**********************************************************************************************************************
 264 int ShhhSeqsCommand::readData(correctDist* correct, seqNoise& noise, vector<string>& seqs, vector<string>& uNames, vector<string>& rNames, vector<int>& freq) {
 265         try {
 266                 map<string, string> nameMap;
 267                 map<string, string>::iterator it;
 268                 m->readNames(namefile, nameMap);
 269                 bool error = false;
 270
 271                 ifstream in;
 272                 m->openInputFile(fastafile, in);
 273
 274                 while (!in.eof()) {
 275
 276                         if (m->control_pressed) { in.close(); return 0; }
 277
 278                         Sequence seq(in); m->gobble(in);
 279
 280                         if (seq.getName() != "") {
 281                                 correct->addSeq(seq.getName(), seq.getAligned());
 282
 283                                 it = nameMap.find(seq.getName());
 284                                 if (it != nameMap.end()) {
 285                                         noise.addSeq(seq.getAligned(), seqs);
 286                                         noise.addRedundantName(it->first, it->second, uNames, rNames, freq);
 287                                 }else {
 288                                         m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file and not in your namefile, please correct.");
 289                                         error = true;
 290                                 }
 291                         }
 292                 }
 293                 in.close();
 294
 295                 if (error) { m->control_pressed = true; }
 296
 297                 return seqs.size();
 298
 299         }catch(exception& e) {
 300                 m->errorOut(e, "ShhhSeqsCommand", "readData");
 301                 exit(1);
 302         }
 303 }
 304 //**********************************************************************************************************************
 305 int ShhhSeqsCommand::loadData(correctDist* correct, seqNoise& noise, vector<string>& seqs, vector<string>& uNames, vector<string>& rNames, vector<int>& freq, map<string, string>& nameMap, vector<Sequence>& sequences) {
 306         try {
 307                 bool error = false;
 308                 map<string, string>::iterator it;
 309
 310                 for (int i = 0; i < sequences.size(); i++) {
 311
 312                         if (m->control_pressed) { return 0; }
 313
 314                         if (sequences[i].getName() != "") {
 315                                 correct->addSeq(sequences[i].getName(), sequences[i].getAligned());
 316
 317                                 it = nameMap.find(sequences[i].getName());
 318                                 if (it != nameMap.end()) {
 319                                         noise.addSeq(sequences[i].getAligned(), seqs);
 320                                         noise.addRedundantName(it->first, it->second, uNames, rNames, freq);
 321                                 }else {
 322                                         m->mothurOut("[ERROR]: " + sequences[i].getName() + " is in your fasta file and not in your namefile, please correct.");
 323                                         error = true;
 324                                 }
 325                         }
 326                 }
 327
 328                 if (error) { m->control_pressed = true; }
 329
 330                 return seqs.size();
 331
 332         }catch(exception& e) {
 333                 m->errorOut(e, "ShhhSeqsCommand", "loadData");
 334                 exit(1);
 335         }
 336 }
 337 /**************************************************************************************************/
 338 int ShhhSeqsCommand::createProcessesGroups(SequenceParser& parser, string newFName, string newNName, string newMName, vector<string> groups) {
 339         try {
 340
 341                 vector<int> processIDS;
 342                 int process = 1;
 343
 344                 //sanity check
 345                 if (groups.size() < processors) { processors = groups.size(); }
 346
 347                 //divide the groups between the processors
 348                 vector<linePair> lines;
 349                 int numGroupsPerProcessor = groups.size() / processors;
 350                 for (int i = 0; i < processors; i++) {
 351                         int startIndex =  i * numGroupsPerProcessor;
 352                         int endIndex = (i+1) * numGroupsPerProcessor;
 353                         if(i == (processors - 1)){      endIndex = groups.size();       }
 354                         lines.push_back(linePair(startIndex, endIndex));
 355                 }
 356
 357 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 358
 359                 //loop through and create all the processes you want
 360                 while (process != processors) {
 361                         int pid = fork();
 362
 363                         if (pid > 0) {
 364                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 365                                 process++;
 366                         }else if (pid == 0){
 367                                 driverGroups(parser, newFName + toString(getpid()) + ".temp", newNName + toString(getpid()) + ".temp", newMName, lines[process].start, lines[process].end, groups);
 368                                 exit(0);
 369                         }else {
 370                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 371                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 372                                 exit(0);
 373                         }
 374                 }
 375
 376                 //do my part
 377                 driverGroups(parser, newFName, newNName, newMName, lines[0].start, lines[0].end, groups);
 378
 379                 //force parent to wait until all the processes are done
 380                 for (int i=0;i<processIDS.size();i++) {
 381                         int temp = processIDS[i];
 382                         wait(&temp);
 383                 }
 384
 385 #else
 386
 387                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 388                 //Windows version shared memory, so be careful when passing variables through the shhhseqsData struct.
 389                 //Above fork() will clone, so memory is separate, but that's not the case with windows,
 390                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 391
 392                 vector<shhhseqsData*> pDataArray;
 393                 DWORD   dwThreadIdArray[processors-1];
 394                 HANDLE  hThreadArray[processors-1];
 395
 396                 //Create processor worker threads.
 397                 for( int i=1; i<processors; i++ ){
 398                         // Allocate memory for thread data.
 399                         string extension = toString(i) + ".temp";
 400
 401                         shhhseqsData* tempShhhseqs = new shhhseqsData(fastafile, namefile, groupfile, (newFName+extension), (newNName+extension), newMName, groups, m, lines[i].start, lines[i].end, sigma, i);
 402                         pDataArray.push_back(tempShhhseqs);
 403                         processIDS.push_back(i);
 404
 405                         //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
 406                         //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
 407                         hThreadArray[i-1] = CreateThread(NULL, 0, MyShhhSeqsThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
 408                 }
 409
 410
 411                 //using the main process as a worker saves time and memory
 412                 driverGroups(parser, newFName, newNName, newMName, lines[0].start, lines[0].end, groups);
 413
 414                 //Wait until all threads have terminated.
 415                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 416
 417                 //Close all thread handles and free memory allocations.
 418                 for(int i=0; i < pDataArray.size(); i++){
 419                         CloseHandle(hThreadArray[i]);
 420                         delete pDataArray[i];
 421                 }
 422
 423 #endif
 424
 425                 //append output files
 426                 for(int i=0;i<processIDS.size();i++){
 427                         m->appendFiles((newFName + toString(processIDS[i]) + ".temp"), newFName);
 428                         m->mothurRemove((newFName + toString(processIDS[i]) + ".temp"));
 429
 430                         m->appendFiles((newNName + toString(processIDS[i]) + ".temp"), newNName);
 431                         m->mothurRemove((newNName + toString(processIDS[i]) + ".temp"));
 432                 }
 433
 434                 return 0;
 435
 436         }
 437         catch(exception& e) {
 438                 m->errorOut(e, "ShhhSeqsCommand", "createProcessesGroups");
 439                 exit(1);
 440         }
 441 }
 442 /**************************************************************************************************/
 443 int ShhhSeqsCommand::driverGroups(SequenceParser& parser, string newFFile, string newNFile, string newMFile, int start, int end, vector<string> groups){
 444         try {
 445
 446                 for (int i = start; i < end; i++) {
 447
 448                         start = time(NULL);
 449
 450                         if (m->control_pressed) {  return 0; }
 451
 452                         m->mothurOutEndLine(); m->mothurOut("Processing group " + groups[i] + ":"); m->mothurOutEndLine();
 453
 454                         map<string, string> thisNameMap;
 455                         thisNameMap = parser.getNameMap(groups[i]);
 456                         vector<Sequence> thisSeqs = parser.getSeqs(groups[i]);
 457
 458                         vector<string> sequences;
 459                         vector<string> uniqueNames;
 460                         vector<string> redundantNames;
 461                         vector<int> seqFreq;
 462
 463                         seqNoise noise;
 464                         correctDist* correct = new correctDist(1); //we use one processor since we already split up the work load.
 465
 466                         //load this groups info in order
 467                         loadData(correct, noise, sequences, uniqueNames, redundantNames, seqFreq, thisNameMap, thisSeqs);
 468                         if (m->control_pressed) { return 0; }
 469
 470                         //calc distances for cluster
 471                         string distFileName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + groups[i] + ".shhh.dist";
 472                         correct->execute(distFileName);
 473                         delete correct;
 474
 475                         if (m->control_pressed) { m->mothurRemove(distFileName); return 0; }
 476
 477                         driver(noise, sequences, uniqueNames, redundantNames, seqFreq, distFileName, newFFile+groups[i], newNFile+groups[i], newMFile+groups[i]+".map");
 478
 479                         if (m->control_pressed) { return 0; }
 480
 481                         m->appendFiles(newFFile+groups[i], newFFile); m->mothurRemove(newFFile+groups[i]);
 482                         m->appendFiles(newNFile+groups[i], newNFile); m->mothurRemove(newNFile+groups[i]);
 483
 484                         m->mothurOut("It took " + toString(time(NULL) - start) + " secs to process group " + groups[i] + "."); m->mothurOutEndLine();
 485                 }
 486
 487                 return 0;
 488         }
 489         catch(exception& e) {
 490                 m->errorOut(e, "ShhhSeqsCommand", "driverGroups");
 491                 exit(1);
 492         }
 493 }
 494 //**********************************************************************************************************************
 495 int ShhhSeqsCommand::driver(seqNoise& noise,
 496                                                         vector<string>& sequences,
 497                                                         vector<string>& uniqueNames,
 498                                                         vector<string>& redundantNames,
 499                                                         vector<int>& seqFreq,
 500                                                         string distFileName, string outputFileName, string nameFileName, string mapFileName) {
 501         try {
 502                 double cutOff = 0.08;
 503                 int minIter = 10;
 504                 int maxIter = 1000;
 505                 double minDelta = 1e-6;
 506                 int numIters = 0;
 507                 double maxDelta = 1e6;
 508                 int numSeqs = sequences.size();
 509
 510                 //run cluster command
 511                 string inputString = "phylip=" + distFileName + ", method=furthest, cutoff=0.08";
 512                 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 513                 m->mothurOut("Running command: cluster(" + inputString + ")"); m->mothurOutEndLine();
 514
 515                 Command* clusterCommand = new ClusterCommand(inputString);
 516                 clusterCommand->execute();
 517
 518                 map<string, vector<string> > filenames = clusterCommand->getOutputFiles();
 519                 string listFileName = filenames["list"][0];
 520                 string rabundFileName = filenames["rabund"][0]; m->mothurRemove(rabundFileName);
 521                 string sabundFileName = filenames["sabund"][0]; m->mothurRemove(sabundFileName);
 522
 523                 delete clusterCommand;
 524                 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 525
 526                 if (m->control_pressed) { m->mothurRemove(distFileName); m->mothurRemove(listFileName); return 0; }
 527
 528                 vector<double> distances(numSeqs * numSeqs);
 529                 noise.getDistanceData(distFileName, distances);
 530                 m->mothurRemove(distFileName);
 531                 if (m->control_pressed) { m->mothurRemove(listFileName); return 0; }
 532
 533                 vector<int> otuData(numSeqs);
 534                 vector<int> otuFreq;
 535                 vector<vector<int> > otuBySeqLookUp;
 536                 noise.getListData(listFileName, cutOff, otuData, otuFreq, otuBySeqLookUp);
 537                 m->mothurRemove(listFileName);
 538                 if (m->control_pressed) { return 0; }
 539
 540                 int numOTUs = otuFreq.size();
 541
 542                 vector<double> weights(numOTUs, 0);
 543                 vector<int> change(numOTUs, 1);
 544                 vector<int> centroids(numOTUs, -1);
 545                 vector<int> cumCount(numOTUs, 0);
 546
 547                 vector<double> tau(numSeqs, 1);
 548                 vector<int> anP(numSeqs, 0);
 549                 vector<int> anI(numSeqs, 0);
 550                 vector<int> anN(numSeqs, 0);
 551                 vector<vector<int> > aanI = otuBySeqLookUp;
 552
 553                 while(numIters < minIter || ((maxDelta > minDelta) && (numIters < maxIter))){
 554
 555                         if (m->control_pressed) { return 0; }
 556
 557                         noise.updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount); if (m->control_pressed) { return 0; }
 558                         maxDelta = noise.calcNewWeights(weights, seqFreq, anI, cumCount, anP, otuFreq, tau);  if (m->control_pressed) { return 0; }
 559
 560                         noise.calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau); if (m->control_pressed) { return 0; }
 561                         noise.checkCentroids(weights, centroids); if (m->control_pressed) { return 0; }
 562
 563                         otuFreq.assign(numOTUs, 0);
 564
 565                         int total = 0;
 566
 567                         for(int i=0;i<numSeqs;i++){
 568                                 if (m->control_pressed) { return 0; }
 569
 570                                 double offset = 1e6;
 571                                 double norm = 0.0000;
 572                                 double minWeight = 0.1;
 573                                 vector<double> currentTau(numOTUs);
 574
 575                                 for(int j=0;j<numOTUs;j++){
 576                                         if (m->control_pressed) { return 0; }
 577                                         if(weights[j] > minWeight && distances[i * numSeqs+centroids[j]] < offset){
 578                                                 offset = distances[i * numSeqs+centroids[j]];
 579                                         }
 580                                 }
 581
 582                                 for(int j=0;j<numOTUs;j++){
 583                                         if (m->control_pressed) { return 0; }
 584                                         if(weights[j] > minWeight){
 585                                                 currentTau[j] = exp(sigma * (-distances[(i * numSeqs + centroids[j])] + offset)) * weights[j];
 586                                                 norm += currentTau[j];
 587                                         }
 588                                         else{
 589                                                 currentTau[j] = 0.0000;
 590                                         }
 591                                 }
 592
 593                                 for(int j=0;j<numOTUs;j++){
 594                                         if (m->control_pressed) { return 0; }
 595                                         currentTau[j] /= norm;
 596                                 }
 597
 598                                 for(int j=0;j<numOTUs;j++){
 599                                         if (m->control_pressed) { return 0; }
 600
 601                                         if(currentTau[j] > 1.0e-4){
 602                                                 int oldTotal = total;
 603                                                 total++;
 604
 605                                                 tau.resize(oldTotal+1);
 606                                                 tau[oldTotal] = currentTau[j];
 607                                                 otuBySeqLookUp[j][otuFreq[j]] = oldTotal;
 608                                                 aanI[j][otuFreq[j]] = i;
 609                                                 otuFreq[j]++;
 610
 611                                         }
 612                                 }
 613
 614                                 anP.resize(total);
 615                                 anI.resize(total);
 616                         }
 617
 618                         numIters++;
 619                 }
 620
 621                 noise.updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);  if (m->control_pressed) { return 0; }
 622
 623                 vector<double> percentage(numSeqs);
 624                 noise.setUpOTUData(otuData, percentage, cumCount, tau, otuFreq, anP, anI);  if (m->control_pressed) { return 0; }
 625                 noise.finishOTUData(otuData, otuFreq, anP, anI, cumCount, otuBySeqLookUp, aanI, tau);  if (m->control_pressed) { return 0; }
 626
 627                 change.assign(numOTUs, 1);
 628                 noise.calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau); if (m->control_pressed) { return 0; }
 629
 630
 631                 vector<int> finalTau(numOTUs, 0);
 632                 for(int i=0;i<numSeqs;i++){
 633                         if (m->control_pressed) { return 0; }
 634                         finalTau[otuData[i]] += int(seqFreq[i]);
 635                 }
 636
 637                 noise.writeOutput(outputFileName, nameFileName, mapFileName, finalTau, centroids, otuData, sequences, uniqueNames, redundantNames, seqFreq, distances);
 638
 639                 return 0;
 640
 641         }catch(exception& e) {
 642                 m->errorOut(e, "ShhhSeqsCommand", "driver");
 643                 exit(1);
 644         }
 645 }
 646 //**********************************************************************************************************************
 647 int ShhhSeqsCommand::deconvoluteResults(string fastaFile, string nameFile){
 648         try {
 649                 m->mothurOutEndLine(); m->mothurOut("Deconvoluting results:"); m->mothurOutEndLine(); m->mothurOutEndLine();
 650
 651                 //use unique.seqs to create new name and fastafile
 652                 string inputString = "fasta=" + fastaFile + ", name=" + nameFile;
 653                 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 654                 m->mothurOut("Running command: unique.seqs(" + inputString + ")"); m->mothurOutEndLine();
 655
 656                 Command* uniqueCommand = new DeconvoluteCommand(inputString);
 657                 uniqueCommand->execute();
 658
 659                 map<string, vector<string> > filenames = uniqueCommand->getOutputFiles();
 660
 661                 delete uniqueCommand;
 662
 663                 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 664
 665                 string newnameFile = filenames["name"][0];
 666                 string newfastaFile = filenames["fasta"][0];
 667
 668                 m->mothurRemove(fastaFile); rename(newfastaFile.c_str(), fastaFile.c_str());
 669                 m->mothurRemove(nameFile); rename(newnameFile.c_str(), nameFile.c_str());
 670
 671                 return 0;
 672         }
 673         catch(exception& e) {
 674                 m->errorOut(e, "ShhhSeqsCommand", "deconvoluteResults");
 675                 exit(1);
 676         }
 677 }
 678 //**********************************************************************************************************************
 679
 680
 681