chimeraperseuscommand.cpp

   1 /*
   2  *  chimeraperseuscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 10/26/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "chimeraperseuscommand.h"
  11 #include "deconvolutecommand.h"
  12 #include "sequence.hpp"
  13 //**********************************************************************************************************************
  14 vector<string> ChimeraPerseusCommand::setParameters(){
  15         try {
  16                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
  17                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pname);
  18                 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
  19                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  20                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  21                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  22                 CommandParameter pcutoff("cutoff", "Number", "", "0.5", "", "", "",false,false); parameters.push_back(pcutoff);
  23                 CommandParameter palpha("alpha", "Number", "", "-5.54", "", "", "",false,false); parameters.push_back(palpha);
  24                 CommandParameter pbeta("beta", "Number", "", "0.33", "", "", "",false,false); parameters.push_back(pbeta);
  25
  26                 vector<string> myArray;
  27                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  28                 return myArray;
  29         }
  30         catch(exception& e) {
  31                 m->errorOut(e, "ChimeraPerseusCommand", "setParameters");
  32                 exit(1);
  33         }
  34 }
  35 //**********************************************************************************************************************
  36 string ChimeraPerseusCommand::getHelpString(){
  37         try {
  38                 string helpString = "";
  39                 helpString += "The chimera.perseus command reads a fastafile and namefile and outputs potentially chimeric sequences.\n";
  40                 helpString += "The chimera.perseus command parameters are fasta, name, group, cutoff, processors, alpha and beta.\n";
  41                 helpString += "The fasta parameter allows you to enter the fasta file containing your potentially chimeric sequences, and is required, unless you have a valid current fasta file. \n";
  42                 helpString += "The name parameter allows you to provide a name file associated with your fasta file. It is required. \n";
  43                 helpString += "You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amazon.fasta \n";
  44                 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";
  45                 helpString += "The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n";
  46                 helpString += "The alpha parameter ....  The default is -5.54. \n";
  47                 helpString += "The beta parameter ....  The default is 0.33. \n";
  48                 helpString += "The cutoff parameter ....  The default is 0.50. \n";
  49                 helpString += "The chimera.perseus command should be in the following format: \n";
  50                 helpString += "chimera.perseus(fasta=yourFastaFile, name=yourNameFile) \n";
  51                 helpString += "Example: chimera.perseus(fasta=AD.align, name=AD.names) \n";
  52                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n";
  53                 return helpString;
  54         }
  55         catch(exception& e) {
  56                 m->errorOut(e, "ChimeraPerseusCommand", "getHelpString");
  57                 exit(1);
  58         }
  59 }
  60 //**********************************************************************************************************************
  61 ChimeraPerseusCommand::ChimeraPerseusCommand(){
  62         try {
  63                 abort = true; calledHelp = true;
  64                 setParameters();
  65                 vector<string> tempOutNames;
  66                 outputTypes["chimera"] = tempOutNames;
  67                 outputTypes["accnos"] = tempOutNames;
  68         }
  69         catch(exception& e) {
  70                 m->errorOut(e, "ChimeraPerseusCommand", "ChimeraPerseusCommand");
  71                 exit(1);
  72         }
  73 }
  74 //***************************************************************************************************************
  75 ChimeraPerseusCommand::ChimeraPerseusCommand(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("chimera.uchime");
  90                         map<string,string>::iterator it;
  91
  92                         //check to make sure all parameters are valid for command
  93                         for (it = parameters.begin(); it != parameters.end(); it++) {
  94                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  95                         }
  96
  97                         vector<string> tempOutNames;
  98                         outputTypes["chimera"] = tempOutNames;
  99                         outputTypes["accnos"] = tempOutNames;
 100
 101                         //if the user changes the input directory command factory will send this info to us in the output parameter
 102                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 103                         if (inputDir == "not found"){   inputDir = "";          }
 104
 105                         //check for required parameters
 106                         fastafile = validParameter.validFile(parameters, "fasta", false);
 107                         if (fastafile == "not found") {
 108                                 //if there is a current fasta file, use it
 109                                 string filename = m->getFastaFile();
 110                                 if (filename != "") { fastaFileNames.push_back(filename); m->mothurOut("Using " + filename + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
 111                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
 112                         }else {
 113                                 m->splitAtDash(fastafile, fastaFileNames);
 114
 115                                 //go through files and make sure they are good, if not, then disregard them
 116                                 for (int i = 0; i < fastaFileNames.size(); i++) {
 117
 118                                         bool ignore = false;
 119                                         if (fastaFileNames[i] == "current") {
 120                                                 fastaFileNames[i] = m->getFastaFile();
 121                                                 if (fastaFileNames[i] != "") {  m->mothurOut("Using " + fastaFileNames[i] + " as input file for the fasta parameter where you had given current."); m->mothurOutEndLine(); }
 122                                                 else {
 123                                                         m->mothurOut("You have no current fastafile, ignoring current."); m->mothurOutEndLine(); ignore=true;
 124                                                         //erase from file list
 125                                                         fastaFileNames.erase(fastaFileNames.begin()+i);
 126                                                         i--;
 127                                                 }
 128                                         }
 129
 130                                         if (!ignore) {
 131
 132                                                 if (inputDir != "") {
 133                                                         string path = m->hasPath(fastaFileNames[i]);
 134                                                         //if the user has not given a path then, add inputdir. else leave path alone.
 135                                                         if (path == "") {       fastaFileNames[i] = inputDir + fastaFileNames[i];               }
 136                                                 }
 137
 138                                                 int ableToOpen;
 139                                                 ifstream in;
 140
 141                                                 ableToOpen = m->openInputFile(fastaFileNames[i], in, "noerror");
 142
 143                                                 //if you can't open it, try default location
 144                                                 if (ableToOpen == 1) {
 145                                                         if (m->getDefaultPath() != "") { //default path is set
 146                                                                 string tryPath = m->getDefaultPath() + m->getSimpleName(fastaFileNames[i]);
 147                                                                 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
 148                                                                 ifstream in2;
 149                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 150                                                                 in2.close();
 151                                                                 fastaFileNames[i] = tryPath;
 152                                                         }
 153                                                 }
 154
 155                                                 if (ableToOpen == 1) {
 156                                                         if (m->getOutputDir() != "") { //default path is set
 157                                                                 string tryPath = m->getOutputDir() + m->getSimpleName(fastaFileNames[i]);
 158                                                                 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
 159                                                                 ifstream in2;
 160                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 161                                                                 in2.close();
 162                                                                 fastaFileNames[i] = tryPath;
 163                                                         }
 164                                                 }
 165
 166                                                 in.close();
 167
 168                                                 if (ableToOpen == 1) {
 169                                                         m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
 170                                                         //erase from file list
 171                                                         fastaFileNames.erase(fastaFileNames.begin()+i);
 172                                                         i--;
 173                                                 }else {
 174                                                         m->setFastaFile(fastaFileNames[i]);
 175                                                 }
 176                                         }
 177                                 }
 178
 179                                 //make sure there is at least one valid file left
 180                                 if (fastaFileNames.size() == 0) { m->mothurOut("[ERROR]: no valid files."); m->mothurOutEndLine(); abort = true; }
 181                         }
 182
 183
 184                         //check for required parameters
 185                         bool hasName = true;
 186                         namefile = validParameter.validFile(parameters, "name", false);
 187                         if (namefile == "not found") {
 188                                 //if there is a current fasta file, use it
 189                                 string filename = m->getNameFile();
 190                                 if (filename != "") { nameFileNames.push_back(filename); m->mothurOut("Using " + filename + " as input file for the name parameter."); m->mothurOutEndLine(); }
 191                                 else {  m->mothurOut("You have no current namefile and the name parameter is required."); m->mothurOutEndLine(); abort = true; }
 192                                 hasName = false;
 193                         }else {
 194                                 m->splitAtDash(namefile, nameFileNames);
 195
 196                                 //go through files and make sure they are good, if not, then disregard them
 197                                 for (int i = 0; i < nameFileNames.size(); i++) {
 198
 199                                         bool ignore = false;
 200                                         if (nameFileNames[i] == "current") {
 201                                                 nameFileNames[i] = m->getNameFile();
 202                                                 if (nameFileNames[i] != "") {  m->mothurOut("Using " + nameFileNames[i] + " as input file for the name parameter where you had given current."); m->mothurOutEndLine(); }
 203                                                 else {
 204                                                         m->mothurOut("You have no current namefile, ignoring current."); m->mothurOutEndLine(); ignore=true;
 205                                                         //erase from file list
 206                                                         nameFileNames.erase(nameFileNames.begin()+i);
 207                                                         i--;
 208                                                 }
 209                                         }
 210
 211                                         if (!ignore) {
 212
 213                                                 if (inputDir != "") {
 214                                                         string path = m->hasPath(nameFileNames[i]);
 215                                                         //if the user has not given a path then, add inputdir. else leave path alone.
 216                                                         if (path == "") {       nameFileNames[i] = inputDir + nameFileNames[i];         }
 217                                                 }
 218
 219                                                 int ableToOpen;
 220                                                 ifstream in;
 221
 222                                                 ableToOpen = m->openInputFile(nameFileNames[i], in, "noerror");
 223
 224                                                 //if you can't open it, try default location
 225                                                 if (ableToOpen == 1) {
 226                                                         if (m->getDefaultPath() != "") { //default path is set
 227                                                                 string tryPath = m->getDefaultPath() + m->getSimpleName(nameFileNames[i]);
 228                                                                 m->mothurOut("Unable to open " + nameFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
 229                                                                 ifstream in2;
 230                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 231                                                                 in2.close();
 232                                                                 nameFileNames[i] = tryPath;
 233                                                         }
 234                                                 }
 235
 236                                                 if (ableToOpen == 1) {
 237                                                         if (m->getOutputDir() != "") { //default path is set
 238                                                                 string tryPath = m->getOutputDir() + m->getSimpleName(nameFileNames[i]);
 239                                                                 m->mothurOut("Unable to open " + nameFileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
 240                                                                 ifstream in2;
 241                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 242                                                                 in2.close();
 243                                                                 nameFileNames[i] = tryPath;
 244                                                         }
 245                                                 }
 246
 247                                                 in.close();
 248
 249                                                 if (ableToOpen == 1) {
 250                                                         m->mothurOut("Unable to open " + nameFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
 251                                                         //erase from file list
 252                                                         nameFileNames.erase(nameFileNames.begin()+i);
 253                                                         i--;
 254                                                 }else {
 255                                                         m->setNameFile(nameFileNames[i]);
 256                                                 }
 257                                         }
 258                                 }
 259
 260                                 //make sure there is at least one valid file left
 261                                 if (nameFileNames.size() == 0) { m->mothurOut("[ERROR]: no valid name files."); m->mothurOutEndLine(); abort = true; }
 262                         }
 263
 264                         if (hasName && (nameFileNames.size() != fastaFileNames.size())) { m->mothurOut("[ERROR]: The number of namefiles does not match the number of fastafiles, please correct."); m->mothurOutEndLine(); abort=true; }
 265
 266                         bool hasGroup = true;
 267                         groupfile = validParameter.validFile(parameters, "group", false);
 268                         if (groupfile == "not found") { groupfile = "";  hasGroup = false; }
 269                         else {
 270                                 m->splitAtDash(groupfile, groupFileNames);
 271
 272                                 //go through files and make sure they are good, if not, then disregard them
 273                                 for (int i = 0; i < groupFileNames.size(); i++) {
 274
 275                                         bool ignore = false;
 276                                         if (groupFileNames[i] == "current") {
 277                                                 groupFileNames[i] = m->getGroupFile();
 278                                                 if (groupFileNames[i] != "") {  m->mothurOut("Using " + groupFileNames[i] + " as input file for the group parameter where you had given current."); m->mothurOutEndLine(); }
 279                                                 else {
 280                                                         m->mothurOut("You have no current namefile, ignoring current."); m->mothurOutEndLine(); ignore=true;
 281                                                         //erase from file list
 282                                                         groupFileNames.erase(groupFileNames.begin()+i);
 283                                                         i--;
 284                                                 }
 285                                         }
 286
 287                                         if (!ignore) {
 288
 289                                                 if (inputDir != "") {
 290                                                         string path = m->hasPath(groupFileNames[i]);
 291                                                         //if the user has not given a path then, add inputdir. else leave path alone.
 292                                                         if (path == "") {       groupFileNames[i] = inputDir + groupFileNames[i];               }
 293                                                 }
 294
 295                                                 int ableToOpen;
 296                                                 ifstream in;
 297
 298                                                 ableToOpen = m->openInputFile(groupFileNames[i], in, "noerror");
 299
 300                                                 //if you can't open it, try default location
 301                                                 if (ableToOpen == 1) {
 302                                                         if (m->getDefaultPath() != "") { //default path is set
 303                                                                 string tryPath = m->getDefaultPath() + m->getSimpleName(groupFileNames[i]);
 304                                                                 m->mothurOut("Unable to open " + groupFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
 305                                                                 ifstream in2;
 306                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 307                                                                 in2.close();
 308                                                                 groupFileNames[i] = tryPath;
 309                                                         }
 310                                                 }
 311
 312                                                 if (ableToOpen == 1) {
 313                                                         if (m->getOutputDir() != "") { //default path is set
 314                                                                 string tryPath = m->getOutputDir() + m->getSimpleName(groupFileNames[i]);
 315                                                                 m->mothurOut("Unable to open " + groupFileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
 316                                                                 ifstream in2;
 317                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 318                                                                 in2.close();
 319                                                                 groupFileNames[i] = tryPath;
 320                                                         }
 321                                                 }
 322
 323                                                 in.close();
 324
 325                                                 if (ableToOpen == 1) {
 326                                                         m->mothurOut("Unable to open " + groupFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
 327                                                         //erase from file list
 328                                                         groupFileNames.erase(groupFileNames.begin()+i);
 329                                                         i--;
 330                                                 }else {
 331                                                         m->setGroupFile(groupFileNames[i]);
 332                                                 }
 333                                         }
 334                                 }
 335
 336                                 //make sure there is at least one valid file left
 337                                 if (groupFileNames.size() == 0) { m->mothurOut("[ERROR]: no valid group files."); m->mothurOutEndLine(); abort = true; }
 338                         }
 339
 340                         if (hasGroup && (groupFileNames.size() != fastaFileNames.size())) { m->mothurOut("[ERROR]: The number of groupfiles does not match the number of fastafiles, please correct."); m->mothurOutEndLine(); abort=true; }
 341
 342
 343                         //if the user changes the output directory command factory will send this info to us in the output parameter
 344                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 345
 346                         string temp = validParameter.validFile(parameters, "processors", false);        if (temp == "not found"){       temp = m->getProcessors();      }
 347                         m->setProcessors(temp);
 348                         m->mothurConvert(temp, processors);
 349
 350                         temp = validParameter.validFile(parameters, "cutoff", false);   if (temp == "not found"){       temp = "0.50";  }
 351                         m->mothurConvert(temp, cutoff);
 352
 353                         temp = validParameter.validFile(parameters, "alpha", false);    if (temp == "not found"){       temp = "-5.54"; }
 354                         m->mothurConvert(temp, alpha);
 355
 356                         temp = validParameter.validFile(parameters, "cutoff", false);   if (temp == "not found"){       temp = "0.33";  }
 357                         m->mothurConvert(temp, beta);
 358                 }
 359         }
 360         catch(exception& e) {
 361                 m->errorOut(e, "ChimeraPerseusCommand", "ChimeraPerseusCommand");
 362                 exit(1);
 363         }
 364 }
 365 //***************************************************************************************************************
 366
 367 int ChimeraPerseusCommand::execute(){
 368         try{
 369                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 370
 371
 372                 //process each file
 373                 for (int s = 0; s < fastaFileNames.size(); s++) {
 374
 375                         m->mothurOut("Checking sequences from " + fastaFileNames[s] + " ..." ); m->mothurOutEndLine();
 376
 377                         int start = time(NULL);
 378                         if (outputDir == "") { outputDir = m->hasPath(fastaFileNames[s]);  }//if user entered a file with a path then preserve it
 379                         string outputFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "perseus.chimera";
 380                         string accnosFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s]))  + "perseus.accnos";
 381                         //string newFasta = m->getRootName(fastaFileNames[s]) + "temp";
 382
 383                         //you provided a groupfile
 384                         string groupFile = "";
 385                         if (groupFileNames.size() != 0) { groupFile = groupFileNames[s]; }
 386
 387                         string nameFile = "";
 388                         if (nameFileNames.size() != 0) { //you provided a namefile and we don't need to create one
 389                                 nameFile = nameFileNames[s];
 390                         }else { nameFile = getNamesFile(fastaFileNames[s]); }
 391
 392                         if (m->control_pressed) {  for (int j = 0; j < outputNames.size(); j++) {       m->mothurRemove(outputNames[j]);        } return 0;     }
 393
 394                         int numSeqs = 0;
 395                         int numChimeras = 0;
 396
 397                         if (groupFile != "") {
 398                                 //Parse sequences by group
 399                                 SequenceParser parser(groupFile, fastaFileNames[s], nameFile);
 400                                 vector<string> groups = parser.getNamesOfGroups();
 401
 402                                 if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) {        m->mothurRemove(outputNames[j]);        }  return 0; }
 403
 404                                 //clears files
 405                                 ofstream out, out1, out2;
 406                                 m->openOutputFile(outputFileName, out); out.close();
 407                                 m->openOutputFile(accnosFileName, out1); out1.close();
 408
 409                                 if(processors == 1)     {       numSeqs = driverGroups(parser, outputFileName, accnosFileName, 0, groups.size(), groups);       }
 410                                 else                            {       numSeqs = createProcessesGroups(parser, outputFileName, accnosFileName, groups, groupFile, fastaFileNames[s], nameFile);                        }
 411
 412                                 if (m->control_pressed) {  for (int j = 0; j < outputNames.size(); j++) {       m->mothurRemove(outputNames[j]);        }  return 0;    }
 413
 414                                 numChimeras = deconvoluteResults(parser, outputFileName, accnosFileName);
 415
 416                                 m->mothurOut("The number of sequences checked may be larger than the number of unique sequences because some sequences are found in several samples."); m->mothurOutEndLine();
 417
 418                                 if (m->control_pressed) {  for (int j = 0; j < outputNames.size(); j++) {       m->mothurRemove(outputNames[j]);        }  return 0;    }
 419
 420                         }else{
 421                                 if (processors != 1) { m->mothurOut("Without a groupfile, mothur can only use 1 processor, continuing."); m->mothurOutEndLine(); processors = 1; }
 422
 423                                 //read sequences and store sorted by frequency
 424                                 vector<seqData> sequences = readFiles(fastaFileNames[s], nameFile);
 425
 426                                 if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) {        m->mothurRemove(outputNames[j]);        } return 0; }
 427
 428                                 numSeqs = driver(outputFileName, sequences, accnosFileName, numChimeras);
 429                         }
 430
 431                         if (m->control_pressed) { for (int j = 0; j < outputNames.size(); j++) {        m->mothurRemove(outputNames[j]);        } return 0; }
 432
 433                         m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences. " + toString(numChimeras) + " chimeras were found.");      m->mothurOutEndLine();
 434                         outputNames.push_back(outputFileName); outputTypes["chimera"].push_back(outputFileName);
 435                         outputNames.push_back(accnosFileName); outputTypes["accnos"].push_back(accnosFileName);
 436                 }
 437
 438                 //set accnos file as new current accnosfile
 439                 string current = "";
 440                 itTypes = outputTypes.find("accnos");
 441                 if (itTypes != outputTypes.end()) {
 442                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
 443                 }
 444
 445                 m->mothurOutEndLine();
 446                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 447                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 448                 m->mothurOutEndLine();
 449
 450                 return 0;
 451
 452         }
 453         catch(exception& e) {
 454                 m->errorOut(e, "ChimeraPerseusCommand", "execute");
 455                 exit(1);
 456         }
 457 }
 458 //**********************************************************************************************************************
 459 string ChimeraPerseusCommand::getNamesFile(string& inputFile){
 460         try {
 461                 string nameFile = "";
 462
 463                 m->mothurOutEndLine(); m->mothurOut("No namesfile given, running unique.seqs command to generate one."); m->mothurOutEndLine(); m->mothurOutEndLine();
 464
 465                 //use unique.seqs to create new name and fastafile
 466                 string inputString = "fasta=" + inputFile;
 467                 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 468                 m->mothurOut("Running command: unique.seqs(" + inputString + ")"); m->mothurOutEndLine();
 469                 m->mothurCalling = true;
 470
 471                 Command* uniqueCommand = new DeconvoluteCommand(inputString);
 472                 uniqueCommand->execute();
 473
 474                 map<string, vector<string> > filenames = uniqueCommand->getOutputFiles();
 475
 476                 delete uniqueCommand;
 477                 m->mothurCalling = false;
 478                 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 479
 480                 nameFile = filenames["name"][0];
 481                 inputFile = filenames["fasta"][0];
 482
 483                 return nameFile;
 484         }
 485         catch(exception& e) {
 486                 m->errorOut(e, "ChimeraPerseusCommand", "getNamesFile");
 487                 exit(1);
 488         }
 489 }
 490 //**********************************************************************************************************************
 491 int ChimeraPerseusCommand::driverGroups(SequenceParser& parser, string outputFName, string accnos, int start, int end, vector<string> groups){
 492         try {
 493
 494                 int totalSeqs = 0;
 495                 int numChimeras = 0;
 496
 497                 for (int i = start; i < end; i++) {
 498
 499                         m->mothurOutEndLine(); m->mothurOut("Checking sequences from group " + groups[i] + "...");      m->mothurOutEndLine();
 500
 501                         int start = time(NULL);  if (m->control_pressed) {  return 0; }
 502
 503                         vector<seqData> sequences = loadSequences(parser, groups[i]);
 504
 505                         if (m->control_pressed) { return 0; }
 506
 507                         int numSeqs = driver((outputFName + groups[i]), sequences, (accnos+groups[i]), numChimeras);
 508                         totalSeqs += numSeqs;
 509
 510                         if (m->control_pressed) { return 0; }
 511
 512                         //append files
 513                         m->appendFiles((outputFName+groups[i]), outputFName); m->mothurRemove((outputFName+groups[i]));
 514                         m->appendFiles((accnos+groups[i]), accnos); m->mothurRemove((accnos+groups[i]));
 515
 516                         m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences from group " + groups[i] + ".");    m->mothurOutEndLine();
 517                 }
 518
 519                 return totalSeqs;
 520
 521         }
 522         catch(exception& e) {
 523                 m->errorOut(e, "ChimeraPerseusCommand", "driverGroups");
 524                 exit(1);
 525         }
 526 }
 527 //**********************************************************************************************************************
 528 vector<seqData> ChimeraPerseusCommand::loadSequences(SequenceParser& parser, string group){
 529         try {
 530
 531                 vector<Sequence> thisGroupsSeqs = parser.getSeqs(group);
 532                 map<string, string> nameMap = parser.getNameMap(group);
 533                 map<string, string>::iterator it;
 534
 535                 vector<seqData> sequences;
 536                 bool error = false;
 537         alignLength = 0;
 538
 539                 for (int i = 0; i < thisGroupsSeqs.size(); i++) {
 540
 541                         if (m->control_pressed) {  return sequences; }
 542
 543                         it = nameMap.find(thisGroupsSeqs[i].getName());
 544                         if (it == nameMap.end()) { error = true; m->mothurOut("[ERROR]: " + thisGroupsSeqs[i].getName() + " is in your fasta file and not in your namefile, please correct."); m->mothurOutEndLine(); }
 545                         else {
 546                                 int num = m->getNumNames(it->second);
 547                                 sequences.push_back(seqData(thisGroupsSeqs[i].getName(), thisGroupsSeqs[i].getUnaligned(), num));
 548                 if (thisGroupsSeqs[i].getUnaligned().length() > alignLength) { alignLength = thisGroupsSeqs[i].getUnaligned().length(); }
 549                         }
 550                 }
 551
 552                 if (error) { m->control_pressed = true; }
 553
 554                 //sort by frequency
 555                 sort(sequences.rbegin(), sequences.rend());
 556
 557                 return sequences;
 558         }
 559         catch(exception& e) {
 560                 m->errorOut(e, "ChimeraPerseusCommand", "loadSequences");
 561                 exit(1);
 562         }
 563 }
 564
 565 //**********************************************************************************************************************
 566 vector<seqData> ChimeraPerseusCommand::readFiles(string inputFile, string name){
 567         try {
 568                 map<string, int>::iterator it;
 569                 map<string, int> nameMap = m->readNames(name);
 570
 571                 //read fasta file and create sequenceData structure - checking for file mismatches
 572                 vector<seqData> sequences;
 573                 bool error = false;
 574                 ifstream in;
 575                 m->openInputFile(inputFile, in);
 576                 alignLength = 0;
 577
 578                 while (!in.eof()) {
 579
 580                         if (m->control_pressed) { in.close(); return sequences; }
 581
 582                         Sequence temp(in); m->gobble(in);
 583
 584                         it = nameMap.find(temp.getName());
 585                         if (it == nameMap.end()) { error = true; m->mothurOut("[ERROR]: " + temp.getName() + " is in your fasta file and not in your namefile, please correct."); m->mothurOutEndLine(); }
 586                         else {
 587                                 sequences.push_back(seqData(temp.getName(), temp.getUnaligned(), it->second));
 588                 if (temp.getUnaligned().length() > alignLength) { alignLength = temp.getUnaligned().length(); }
 589                         }
 590                 }
 591                 in.close();
 592
 593                 if (error) { m->control_pressed = true; }
 594
 595                 //sort by frequency
 596                 sort(sequences.rbegin(), sequences.rend());
 597
 598                 return sequences;
 599         }
 600         catch(exception& e) {
 601                 m->errorOut(e, "ChimeraPerseusCommand", "getNamesFile");
 602                 exit(1);
 603         }
 604 }
 605 //**********************************************************************************************************************
 606 int ChimeraPerseusCommand::driver(string chimeraFileName, vector<seqData>& sequences, string accnosFileName, int& numChimeras){
 607         try {
 608
 609                 vector<vector<double> > correctModel(4);        //could be an option in the future to input own model matrix
 610                 for(int i=0;i<4;i++){   correctModel[i].resize(4);      }
 611
 612                 correctModel[0][0] = 0.000000;  //AA
 613                 correctModel[1][0] = 11.619259; //CA
 614                 correctModel[2][0] = 11.694004; //TA
 615                 correctModel[3][0] = 7.748623;  //GA
 616
 617                 correctModel[1][1] = 0.000000;  //CC
 618                 correctModel[2][1] = 7.619657;  //TC
 619                 correctModel[3][1] = 12.852562; //GC
 620
 621                 correctModel[2][2] = 0.000000;  //TT
 622                 correctModel[3][2] = 10.964048; //TG
 623
 624                 correctModel[3][3] = 0.000000;  //GG
 625
 626                 for(int i=0;i<4;i++){
 627                         for(int j=0;j<i;j++){
 628                                 correctModel[j][i] = correctModel[i][j];
 629                         }
 630                 }
 631
 632                 int numSeqs = sequences.size();
 633                 //int alignLength = sequences[0].sequence.size();
 634
 635                 ofstream chimeraFile;
 636                 ofstream accnosFile;
 637                 m->openOutputFile(chimeraFileName, chimeraFile);
 638                 m->openOutputFile(accnosFileName, accnosFile);
 639
 640                 Perseus myPerseus;
 641                 vector<vector<double> > binMatrix = myPerseus.binomial(alignLength);
 642
 643                 chimeraFile << "SequenceIndex\tName\tDiffsToBestMatch\tBestMatchIndex\tBestMatchName\tDiffstToChimera\tIndexofLeftParent\tIndexOfRightParent\tNameOfLeftParent\tNameOfRightParent\tDistanceToBestMatch\tcIndex\t(cIndex - singleDist)\tloonIndex\tMismatchesToChimera\tMismatchToTrimera\tChimeraBreakPoint\tLogisticProbability\tTypeOfSequence\n";
 644
 645                 vector<bool> chimeras(numSeqs, 0);
 646
 647                 for(int i=0;i<numSeqs;i++){
 648                         if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 649
 650                         vector<bool> restricted = chimeras;
 651
 652                         vector<vector<int> > leftDiffs(numSeqs);
 653                         vector<vector<int> > leftMaps(numSeqs);
 654                         vector<vector<int> > rightDiffs(numSeqs);
 655                         vector<vector<int> > rightMaps(numSeqs);
 656
 657                         vector<int> singleLeft, bestLeft;
 658                         vector<int> singleRight, bestRight;
 659
 660                         int bestSingleIndex, bestSingleDiff;
 661                         vector<pwAlign> alignments(numSeqs);
 662
 663                         int comparisons = myPerseus.getAlignments(i, sequences, alignments, leftDiffs, leftMaps, rightDiffs, rightMaps, bestSingleIndex, bestSingleDiff, restricted);
 664                         if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 665
 666                         int minMismatchToChimera, leftParentBi, rightParentBi, breakPointBi;
 667
 668                         string dummyA, dummyB;
 669
 670             if (sequences[i].sequence.size() < 3) {
 671                 chimeraFile << i << '\t' << sequences[i].seqName << "\t0\t0\tNull\t0\t0\t0\tNull\tNull\t0.0\t0.0\t0.0\t0\t0\t0\t0.0\t0.0\tgood" << endl;
 672             }else if(comparisons >= 2){
 673                                 minMismatchToChimera = myPerseus.getChimera(sequences, leftDiffs, rightDiffs, leftParentBi, rightParentBi, breakPointBi, singleLeft, bestLeft, singleRight, bestRight, restricted);
 674                                 if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 675
 676                                 int minMismatchToTrimera = numeric_limits<int>::max();
 677                                 int leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB;
 678
 679                                 if(minMismatchToChimera >= 3 && comparisons >= 3){
 680                                         minMismatchToTrimera = myPerseus.getTrimera(sequences, leftDiffs, leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB, singleLeft, bestLeft, singleRight, bestRight, restricted);
 681                                         if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 682                                 }
 683
 684                                 double singleDist = myPerseus.modeledPairwiseAlignSeqs(sequences[i].sequence, sequences[bestSingleIndex].sequence, dummyA, dummyB, correctModel);
 685
 686                                 if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 687
 688                                 string type;
 689                                 string chimeraRefSeq;
 690
 691                                 if(minMismatchToChimera - minMismatchToTrimera >= 3){
 692                                         type = "trimera";
 693                                         chimeraRefSeq = myPerseus.stitchTrimera(alignments, leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB, leftMaps, rightMaps);
 694                                 }
 695                                 else{
 696                                         type = "chimera";
 697                                         chimeraRefSeq = myPerseus.stitchBimera(alignments, leftParentBi, rightParentBi, breakPointBi, leftMaps, rightMaps);
 698                                 }
 699                                 ;
 700                                 if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 701
 702                                 double chimeraDist = myPerseus.modeledPairwiseAlignSeqs(sequences[i].sequence, chimeraRefSeq, dummyA, dummyB, correctModel);
 703
 704                                 if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 705
 706                                 double cIndex = chimeraDist;//modeledPairwiseAlignSeqs(sequences[i].sequence, chimeraRefSeq);
 707                                 double loonIndex = myPerseus.calcLoonIndex(sequences[i].sequence, sequences[leftParentBi].sequence, sequences[rightParentBi].sequence, breakPointBi, binMatrix);
 708
 709                                 if (m->control_pressed) { chimeraFile.close(); m->mothurRemove(chimeraFileName); accnosFile.close(); m->mothurRemove(accnosFileName); return 0; }
 710
 711                                 chimeraFile << i << '\t' << sequences[i].seqName << '\t' << bestSingleDiff << '\t' << bestSingleIndex << '\t' << sequences[bestSingleIndex].seqName << '\t';
 712                                 chimeraFile << minMismatchToChimera << '\t' << leftParentBi << '\t' << rightParentBi << '\t' << sequences[leftParentBi].seqName << '\t' << sequences[rightParentBi].seqName << '\t';
 713                                 chimeraFile << singleDist << '\t' << cIndex << '\t' << (cIndex - singleDist) << '\t' << loonIndex << '\t';
 714                                 chimeraFile << minMismatchToChimera << '\t' << minMismatchToTrimera << '\t' << breakPointBi << '\t';
 715
 716                                 double probability = myPerseus.classifyChimera(singleDist, cIndex, loonIndex, alpha, beta);
 717
 718                                 chimeraFile << probability << '\t';
 719
 720                                 if(probability > cutoff){
 721                                         chimeraFile << type << endl;
 722                                         accnosFile << sequences[i].seqName << endl;
 723                                         chimeras[i] = 1;
 724                                         numChimeras++;
 725                                 }
 726                                 else{
 727                                         chimeraFile << "good" << endl;
 728                                 }
 729
 730                         }
 731                         else{
 732                                 chimeraFile << i << '\t' << sequences[i].seqName << "\t0\t0\tNull\t0\t0\t0\tNull\tNull\t0.0\t0.0\t0.0\t0\t0\t0\t0.0\t0.0\tgood" << endl;
 733                         }
 734
 735                         //report progress
 736                         if((i+1) % 100 == 0){   m->mothurOut("Processing sequence: " + toString(i+1) + "\n");           }
 737                 }
 738
 739                 if((numSeqs) % 100 != 0){       m->mothurOut("Processing sequence: " + toString(numSeqs) + "\n");               }
 740
 741                 chimeraFile.close();
 742                 accnosFile.close();
 743
 744                 return numSeqs;
 745         }
 746         catch(exception& e) {
 747                 m->errorOut(e, "ChimeraPerseusCommand", "driver");
 748                 exit(1);
 749         }
 750 }
 751 /**************************************************************************************************/
 752 int ChimeraPerseusCommand::createProcessesGroups(SequenceParser& parser, string outputFName, string accnos, vector<string> groups, string group, string fasta, string name) {
 753         try {
 754
 755                 vector<int> processIDS;
 756                 int process = 1;
 757                 int num = 0;
 758
 759                 //sanity check
 760                 if (groups.size() < processors) { processors = groups.size(); }
 761
 762                 //divide the groups between the processors
 763                 vector<linePair> lines;
 764                 int numGroupsPerProcessor = groups.size() / processors;
 765                 for (int i = 0; i < processors; i++) {
 766                         int startIndex =  i * numGroupsPerProcessor;
 767                         int endIndex = (i+1) * numGroupsPerProcessor;
 768                         if(i == (processors - 1)){      endIndex = groups.size();       }
 769                         lines.push_back(linePair(startIndex, endIndex));
 770                 }
 771
 772 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 773
 774                 //loop through and create all the processes you want
 775                 while (process != processors) {
 776                         int pid = fork();
 777
 778                         if (pid > 0) {
 779                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 780                                 process++;
 781                         }else if (pid == 0){
 782                                 num = driverGroups(parser, outputFName + toString(getpid()) + ".temp", accnos + toString(getpid()) + ".temp", lines[process].start, lines[process].end, groups);
 783
 784                                 //pass numSeqs to parent
 785                                 ofstream out;
 786                                 string tempFile = outputFName + toString(getpid()) + ".num.temp";
 787                                 m->openOutputFile(tempFile, out);
 788                                 out << num << endl;
 789                                 out.close();
 790
 791                                 exit(0);
 792                         }else {
 793                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 794                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 795                                 exit(0);
 796                         }
 797                 }
 798
 799                 //do my part
 800                 num = driverGroups(parser, outputFName, accnos, lines[0].start, lines[0].end, groups);
 801
 802                 //force parent to wait until all the processes are done
 803                 for (int i=0;i<processIDS.size();i++) {
 804                         int temp = processIDS[i];
 805                         wait(&temp);
 806                 }
 807
 808                 for (int i = 0; i < processIDS.size(); i++) {
 809                         ifstream in;
 810                         string tempFile =  outputFName + toString(processIDS[i]) + ".num.temp";
 811                         m->openInputFile(tempFile, in);
 812                         if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
 813                         in.close(); m->mothurRemove(tempFile);
 814                 }
 815
 816 #else
 817                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 818                 //Windows version shared memory, so be careful when passing variables through the preClusterData struct.
 819                 //Above fork() will clone, so memory is separate, but that's not the case with windows,
 820                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 821
 822                 vector<perseusData*> pDataArray;
 823                 DWORD   dwThreadIdArray[processors-1];
 824                 HANDLE  hThreadArray[processors-1];
 825
 826                 //Create processor worker threads.
 827                 for( int i=1; i<processors; i++ ){
 828                         // Allocate memory for thread data.
 829                         string extension = toString(i) + ".temp";
 830
 831                         perseusData* tempPerseus = new perseusData(alpha, beta, cutoff, outputFName+extension, fasta, name, group, accnos+extension, groups, m, lines[i].start, lines[i].end, i);
 832
 833                         pDataArray.push_back(tempPerseus);
 834                         processIDS.push_back(i);
 835
 836                         //MyPerseusThreadFunction is in header. It must be global or static to work with the threads.
 837                         //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
 838                         hThreadArray[i-1] = CreateThread(NULL, 0, MyPerseusThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
 839                 }
 840
 841
 842                 //using the main process as a worker saves time and memory
 843                 num = driverGroups(parser, outputFName, accnos, lines[0].start, lines[0].end, groups);
 844
 845                 //Wait until all threads have terminated.
 846                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 847
 848                 //Close all thread handles and free memory allocations.
 849                 for(int i=0; i < pDataArray.size(); i++){
 850                         num += pDataArray[i]->count;
 851                         CloseHandle(hThreadArray[i]);
 852                         delete pDataArray[i];
 853                 }
 854 #endif
 855
 856
 857                 //append output files
 858                 for(int i=0;i<processIDS.size();i++){
 859                         m->appendFiles((outputFName + toString(processIDS[i]) + ".temp"), outputFName);
 860                         m->mothurRemove((outputFName + toString(processIDS[i]) + ".temp"));
 861
 862                         m->appendFiles((accnos + toString(processIDS[i]) + ".temp"), accnos);
 863                         m->mothurRemove((accnos + toString(processIDS[i]) + ".temp"));
 864                 }
 865
 866                 return num;
 867
 868         }
 869         catch(exception& e) {
 870                 m->errorOut(e, "ChimeraPerseusCommand", "createProcessesGroups");
 871                 exit(1);
 872         }
 873 }
 874 //**********************************************************************************************************************
 875 int ChimeraPerseusCommand::deconvoluteResults(SequenceParser& parser, string outputFileName, string accnosFileName){
 876         try {
 877                 map<string, string> uniqueNames = parser.getAllSeqsMap();
 878                 map<string, string>::iterator itUnique;
 879                 int total = 0;
 880
 881                 //edit accnos file
 882                 ifstream in2;
 883                 m->openInputFile(accnosFileName, in2);
 884
 885                 ofstream out2;
 886                 m->openOutputFile(accnosFileName+".temp", out2);
 887
 888                 string name;
 889                 set<string> namesInFile; //this is so if a sequence is found to be chimera in several samples we dont write it to the results file more than once
 890                 set<string>::iterator itNames;
 891                 set<string> chimerasInFile;
 892                 set<string>::iterator itChimeras;
 893
 894
 895                 while (!in2.eof()) {
 896                         if (m->control_pressed) { in2.close(); out2.close(); m->mothurRemove(outputFileName); m->mothurRemove((accnosFileName+".temp")); return 0; }
 897
 898                         in2 >> name; m->gobble(in2);
 899
 900                         //find unique name
 901                         itUnique = uniqueNames.find(name);
 902
 903                         if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing accnos results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; }
 904                         else {
 905                                 itChimeras = chimerasInFile.find((itUnique->second));
 906
 907                                 if (itChimeras == chimerasInFile.end()) {
 908                                         out2 << itUnique->second << endl;
 909                                         chimerasInFile.insert((itUnique->second));
 910                                         total++;
 911                                 }
 912                         }
 913                 }
 914                 in2.close();
 915                 out2.close();
 916
 917                 m->mothurRemove(accnosFileName);
 918                 rename((accnosFileName+".temp").c_str(), accnosFileName.c_str());
 919
 920                 //edit chimera file
 921                 ifstream in;
 922                 m->openInputFile(outputFileName, in);
 923
 924                 ofstream out;
 925                 m->openOutputFile(outputFileName+".temp", out); out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 926
 927                 int DiffsToBestMatch, BestMatchIndex, DiffstToChimera, IndexofLeftParent, IndexOfRightParent;
 928                 float temp1,temp2, temp3, temp4, temp5, temp6, temp7, temp8;
 929                 string index, BestMatchName, parent1, parent2, flag;
 930                 name = "";
 931                 namesInFile.clear();
 932                 //assumptions - in file each read will always look like
 933                 /*
 934                  SequenceIndex  Name    DiffsToBestMatch        BestMatchIndex  BestMatchName   DiffstToChimera IndexofLeftParent       IndexOfRightParent      NameOfLeftParent        NameOfRightParent       DistanceToBestMatch     cIndex  (cIndex - singleDist)   loonIndex       MismatchesToChimera     MismatchToTrimera       ChimeraBreakPoint       LogisticProbability     TypeOfSequence
 935                  0      F01QG4L02JVBQY  0       0       Null    0       0       0       Null    Null    0.0     0.0     0.0     0.0     0       0       0       0.0     0.0     good
 936                  1      F01QG4L02ICTC6  0       0       Null    0       0       0       Null    Null    0.0     0.0     0.0     0.0     0       0       0       0.0     0.0     good
 937                  2      F01QG4L02JZOEC  48      0       F01QG4L02JVBQY  47      0       0       F01QG4L02JVBQY  F01QG4L02JVBQY  2.0449  2.03545 -0.00944493     0       47      2147483647      138     0       good
 938                  3      F01QG4L02G7JEC  42      0       F01QG4L02JVBQY  40      1       0       F01QG4L02ICTC6  F01QG4L02JVBQY  1.87477 1.81113 -0.0636404      5.80145 40      2147483647      25      0       good
 939                  */
 940
 941                 //get and print headers
 942                 BestMatchName = m->getline(in); m->gobble(in);
 943                 out << BestMatchName << endl;
 944
 945                 while (!in.eof()) {
 946
 947                         if (m->control_pressed) { in.close(); out.close(); m->mothurRemove((outputFileName+".temp")); return 0; }
 948
 949                         bool print = false;
 950                         in >> index;    m->gobble(in);
 951
 952                         if (index != "SequenceIndex") { //if you are not a header line, there will be a header line for each group if group file is given
 953                                 in >> name;             m->gobble(in);
 954                                 in >> DiffsToBestMatch; m->gobble(in);
 955                                 in >> BestMatchIndex; m->gobble(in);
 956                                 in >> BestMatchName; m->gobble(in);
 957                                 in >> DiffstToChimera; m->gobble(in);
 958                                 in >> IndexofLeftParent; m->gobble(in);
 959                                 in >> IndexOfRightParent; m->gobble(in);
 960                                 in >> parent1;  m->gobble(in);
 961                                 in >> parent2;  m->gobble(in);
 962                                 in >> temp1 >> temp2 >> temp3 >> temp4 >> temp5 >> temp6 >> temp7 >> temp8 >> flag; m->gobble(in);
 963
 964                                 //find unique name
 965                                 itUnique = uniqueNames.find(name);
 966
 967                                 if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; }
 968                                 else {
 969                                         name = itUnique->second;
 970                                         //is this name already in the file
 971                                         itNames = namesInFile.find((name));
 972
 973                                         if (itNames == namesInFile.end()) { //no not in file
 974                                                 if (flag == "good") { //are you really a no??
 975                                                         //is this sequence really not chimeric??
 976                                                         itChimeras = chimerasInFile.find(name);
 977
 978                                                         //then you really are a no so print, otherwise skip
 979                                                         if (itChimeras == chimerasInFile.end()) { print = true; }
 980                                                 }else{ print = true; }
 981                                         }
 982                                 }
 983
 984                                 if (print) {
 985                                         out << index << '\t' << name  << '\t' << DiffsToBestMatch << '\t' << BestMatchIndex << '\t';
 986                                         namesInFile.insert(name);
 987
 988                                         if (BestMatchName != "Null") {
 989                                                 itUnique = uniqueNames.find(BestMatchName);
 990                                                 if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find BestMatchName "+ BestMatchName + "."); m->mothurOutEndLine(); m->control_pressed = true; }
 991                                                 else {  out << itUnique->second << '\t';        }
 992                                         }else { out << "Null" << '\t'; }
 993
 994                                         out << DiffstToChimera << '\t' << IndexofLeftParent << '\t' << IndexOfRightParent << '\t';
 995
 996                                         if (parent1 != "Null") {
 997                                                 itUnique = uniqueNames.find(parent1);
 998                                                 if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parent1 "+ parent1 + "."); m->mothurOutEndLine(); m->control_pressed = true; }
 999                                                 else {  out << itUnique->second << '\t';        }
1000                                         }else { out << "Null" << '\t'; }
1001
1002                                         if (parent1 != "Null") {
1003                                                 itUnique = uniqueNames.find(parent2);
1004                                                 if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parent2 "+ parent2 + "."); m->mothurOutEndLine(); m->control_pressed = true; }
1005                                                 else {  out << itUnique->second << '\t';        }
1006                                         }else { out << "Null" << '\t'; }
1007
1008                                         out << temp1 << '\t' << temp2 << '\t' << temp3 << '\t' << temp4 << '\t' << temp5 << '\t' << temp6 << '\t' << temp7 << '\t' << temp8 << '\t' << flag << endl;
1009                                 }
1010                         }else { index = m->getline(in); m->gobble(in); }
1011                 }
1012                 in.close();
1013                 out.close();
1014
1015                 m->mothurRemove(outputFileName);
1016                 rename((outputFileName+".temp").c_str(), outputFileName.c_str());
1017
1018                 return total;
1019         }
1020         catch(exception& e) {
1021                 m->errorOut(e, "ChimeraPerseusCommand", "deconvoluteResults");
1022                 exit(1);
1023         }
1024 }
1025 //**********************************************************************************************************************
1026
1027