chimerapintailcommand.cpp

   1 /*
   2  *  chimerapintailcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 4/1/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "chimerapintailcommand.h"
  11 #include "pintail.h"
  12
  13 //***************************************************************************************************************
  14
  15 ChimeraPintailCommand::ChimeraPintailCommand(string option)  {
  16         try {
  17                 abort = false;
  18
  19                 //allow user to run help
  20                 if(option == "help") { help(); abort = true; }
  21
  22                 else {
  23                         //valid paramters for this command
  24                         string Array[] =  {"fasta","filter","processors","window" ,"increment","template","conservation","quantile","mask","outputdir","inputdir"};
  25                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  26
  27                         OptionParser parser(option);
  28                         map<string,string> parameters = parser.getParameters();
  29
  30                         ValidParameters validParameter("chimera.pintail");
  31                         map<string,string>::iterator it;
  32
  33                         //check to make sure all parameters are valid for command
  34                         for (it = parameters.begin(); it != parameters.end(); it++) {
  35                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  36                         }
  37
  38                         //if the user changes the input directory command factory will send this info to us in the output parameter
  39                         inputDir = validParameter.validFile(parameters, "inputdir", false);
  40                         if (inputDir == "not found"){   inputDir = "";          }
  41                         else {
  42                                 string path;
  43                                 it = parameters.find("template");
  44                                 //user has given a template file
  45                                 if(it != parameters.end()){
  46                                         path = m->hasPath(it->second);
  47                                         //if the user has not given a path then, add inputdir. else leave path alone.
  48                                         if (path == "") {       parameters["template"] = inputDir + it->second;         }
  49                                 }
  50
  51                                 it = parameters.find("conservation");
  52                                 //user has given a template file
  53                                 if(it != parameters.end()){
  54                                         path = m->hasPath(it->second);
  55                                         //if the user has not given a path then, add inputdir. else leave path alone.
  56                                         if (path == "") {       parameters["conservation"] = inputDir + it->second;             }
  57                                 }
  58
  59                                 it = parameters.find("quantile");
  60                                 //user has given a template file
  61                                 if(it != parameters.end()){
  62                                         path = m->hasPath(it->second);
  63                                         //if the user has not given a path then, add inputdir. else leave path alone.
  64                                         if (path == "") {       parameters["quantile"] = inputDir + it->second;         }
  65                                 }
  66                         }
  67
  68
  69                         //check for required parameters
  70                         fastafile = validParameter.validFile(parameters, "fasta", false);
  71                         if (fastafile == "not found") { fastafile = ""; m->mothurOut("fasta is a required parameter for the chimera.pintail command."); m->mothurOutEndLine(); abort = true;  }
  72                         else {
  73                                 m->splitAtDash(fastafile, fastaFileNames);
  74
  75                                 //go through files and make sure they are good, if not, then disregard them
  76                                 for (int i = 0; i < fastaFileNames.size(); i++) {
  77                                         if (inputDir != "") {
  78                                                 string path = m->hasPath(fastaFileNames[i]);
  79                                                 //if the user has not given a path then, add inputdir. else leave path alone.
  80                                                 if (path == "") {       fastaFileNames[i] = inputDir + fastaFileNames[i];               }
  81                                         }
  82
  83                                         int ableToOpen;
  84                                         ifstream in;
  85
  86                                         ableToOpen = m->openInputFile(fastaFileNames[i], in, "noerror");
  87
  88                                         //if you can't open it, try default location
  89                                         if (ableToOpen == 1) {
  90                                                 if (m->getDefaultPath() != "") { //default path is set
  91                                                         string tryPath = m->getDefaultPath() + m->getSimpleName(fastaFileNames[i]);
  92                                                         m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
  93                                                         ableToOpen = m->openInputFile(tryPath, in, "noerror");
  94                                                         fastaFileNames[i] = tryPath;
  95                                                 }
  96                                         }
  97                                         in.close();
  98
  99                                         if (ableToOpen == 1) {
 100                                                 m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
 101                                                 //erase from file list
 102                                                 fastaFileNames.erase(fastaFileNames.begin()+i);
 103                                                 i--;
 104                                         }
 105                                 }
 106
 107                                 //make sure there is at least one valid file left
 108                                 if (fastaFileNames.size() == 0) { m->mothurOut("no valid files."); m->mothurOutEndLine(); abort = true; }
 109                         }
 110
 111                         string temp;
 112                         temp = validParameter.validFile(parameters, "filter", false);                   if (temp == "not found") { temp = "F"; }
 113                         filter = m->isTrue(temp);
 114
 115                         temp = validParameter.validFile(parameters, "processors", false);               if (temp == "not found") { temp = "1"; }
 116                         convert(temp, processors);
 117
 118                         temp = validParameter.validFile(parameters, "window", false);                   if (temp == "not found") { temp = "0"; }
 119                         convert(temp, window);
 120
 121                         temp = validParameter.validFile(parameters, "increment", false);                if (temp == "not found") { temp = "25"; }
 122                         convert(temp, increment);
 123
 124                         maskfile = validParameter.validFile(parameters, "mask", false);
 125                         if (maskfile == "not found") { maskfile = "";  }
 126                         else if (maskfile != "default")  {
 127                                 if (inputDir != "") {
 128                                         string path = m->hasPath(maskfile);
 129                                         //if the user has not given a path then, add inputdir. else leave path alone.
 130                                         if (path == "") {       maskfile = inputDir + maskfile;         }
 131                                 }
 132
 133                                 ifstream in;
 134                                 int     ableToOpen = m->openInputFile(maskfile, in, "no error");
 135                                 if (ableToOpen == 1) {
 136                                         if (m->getDefaultPath() != "") { //default path is set
 137                                                         string tryPath = m->getDefaultPath() + m->getSimpleName(maskfile);
 138                                                         m->mothurOut("Unable to open " + maskfile + ". Trying default " + tryPath); m->mothurOutEndLine();
 139                                                         ableToOpen = m->openInputFile(tryPath, in, "noerror");
 140                                                         maskfile = tryPath;
 141                                         }
 142                                 }
 143                                         in.close();
 144
 145                                 if (ableToOpen == 1) {
 146                                                 m->mothurOut("Unable to open " + maskfile + "."); m->mothurOutEndLine();
 147                                                 abort = true;
 148                                 }
 149                         }
 150
 151
 152                         //if the user changes the output directory command factory will send this info to us in the output parameter
 153                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 154                                 outputDir = "";
 155                                 outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it
 156                         }
 157
 158                         templatefile = validParameter.validFile(parameters, "template", true);
 159                         if (templatefile == "not open") { abort = true; }
 160                         else if (templatefile == "not found") { templatefile = "";  m->mothurOut("template is a required parameter for the chimera.pintail command."); m->mothurOutEndLine(); abort = true;  }
 161
 162                         consfile = validParameter.validFile(parameters, "conservation", true);
 163                         if (consfile == "not open") { abort = true; }
 164                         else if (consfile == "not found") {
 165                                 consfile = "";
 166                                 //check for consfile
 167                                 string tempConsFile = m->getRootName(inputDir + m->getSimpleName(templatefile)) + "freq";
 168                                 ifstream FileTest(tempConsFile.c_str());
 169                                 if(FileTest){
 170                                         bool GoodFile = m->checkReleaseVersion(FileTest, m->getVersion());
 171                                         if (GoodFile) {
 172                                                 m->mothurOut("I found " + tempConsFile + " in your input file directory. I will use it to save time."); m->mothurOutEndLine();  consfile = tempConsFile;  FileTest.close();
 173                                         }
 174                                 }else {
 175                                         string tempConsFile = m->getDefaultPath() + m->getRootName(m->getSimpleName(templatefile)) + "freq";
 176                                         ifstream FileTest2(tempConsFile.c_str());
 177                                         if(FileTest2){
 178                                                 bool GoodFile = m->checkReleaseVersion(FileTest2, m->getVersion());
 179                                                 if (GoodFile) {
 180                                                         m->mothurOut("I found " + tempConsFile + " in your input file directory. I will use it to save time."); m->mothurOutEndLine();  consfile = tempConsFile;  FileTest2.close();
 181                                                 }
 182                                         }
 183                                 }
 184                         }
 185
 186                         quanfile = validParameter.validFile(parameters, "quantile", true);
 187                         if (quanfile == "not open") { abort = true; }
 188                         else if (quanfile == "not found") { quanfile = ""; }
 189                 }
 190         }
 191         catch(exception& e) {
 192                 m->errorOut(e, "ChimeraPintailCommand", "ChimeraPintailCommand");
 193                 exit(1);
 194         }
 195 }
 196 //**********************************************************************************************************************
 197
 198 void ChimeraPintailCommand::help(){
 199         try {
 200
 201                 m->mothurOut("The chimera.pintail command reads a fastafile and templatefile and outputs potentially chimeric sequences.\n");
 202                 m->mothurOut("This command was created using the algorythms described in the 'At Least 1 in 20 16S rRNA Sequence Records Currently Held in the Public Repositories is Estimated To Contain Substantial Anomalies' paper by Kevin E. Ashelford 1, Nadia A. Chuzhanova 3, John C. Fry 1, Antonia J. Jones 2 and Andrew J. Weightman 1.\n");
 203                 m->mothurOut("The chimera.pintail command parameters are fasta, template, filter, mask, processors, window, increment, conservation and quantile.\n");
 204                 m->mothurOut("The fasta parameter allows you to enter the fasta file containing your potentially chimeric sequences, and is required. \n");
 205                 m->mothurOut("You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amzon.fasta \n");
 206                 m->mothurOut("The template parameter allows you to enter a template file containing known non-chimeric sequences, and is required. \n");
 207                 m->mothurOut("The filter parameter allows you to specify if you would like to apply a vertical and 50% soft filter. \n");
 208                 m->mothurOut("The mask parameter allows you to specify a file containing one sequence you wish to use as a mask for the your sequences, by default no mask is applied.  You can apply an ecoli mask by typing, mask=default. \n");
 209                 m->mothurOut("The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n");
 210                 #ifdef USE_MPI
 211                 m->mothurOut("When using MPI, the processors parameter is set to the number of MPI processes running. \n");
 212                 #endif
 213                 m->mothurOut("The window parameter allows you to specify the window size for searching for chimeras, default=300. \n");
 214                 m->mothurOut("The increment parameter allows you to specify how far you move each window while finding chimeric sequences, default=25.\n");
 215                 m->mothurOut("The conservation parameter allows you to enter a frequency file containing the highest bases frequency at each place in the alignment.\n");
 216                 m->mothurOut("The quantile parameter allows you to enter a file containing quantiles for a template files sequences, if you use the filter the quantile file generated becomes unique to the fasta file you used.\n");
 217                 m->mothurOut("The chimera.pintail command should be in the following format: \n");
 218                 m->mothurOut("chimera.pintail(fasta=yourFastaFile, template=yourTemplate) \n");
 219                 m->mothurOut("Example: chimera.pintail(fasta=AD.align, template=silva.bacteria.fasta) \n");
 220                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n\n");
 221         }
 222         catch(exception& e) {
 223                 m->errorOut(e, "ChimeraPintailCommand", "help");
 224                 exit(1);
 225         }
 226 }
 227
 228 //***************************************************************************************************************
 229
 230 ChimeraPintailCommand::~ChimeraPintailCommand(){        /*      do nothing      */      }
 231
 232 //***************************************************************************************************************
 233
 234 int ChimeraPintailCommand::execute(){
 235         try{
 236
 237                 if (abort == true) { return 0; }
 238
 239                 for (int s = 0; s < fastaFileNames.size(); s++) {
 240
 241                         m->mothurOut("Checking sequences from " + fastaFileNames[s] + " ..." ); m->mothurOutEndLine();
 242
 243                         int start = time(NULL);
 244
 245                         //set user options
 246                         if (maskfile == "default") { m->mothurOut("I am using the default 236627 EU009184.1 Shigella dysenteriae str. FBD013."); m->mothurOutEndLine();  }
 247
 248                         //check for quantile to save the time
 249                         string tempQuan = "";
 250                         if ((!filter) && (maskfile == "")) {
 251                                 tempQuan = inputDir + m->getRootName(m->getSimpleName(templatefile)) + "pintail.quan";
 252                         }else if ((!filter) && (maskfile != "")) {
 253                                 tempQuan = inputDir + m->getRootName(m->getSimpleName(templatefile)) + "pintail.masked.quan";
 254                         }else if ((filter) && (maskfile != "")) {
 255                                 tempQuan = inputDir + m->getRootName(m->getSimpleName(templatefile)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastaFileNames[s])) + "masked.quan";
 256                         }else if ((filter) && (maskfile == "")) {
 257                                 tempQuan = inputDir + m->getRootName(m->getSimpleName(templatefile)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastaFileNames[s])) + "quan";
 258                         }
 259
 260                         ifstream FileTest(tempQuan.c_str());
 261                         if(FileTest){
 262                                 bool GoodFile = m->checkReleaseVersion(FileTest, m->getVersion());
 263                                 if (GoodFile) {
 264                                         m->mothurOut("I found " + tempQuan + " in your input file directory. I will use it to save time."); m->mothurOutEndLine();  quanfile = tempQuan;  FileTest.close();
 265                                 }
 266                         }else {
 267                                 string tryPath = m->getDefaultPath();
 268                                 string tempQuan = "";
 269                                 if ((!filter) && (maskfile == "")) {
 270                                         tempQuan = tryPath + m->getRootName(m->getSimpleName(templatefile)) + "pintail.quan";
 271                                 }else if ((!filter) && (maskfile != "")) {
 272                                         tempQuan = tryPath + m->getRootName(m->getSimpleName(templatefile)) + "pintail.masked.quan";
 273                                 }else if ((filter) && (maskfile != "")) {
 274                                         tempQuan = tryPath + m->getRootName(m->getSimpleName(templatefile)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastaFileNames[s])) + "masked.quan";
 275                                 }else if ((filter) && (maskfile == "")) {
 276                                         tempQuan = tryPath + m->getRootName(m->getSimpleName(templatefile)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastaFileNames[s])) + "quan";
 277                                 }
 278
 279                                 ifstream FileTest2(tempQuan.c_str());
 280                                 if(FileTest2){
 281                                         bool GoodFile = m->checkReleaseVersion(FileTest2, m->getVersion());
 282                                         if (GoodFile) {
 283                                                 m->mothurOut("I found " + tempQuan + " in your input file directory. I will use it to save time."); m->mothurOutEndLine();  quanfile = tempQuan;  FileTest2.close();
 284                                         }
 285                                 }
 286                         }
 287
 288                         chimera = new Pintail(fastaFileNames[s], templatefile, filter, processors, maskfile, consfile, quanfile, window, increment, outputDir);
 289
 290                         string outputFileName, accnosFileName;
 291                         if (maskfile != "") {
 292                                 outputFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + m->getSimpleName(m->getRootName(maskfile)) + ".pintail.chimeras";
 293                                 accnosFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + m->getSimpleName(m->getRootName(maskfile)) + ".pintail.accnos";
 294                         }else {
 295                                 outputFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s]))  + "pintail.chimeras";
 296                                 accnosFileName = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s]))  + "pintail.accnos";
 297                         }
 298
 299                         if (m->control_pressed) { delete chimera; for (int j = 0; j < outputNames.size(); j++) {        remove(outputNames[j].c_str()); }  return 0;    }
 300
 301                         if (chimera->getUnaligned()) {
 302                                 m->mothurOut("Your template sequences are different lengths, please correct."); m->mothurOutEndLine();
 303                                 delete chimera;
 304                                 return 0;
 305                         }
 306                         templateSeqsLength = chimera->getLength();
 307
 308                 #ifdef USE_MPI
 309                         int pid, end, numSeqsPerProcessor;
 310                                 int tag = 2001;
 311                                 vector<unsigned long int> MPIPos;
 312
 313                                 MPI_Status status;
 314                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 315                                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
 316
 317                                 MPI_File inMPI;
 318                                 MPI_File outMPI;
 319                                 MPI_File outMPIAccnos;
 320
 321                                 int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
 322                                 int inMode=MPI_MODE_RDONLY;
 323
 324                                 char outFilename[1024];
 325                                 strcpy(outFilename, outputFileName.c_str());
 326
 327                                 char outAccnosFilename[1024];
 328                                 strcpy(outAccnosFilename, accnosFileName.c_str());
 329
 330                                 char inFileName[1024];
 331                                 strcpy(inFileName, fastaFileNames[s].c_str());
 332
 333                                 MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 334                                 MPI_File_open(MPI_COMM_WORLD, outFilename, outMode, MPI_INFO_NULL, &outMPI);
 335                                 MPI_File_open(MPI_COMM_WORLD, outAccnosFilename, outMode, MPI_INFO_NULL, &outMPIAccnos);
 336
 337                                 if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPI);   MPI_File_close(&outMPIAccnos);  for (int j = 0; j < outputNames.size(); j++) {   remove(outputNames[j].c_str()); }  delete chimera; return 0;  }
 338
 339                                 if (pid == 0) { //you are the root process
 340
 341                                         MPIPos = m->setFilePosFasta(fastaFileNames[s], numSeqs); //fills MPIPos, returns numSeqs
 342
 343                                         //send file positions to all processes
 344                                         for(int i = 1; i < processors; i++) {
 345                                                 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 346                                                 MPI_Send(&MPIPos[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
 347                                         }
 348
 349                                         //figure out how many sequences you have to align
 350                                         numSeqsPerProcessor = numSeqs / processors;
 351                                         int startIndex =  pid * numSeqsPerProcessor;
 352                                         if(pid == (processors - 1)){    numSeqsPerProcessor = numSeqs - pid * numSeqsPerProcessor;      }
 353
 354                                         //do your part
 355                                         driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPI, outMPIAccnos, MPIPos);
 356
 357                                         if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPI);   MPI_File_close(&outMPIAccnos);  remove(outputFileName.c_str());  remove(accnosFileName.c_str());  for (int j = 0; j < outputNames.size(); j++) { remove(outputNames[j].c_str()); }  delete chimera; return 0;  }
 358
 359                                 }else{ //you are a child process
 360                                         MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 361                                         MPIPos.resize(numSeqs+1);
 362                                         MPI_Recv(&MPIPos[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
 363
 364                                         //figure out how many sequences you have to align
 365                                         numSeqsPerProcessor = numSeqs / processors;
 366                                         int startIndex =  pid * numSeqsPerProcessor;
 367                                         if(pid == (processors - 1)){    numSeqsPerProcessor = numSeqs - pid * numSeqsPerProcessor;      }
 368
 369                                         //do your part
 370                                         driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPI, outMPIAccnos, MPIPos);
 371
 372                                         if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPI);   MPI_File_close(&outMPIAccnos);  for (int j = 0; j < outputNames.size(); j++) {   remove(outputNames[j].c_str()); }  delete chimera; return 0;  }
 373                                 }
 374
 375                                 //close files
 376                                 MPI_File_close(&inMPI);
 377                                 MPI_File_close(&outMPI);
 378                                 MPI_File_close(&outMPIAccnos);
 379                                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 380                 #else
 381                         vector<unsigned long int> positions = m->divideFile(fastaFileNames[s], processors);
 382
 383                         for (int i = 0; i < (positions.size()-1); i++) {
 384                                 lines.push_back(new linePair(positions[i], positions[(i+1)]));
 385                         }
 386
 387                         //break up file
 388                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 389                                 if(processors == 1){
 390
 391                                         numSeqs = driver(lines[0], outputFileName, fastaFileNames[s], accnosFileName);
 392
 393                                         if (m->control_pressed) { remove(outputFileName.c_str()); remove(accnosFileName.c_str()); for (int j = 0; j < outputNames.size(); j++) {        remove(outputNames[j].c_str()); } for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear(); delete chimera; return 0; }
 394
 395                                 }else{
 396                                         processIDS.resize(0);
 397
 398                                         numSeqs = createProcesses(outputFileName, fastaFileNames[s], accnosFileName);
 399
 400                                         rename((outputFileName + toString(processIDS[0]) + ".temp").c_str(), outputFileName.c_str());
 401                                         rename((accnosFileName + toString(processIDS[0]) + ".temp").c_str(), accnosFileName.c_str());
 402
 403                                         //append output files
 404                                         for(int i=1;i<processors;i++){
 405                                                 m->appendFiles((outputFileName + toString(processIDS[i]) + ".temp"), outputFileName);
 406                                                 remove((outputFileName + toString(processIDS[i]) + ".temp").c_str());
 407                                         }
 408
 409                                         //append output files
 410                                         for(int i=1;i<processors;i++){
 411                                                 m->appendFiles((accnosFileName + toString(processIDS[i]) + ".temp"), accnosFileName);
 412                                                 remove((accnosFileName + toString(processIDS[i]) + ".temp").c_str());
 413                                         }
 414
 415                                         if (m->control_pressed) {
 416                                                 remove(outputFileName.c_str());
 417                                                 remove(accnosFileName.c_str());
 418                                                 for (int j = 0; j < outputNames.size(); j++) {  remove(outputNames[j].c_str()); }
 419                                                 for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 420                                                 delete chimera;
 421                                                 return 0;
 422                                         }
 423                                 }
 424
 425                         #else
 426                                 numSeqs = driver(lines[0], outputFileName, fastaFileNames[s], accnosFileName);
 427
 428                                 if (m->control_pressed) { remove(outputFileName.c_str()); remove(accnosFileName.c_str()); for (int j = 0; j < outputNames.size(); j++) {        remove(outputNames[j].c_str()); } for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear(); delete chimera; return 0; }
 429                         #endif
 430
 431                 #endif
 432
 433                         delete chimera;
 434                         for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 435
 436                         outputNames.push_back(outputFileName);
 437                         outputNames.push_back(accnosFileName);
 438
 439                         m->mothurOutEndLine();
 440                         m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences."); m->mothurOutEndLine();
 441                 }
 442
 443                 m->mothurOutEndLine();
 444                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 445                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 446                 m->mothurOutEndLine();
 447
 448                 return 0;
 449
 450         }
 451         catch(exception& e) {
 452                 m->errorOut(e, "ChimeraPintailCommand", "execute");
 453                 exit(1);
 454         }
 455 }
 456 //**********************************************************************************************************************
 457
 458 int ChimeraPintailCommand::driver(linePair* filePos, string outputFName, string filename, string accnos){
 459         try {
 460                 ofstream out;
 461                 m->openOutputFile(outputFName, out);
 462
 463                 ofstream out2;
 464                 m->openOutputFile(accnos, out2);
 465
 466                 ifstream inFASTA;
 467                 m->openInputFile(filename, inFASTA);
 468
 469                 inFASTA.seekg(filePos->start);
 470
 471                 bool done = false;
 472                 int count = 0;
 473
 474                 while (!done) {
 475
 476                         if (m->control_pressed) {       return 1;       }
 477
 478                         Sequence* candidateSeq = new Sequence(inFASTA);  m->gobble(inFASTA);
 479
 480                         if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
 481
 482                                 if (candidateSeq->getAligned().length() != templateSeqsLength)  {  //chimeracheck does not require seqs to be aligned
 483                                         m->mothurOut(candidateSeq->getName() + " is not the same length as the template sequences. Skipping."); m->mothurOutEndLine();
 484                                 }else{
 485                                         //find chimeras
 486                                         chimera->getChimeras(candidateSeq);
 487
 488                                         if (m->control_pressed) {       delete candidateSeq; return 1;  }
 489
 490                                         //print results
 491                                         chimera->print(out, out2);
 492                                 }
 493                                 count++;
 494                         }
 495                         delete candidateSeq;
 496
 497                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 498                                 unsigned long int pos = inFASTA.tellg();
 499                                 if ((pos == -1) || (pos >= filePos->end)) { break; }
 500                         #else
 501                                 if (inFASTA.eof()) { break; }
 502                         #endif
 503
 504                         //report progress
 505                         if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine();         }
 506                 }
 507                 //report progress
 508                 if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine();         }
 509
 510                 out.close();
 511                 out2.close();
 512                 inFASTA.close();
 513
 514                 return count;
 515         }
 516         catch(exception& e) {
 517                 m->errorOut(e, "ChimeraPintailCommand", "driver");
 518                 exit(1);
 519         }
 520 }
 521 //**********************************************************************************************************************
 522 #ifdef USE_MPI
 523 int ChimeraPintailCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& outMPI, MPI_File& outAccMPI, vector<unsigned long int>& MPIPos){
 524         try {
 525
 526                 MPI_Status status;
 527                 int pid;
 528                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 529
 530                 for(int i=0;i<num;i++){
 531
 532                         if (m->control_pressed) {       return 1;       }
 533
 534                         //read next sequence
 535                         int length = MPIPos[start+i+1] - MPIPos[start+i];
 536
 537                         char* buf4 = new char[length];
 538                         MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
 539
 540                         string tempBuf = buf4;
 541                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
 542                         istringstream iss (tempBuf,istringstream::in);
 543                         delete buf4;
 544
 545                         Sequence* candidateSeq = new Sequence(iss);  m->gobble(iss);
 546
 547                         if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
 548
 549                                 if      (candidateSeq->getAligned().length() != templateSeqsLength) {  //chimeracheck does not require seqs to be aligned
 550                                         m->mothurOut(candidateSeq->getName() + " is not the same length as the template sequences. Skipping."); m->mothurOutEndLine();
 551                                 }else{
 552                                         //find chimeras
 553                                         chimera->getChimeras(candidateSeq);
 554
 555                                         if (m->control_pressed) {       delete candidateSeq; return 1;  }
 556
 557                                         //print results
 558                                         bool isChimeric = chimera->print(outMPI, outAccMPI);
 559                                 }
 560                         }
 561                         delete candidateSeq;
 562
 563                         //report progress
 564                         if((i+1) % 100 == 0){  cout << "Processing sequence: " << (i+1) << endl;        m->mothurOutJustToLog("Processing sequence: " + toString(i+1) + "\n");          }
 565                 }
 566                 //report progress
 567                 if(num % 100 != 0){             cout << "Processing sequence: " << num << endl; m->mothurOutJustToLog("Processing sequence: " + toString(num) + "\n");  }
 568
 569
 570                 return 0;
 571         }
 572         catch(exception& e) {
 573                 m->errorOut(e, "ChimeraPintailCommand", "driverMPI");
 574                 exit(1);
 575         }
 576 }
 577 #endif
 578
 579 /**************************************************************************************************/
 580
 581 int ChimeraPintailCommand::createProcesses(string outputFileName, string filename, string accnos) {
 582         try {
 583 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 584                 int process = 0;
 585                 int num = 0;
 586
 587                 //loop through and create all the processes you want
 588                 while (process != processors) {
 589                         int pid = fork();
 590
 591                         if (pid > 0) {
 592                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 593                                 process++;
 594                         }else if (pid == 0){
 595                                 num = driver(lines[process], outputFileName + toString(getpid()) + ".temp", filename, accnos + toString(getpid()) + ".temp");
 596
 597                                 //pass numSeqs to parent
 598                                 ofstream out;
 599                                 string tempFile = outputFileName + toString(getpid()) + ".num.temp";
 600                                 m->openOutputFile(tempFile, out);
 601                                 out << num << endl;
 602                                 out.close();
 603
 604                                 exit(0);
 605                         }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
 606                 }
 607
 608                 //force parent to wait until all the processes are done
 609                 for (int i=0;i<processors;i++) {
 610                         int temp = processIDS[i];
 611                         wait(&temp);
 612                 }
 613
 614                 for (int i = 0; i < processIDS.size(); i++) {
 615                         ifstream in;
 616                         string tempFile =  outputFileName + toString(processIDS[i]) + ".num.temp";
 617                         m->openInputFile(tempFile, in);
 618                         if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
 619                         in.close(); remove(tempFile.c_str());
 620                 }
 621
 622                 return num;
 623 #endif
 624         }
 625         catch(exception& e) {
 626                 m->errorOut(e, "ChimeraPintailCommand", "createProcesses");
 627                 exit(1);
 628         }
 629 }
 630
 631 /**************************************************************************************************/
 632
 633