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