]> git.donarmstrong.com Git - mothur.git/blob - screenseqscommand.cpp
working on chimera change to add trim feature, fixed bug in print of distance file...
[mothur.git] / screenseqscommand.cpp
1 /*
2  *  screenseqscommand.cpp
3  *  Mothur
4  *
5  *  Created by Pat Schloss on 6/3/09.
6  *  Copyright 2009 Patrick D. Schloss. All rights reserved.
7  *
8  */
9
10 #include "screenseqscommand.h"
11 #include "sequence.hpp"
12
13 //**********************************************************************************************************************
14 vector<string> ScreenSeqsCommand::getValidParameters(){ 
15         try {
16                 string Array[] =  {"fasta", "start", "end", "maxambig", "maxhomop","optimize","criteria", "minlength", "maxlength",
17                                                                         "name", "group", "alignreport","processors","outputdir","inputdir"};
18                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
19                 return myArray;
20         }
21         catch(exception& e) {
22                 m->errorOut(e, "ScreenSeqsCommand", "getValidParameters");
23                 exit(1);
24         }
25 }
26 //**********************************************************************************************************************
27 ScreenSeqsCommand::ScreenSeqsCommand(){ 
28         try {
29                 abort = true;
30                 //initialize outputTypes
31                 vector<string> tempOutNames;
32                 outputTypes["fasta"] = tempOutNames;
33                 outputTypes["name"] = tempOutNames;
34                 outputTypes["group"] = tempOutNames;
35                 outputTypes["alignreport"] = tempOutNames;
36                 outputTypes["accnos"] = tempOutNames;
37         }
38         catch(exception& e) {
39                 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
40                 exit(1);
41         }
42 }
43 //**********************************************************************************************************************
44 vector<string> ScreenSeqsCommand::getRequiredParameters(){      
45         try {
46                 string Array[] =  {"fasta"};
47                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
48                 return myArray;
49         }
50         catch(exception& e) {
51                 m->errorOut(e, "ScreenSeqsCommand", "getRequiredParameters");
52                 exit(1);
53         }
54 }
55 //**********************************************************************************************************************
56 vector<string> ScreenSeqsCommand::getRequiredFiles(){   
57         try {
58                 vector<string> myArray;
59                 return myArray;
60         }
61         catch(exception& e) {
62                 m->errorOut(e, "ScreenSeqsCommand", "getRequiredFiles");
63                 exit(1);
64         }
65 }
66 //***************************************************************************************************************
67
68 ScreenSeqsCommand::ScreenSeqsCommand(string option)  {
69         try {
70                 abort = false;
71                 
72                 //allow user to run help
73                 if(option == "help") { help(); abort = true; }
74                 
75                 else {
76                         //valid paramters for this command
77                         string AlignArray[] =  {"fasta", "start", "end", "maxambig", "maxhomop","optimize","criteria", "minlength", "maxlength",
78                                                                         "name", "group", "alignreport","processors","outputdir","inputdir"};
79                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
80                         
81                         OptionParser parser(option);
82                         map<string,string> parameters = parser.getParameters();
83                         
84                         ValidParameters validParameter("screen.seqs");
85                         map<string,string>::iterator it;
86                         
87                         //check to make sure all parameters are valid for command
88                         for (it = parameters.begin(); it != parameters.end(); it++) { 
89                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
90                         }
91                         
92                         //initialize outputTypes
93                         vector<string> tempOutNames;
94                         outputTypes["fasta"] = tempOutNames;
95                         outputTypes["name"] = tempOutNames;
96                         outputTypes["group"] = tempOutNames;
97                         outputTypes["alignreport"] = tempOutNames;
98                         outputTypes["accnos"] = tempOutNames;
99                         
100                         //if the user changes the input directory command factory will send this info to us in the output parameter 
101                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
102                         if (inputDir == "not found"){   inputDir = "";          }
103                         else {
104                                 string path;
105                                 it = parameters.find("fasta");
106                                 //user has given a template file
107                                 if(it != parameters.end()){ 
108                                         path = m->hasPath(it->second);
109                                         //if the user has not given a path then, add inputdir. else leave path alone.
110                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
111                                 }
112                                 
113                                 it = parameters.find("group");
114                                 //user has given a template file
115                                 if(it != parameters.end()){ 
116                                         path = m->hasPath(it->second);
117                                         //if the user has not given a path then, add inputdir. else leave path alone.
118                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
119                                 }
120                                 
121                                 it = parameters.find("name");
122                                 //user has given a template file
123                                 if(it != parameters.end()){ 
124                                         path = m->hasPath(it->second);
125                                         //if the user has not given a path then, add inputdir. else leave path alone.
126                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
127                                 }
128                                 
129                                 it = parameters.find("alignreport");
130                                 //user has given a template file
131                                 if(it != parameters.end()){ 
132                                         path = m->hasPath(it->second);
133                                         //if the user has not given a path then, add inputdir. else leave path alone.
134                                         if (path == "") {       parameters["alignreport"] = inputDir + it->second;              }
135                                 }
136                         }
137
138                         //check for required parameters
139                         fastafile = validParameter.validFile(parameters, "fasta", true);
140                         if (fastafile == "not found") { m->mothurOut("fasta is a required parameter for the screen.seqs command."); m->mothurOutEndLine(); abort = true; }
141                         else if (fastafile == "not open") { abort = true; }     
142         
143                         groupfile = validParameter.validFile(parameters, "group", true);
144                         if (groupfile == "not open") { abort = true; }  
145                         else if (groupfile == "not found") { groupfile = ""; }
146                         
147                         namefile = validParameter.validFile(parameters, "name", true);
148                         if (namefile == "not open") { abort = true; }
149                         else if (namefile == "not found") { namefile = ""; }    
150
151                         alignreport = validParameter.validFile(parameters, "alignreport", true);
152                         if (alignreport == "not open") { abort = true; }
153                         else if (alignreport == "not found") { alignreport = ""; }      
154                         
155                         //if the user changes the output directory command factory will send this info to us in the output parameter 
156                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  
157                                 outputDir = ""; 
158                                 outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it       
159                         }
160
161                         //check for optional parameter and set defaults
162                         // ...at some point should added some additional type checking...
163                         string temp;
164                         temp = validParameter.validFile(parameters, "start", false);            if (temp == "not found") { temp = "-1"; }
165                         convert(temp, startPos); 
166                 
167                         temp = validParameter.validFile(parameters, "end", false);                      if (temp == "not found") { temp = "-1"; }
168                         convert(temp, endPos);  
169
170                         temp = validParameter.validFile(parameters, "maxambig", false);         if (temp == "not found") { temp = "-1"; }
171                         convert(temp, maxAmbig);  
172
173                         temp = validParameter.validFile(parameters, "maxhomop", false);         if (temp == "not found") { temp = "-1"; }
174                         convert(temp, maxHomoP);  
175
176                         temp = validParameter.validFile(parameters, "minlength", false);        if (temp == "not found") { temp = "-1"; }
177                         convert(temp, minLength); 
178                         
179                         temp = validParameter.validFile(parameters, "maxlength", false);        if (temp == "not found") { temp = "-1"; }
180                         convert(temp, maxLength); 
181                         
182                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
183                         convert(temp, processors); 
184                         
185                         temp = validParameter.validFile(parameters, "optimize", false); //optimizing trumps the optimized values original value
186                         if (temp == "not found"){       temp = "";              }
187                         else {  m->splitAtDash(temp, optimize);         }
188                         
189                         //check for invalid optimize options
190                         set<string> validOptimizers;
191                         validOptimizers.insert("start"); validOptimizers.insert("end"); validOptimizers.insert("maxambig"); validOptimizers.insert("maxhomop"); validOptimizers.insert("minlength"); validOptimizers.insert("maxlength");
192                         
193                         for (int i = 0; i < optimize.size(); i++) { 
194                                 if (validOptimizers.count(optimize[i]) == 0) { 
195                                         m->mothurOut(optimize[i] + " is not a valid optimizer. Valid options are start, end, maxambig, maxhomop, minlength and maxlength."); m->mothurOutEndLine();
196                                         optimize.erase(optimize.begin()+i);
197                                         i--;
198                                 }
199                         }
200                         
201                         temp = validParameter.validFile(parameters, "criteria", false); if (temp == "not found"){       temp = "90";                            }
202                         convert(temp, criteria); 
203                 }
204
205         }
206         catch(exception& e) {
207                 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
208                 exit(1);
209         }
210 }
211 //**********************************************************************************************************************
212
213 void ScreenSeqsCommand::help(){
214         try {
215                 m->mothurOut("The screen.seqs command reads a fastafile and creates .....\n");
216                 m->mothurOut("The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group, optimize, criteria and processors.\n");
217                 m->mothurOut("The fasta parameter is required.\n");
218                 m->mothurOut("The start parameter .... The default is -1.\n");
219                 m->mothurOut("The end parameter .... The default is -1.\n");
220                 m->mothurOut("The maxambig parameter .... The default is -1.\n");
221                 m->mothurOut("The maxhomop parameter .... The default is -1.\n");
222                 m->mothurOut("The minlength parameter .... The default is -1.\n");
223                 m->mothurOut("The maxlength parameter .... The default is -1.\n");
224                 m->mothurOut("The processors parameter allows you to specify the number of processors to use while running the command. The default is 1.\n");
225                 m->mothurOut("The optimize and criteria parameters allow you set the start, end, maxabig, maxhomop, minlength and maxlength parameters relative to your set of sequences .\n");
226                 m->mothurOut("For example optimize=start-end, criteria=90, would set the start and end values to the position 90% of your sequences started and ended.\n");
227                 m->mothurOut("The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n");
228                 m->mothurOut("The screen.seqs command should be in the following format: \n");
229                 m->mothurOut("screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig,  \n");
230                 m->mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength)  \n");    
231                 m->mothurOut("Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n");
232                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
233
234         }
235         catch(exception& e) {
236                 m->errorOut(e, "ScreenSeqsCommand", "help");
237                 exit(1);
238         }
239 }
240
241 //***************************************************************************************************************
242
243 ScreenSeqsCommand::~ScreenSeqsCommand(){        /*      do nothing      */      }
244
245 //***************************************************************************************************************
246
247 int ScreenSeqsCommand::execute(){
248         try{
249                 
250                 if (abort == true) { return 0; }
251                 
252                 //if the user want to optimize we need to no the 90% mark
253                 vector<unsigned long int> positions;
254                 if (optimize.size() != 0) {  getSummary(positions); } //get summary is paralellized so we need to divideFile, no need to do this step twice so I moved it here
255                 else { 
256                         positions = m->divideFile(fastafile, processors);
257                         for (int i = 0; i < (positions.size()-1); i++) {
258                                 lines.push_back(new linePair(positions[i], positions[(i+1)]));
259                         }       
260                 }
261                                 
262                 string goodSeqFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "good" + m->getExtension(fastafile);
263                 string badAccnosFile =  outputDir + m->getRootName(m->getSimpleName(fastafile)) + "bad.accnos";
264                 
265                 int numFastaSeqs = 0;
266                 set<string> badSeqNames;
267                 int start = time(NULL);
268                 
269 #ifdef USE_MPI  
270                         int pid, numSeqsPerProcessor; 
271                         int tag = 2001;
272                         vector<unsigned long int> MPIPos;
273                         
274                         MPI_Status status; 
275                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
276                         MPI_Comm_size(MPI_COMM_WORLD, &processors); 
277
278                         MPI_File inMPI;
279                         MPI_File outMPIGood;
280                         MPI_File outMPIBadAccnos;
281                         
282                         int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
283                         int inMode=MPI_MODE_RDONLY; 
284                         
285                         char outGoodFilename[1024];
286                         strcpy(outGoodFilename, goodSeqFile.c_str());
287
288                         char outBadAccnosFilename[1024];
289                         strcpy(outBadAccnosFilename, badAccnosFile.c_str());
290
291                         char inFileName[1024];
292                         strcpy(inFileName, fastafile.c_str());
293                         
294                         MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
295                         MPI_File_open(MPI_COMM_WORLD, outGoodFilename, outMode, MPI_INFO_NULL, &outMPIGood);
296                         MPI_File_open(MPI_COMM_WORLD, outBadAccnosFilename, outMode, MPI_INFO_NULL, &outMPIBadAccnos);
297                         
298                         if (m->control_pressed) { MPI_File_close(&inMPI);  MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
299                         
300                         if (pid == 0) { //you are the root process 
301                                 
302                                 MPIPos = m->setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
303                                 
304                                 //send file positions to all processes
305                                 for(int i = 1; i < processors; i++) { 
306                                         MPI_Send(&numFastaSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
307                                         MPI_Send(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
308                                 }
309                                 
310                                 //figure out how many sequences you have to align
311                                 numSeqsPerProcessor = numFastaSeqs / processors;
312                                 int startIndex =  pid * numSeqsPerProcessor;
313                                 if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
314                 //      cout << pid << '\t' << numSeqsPerProcessor << '\t' <<   startIndex << endl;
315                                 //align your part
316                                 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
317                 //cout << pid << " done" << endl;
318                                 if (m->control_pressed) { MPI_File_close(&inMPI);  MPI_File_close(&outMPIGood);  MPI_File_close(&outMPIBadAccnos);  return 0; }
319
320                                 for (int i = 1; i < processors; i++) {
321                                 
322                                         //get bad lists
323                                         int badSize;
324                                         MPI_Recv(&badSize, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
325                                         /*for (int j = 0; j < badSize; j++) {
326                                                 int length;
327                                                 MPI_Recv(&length, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);  //recv the length of the name
328                                                 char* buf2 = new char[length];                                                                          //make space to recieve it
329                                                 MPI_Recv(buf2, length, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status);  //get name
330                                                 
331                                                 string tempBuf = buf2;
332                                                 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
333                                                 delete buf2;
334                                                 
335                                                 badSeqNames.insert(tempBuf);
336                                         }*/
337                                 }
338                         }else{ //you are a child process
339                                 MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
340                                 MPIPos.resize(numFastaSeqs+1);
341                                 MPI_Recv(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
342
343                                 //figure out how many sequences you have to align
344                                 numSeqsPerProcessor = numFastaSeqs / processors;
345                                 int startIndex =  pid * numSeqsPerProcessor;
346                                 if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
347                 //cout << pid << '\t' << numSeqsPerProcessor << '\t' <<         startIndex << endl;             
348                                 //align your part
349                                 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
350 //cout << pid << " done" << endl;
351                                 if (m->control_pressed) { MPI_File_close(&inMPI);  MPI_File_close(&outMPIGood);  MPI_File_close(&outMPIBadAccnos); return 0; }
352                                 
353                                 //send bad list 
354                                 int badSize = badSeqNames.size();
355                                 MPI_Send(&badSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
356                                 
357                                 /*
358                                 set<string>::iterator it;
359                                 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
360                                         string name = *it;
361                                         int length = name.length();
362                                         char* buf2 = new char[length];
363                                         memcpy(buf2, name.c_str(), length);
364                                         
365                                         MPI_Send(&length, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
366                                         MPI_Send(buf2, length, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
367                                 }*/
368                         }
369                         
370                         //close files 
371                         MPI_File_close(&inMPI);
372                         MPI_File_close(&outMPIGood);
373                         MPI_File_close(&outMPIBadAccnos);
374                         MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
375                                         
376 #else
377                                                 
378         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
379                         if(processors == 1){
380                                 numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
381                                 
382                                 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
383                                 
384                         }else{
385                                 processIDS.resize(0);
386                                 
387                                 numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames); 
388                                 
389                                 rename((goodSeqFile + toString(processIDS[0]) + ".temp").c_str(), goodSeqFile.c_str());
390                                 rename((badAccnosFile + toString(processIDS[0]) + ".temp").c_str(), badAccnosFile.c_str());
391                                 
392                                 //append alignment and report files
393                                 for(int i=1;i<processors;i++){
394                                         m->appendFiles((goodSeqFile + toString(processIDS[i]) + ".temp"), goodSeqFile);
395                                         remove((goodSeqFile + toString(processIDS[i]) + ".temp").c_str());
396                         
397                                         m->appendFiles((badAccnosFile + toString(processIDS[i]) + ".temp"), badAccnosFile);
398                                         remove((badAccnosFile + toString(processIDS[i]) + ".temp").c_str());
399                                 }
400                                 
401                                 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
402                                 
403                                 //read badSeqs in because root process doesnt know what other "bad" seqs the children found
404                                 ifstream inBad;
405                                 int ableToOpen = m->openInputFile(badAccnosFile, inBad, "no error");
406                                 
407                                 if (ableToOpen == 0) {
408                                         badSeqNames.clear();
409                                         string tempName;
410                                         while (!inBad.eof()) {
411                                                 inBad >> tempName; m->gobble(inBad);
412                                                 badSeqNames.insert(tempName);
413                                         }
414                                         inBad.close();
415                                 }
416                         }
417         #else
418                         numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
419                         
420                         if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
421                         
422         #endif
423
424 #endif          
425
426                 #ifdef USE_MPI
427                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
428                                         
429                         if (pid == 0) { //only one process should fix files
430                         
431                                 //read accnos file with all names in it, process 0 just has its names
432                                 MPI_File inMPIAccnos;
433                                 MPI_Offset size;
434                         
435                                 char inFileName[1024];
436                                 strcpy(inFileName, badAccnosFile.c_str());
437                         
438                                 MPI_File_open(MPI_COMM_SELF, inFileName, inMode, MPI_INFO_NULL, &inMPIAccnos);  //comm, filename, mode, info, filepointer
439                                 MPI_File_get_size(inMPIAccnos, &size);
440                         
441                                 char* buffer = new char[size];
442                                 MPI_File_read(inMPIAccnos, buffer, size, MPI_CHAR, &status);
443                         
444                                 string tempBuf = buffer;
445                                 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size);  }
446                                 istringstream iss (tempBuf,istringstream::in);
447
448                                 delete buffer;
449                                 MPI_File_close(&inMPIAccnos);
450                                 
451                                 badSeqNames.clear();
452                                 string tempName;
453                                 while (!iss.eof()) {
454                                         iss >> tempName; m->gobble(iss);
455                                         badSeqNames.insert(tempName);
456                                 }
457                 #endif
458                                                                                                                                                                         
459                 if(namefile != "" && groupfile != "")   {       
460                         screenNameGroupFile(badSeqNames);       
461                         if (m->control_pressed) {  remove(goodSeqFile.c_str()); return 0; }
462                 }else if(namefile != "")        {       
463                         screenNameGroupFile(badSeqNames);
464                         if (m->control_pressed) {  remove(goodSeqFile.c_str());  return 0; }    
465                 }else if(groupfile != "")                               {       screenGroupFile(badSeqNames);           }       // this screens just the group
466                 
467                 if (m->control_pressed) { remove(goodSeqFile.c_str());  return 0; }
468
469                 if(alignreport != "")                                   {       screenAlignReport(badSeqNames);         }
470                 
471                 if (m->control_pressed) { remove(goodSeqFile.c_str());  return 0; }
472                 
473                 #ifdef USE_MPI
474                         }
475                 #endif
476
477                 m->mothurOutEndLine();
478                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
479                 m->mothurOut(goodSeqFile); m->mothurOutEndLine();       outputTypes["fasta"].push_back(goodSeqFile);
480                 m->mothurOut(badAccnosFile); m->mothurOutEndLine();      outputTypes["accnos"].push_back(badAccnosFile);
481                 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
482                 m->mothurOutEndLine();
483                 m->mothurOutEndLine();
484
485                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
486                 m->mothurOutEndLine();
487
488                 return 0;
489         }
490         catch(exception& e) {
491                 m->errorOut(e, "ScreenSeqsCommand", "execute");
492                 exit(1);
493         }
494 }
495
496 //***************************************************************************************************************
497
498 int ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
499         try {
500                 ifstream inputNames;
501                 m->openInputFile(namefile, inputNames);
502                 set<string> badSeqGroups;
503                 string seqName, seqList, group;
504                 set<string>::iterator it;
505
506                 string goodNameFile = outputDir + m->getRootName(m->getSimpleName(namefile)) + "good" + m->getExtension(namefile);
507                 outputNames.push_back(goodNameFile);  outputTypes["name"].push_back(goodNameFile);
508                 
509                 ofstream goodNameOut;   m->openOutputFile(goodNameFile, goodNameOut);
510                 
511                 while(!inputNames.eof()){
512                         if (m->control_pressed) { goodNameOut.close();  inputNames.close(); remove(goodNameFile.c_str());  return 0; }
513
514                         inputNames >> seqName >> seqList;
515                         it = badSeqNames.find(seqName);
516                                 
517                         if(it != badSeqNames.end()){
518                                 badSeqNames.erase(it);
519                                 
520                                 if(namefile != ""){
521                                         int start = 0;
522                                         for(int i=0;i<seqList.length();i++){
523                                                 if(seqList[i] == ','){
524                                                         badSeqGroups.insert(seqList.substr(start,i-start));
525                                                         start = i+1;
526                                                 }                                       
527                                         }
528                                         badSeqGroups.insert(seqList.substr(start,seqList.length()-start));
529                                 }
530                         }
531                         else{
532                                 goodNameOut << seqName << '\t' << seqList << endl;
533                         }
534                         m->gobble(inputNames);
535                 }
536                 inputNames.close();
537                 goodNameOut.close();
538         
539                 //we were unable to remove some of the bad sequences
540                 if (badSeqNames.size() != 0) {
541                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
542                                 m->mothurOut("Your namefile does not include the sequence " + *it + " please correct."); 
543                                 m->mothurOutEndLine();
544                         }
545                 }
546
547                 if(groupfile != ""){
548                         
549                         ifstream inputGroups;
550                         m->openInputFile(groupfile, inputGroups);
551
552                         string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
553                         outputNames.push_back(goodGroupFile);   outputTypes["group"].push_back(goodGroupFile);
554                         
555                         ofstream goodGroupOut;  m->openOutputFile(goodGroupFile, goodGroupOut);
556                         
557                         while(!inputGroups.eof()){
558                                 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodNameFile.c_str());  remove(goodGroupFile.c_str()); return 0; }
559
560                                 inputGroups >> seqName >> group;
561                                 
562                                 it = badSeqGroups.find(seqName);
563                                 
564                                 if(it != badSeqGroups.end()){
565                                         badSeqGroups.erase(it);
566                                 }
567                                 else{
568                                         goodGroupOut << seqName << '\t' << group << endl;
569                                 }
570                                 m->gobble(inputGroups);
571                         }
572                         inputGroups.close();
573                         goodGroupOut.close();
574                         
575                         //we were unable to remove some of the bad sequences
576                         if (badSeqGroups.size() != 0) {
577                                 for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {  
578                                         m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct."); 
579                                         m->mothurOutEndLine();
580                                 }
581                         }
582                 }
583                 
584                 
585                 return 0;
586         
587         }
588         catch(exception& e) {
589                 m->errorOut(e, "ScreenSeqsCommand", "screenNameGroupFile");
590                 exit(1);
591         }
592 }
593 //***************************************************************************************************************
594 int ScreenSeqsCommand::getSummary(vector<unsigned long int>& positions){
595         try {
596                 
597                 vector<int> startPosition;
598                 vector<int> endPosition;
599                 vector<int> seqLength;
600                 vector<int> ambigBases;
601                 vector<int> longHomoPolymer;
602                 
603                 vector<unsigned long int> positions = m->divideFile(fastafile, processors);
604                                 
605                 for (int i = 0; i < (positions.size()-1); i++) {
606                         lines.push_back(new linePair(positions[i], positions[(i+1)]));
607                 }       
608                 
609                 int numSeqs = 0;
610                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
611                         if(processors == 1){
612                                 numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
613                         }else{
614                                 numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile); 
615                         }
616                                 
617                         if (m->control_pressed) {  return 0; }
618                 #else
619                         numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
620                         if (m->control_pressed) {  return 0; }
621                 #endif
622
623                 sort(startPosition.begin(), startPosition.end());
624                 sort(endPosition.begin(), endPosition.end());
625                 sort(seqLength.begin(), seqLength.end());
626                 sort(ambigBases.begin(), ambigBases.end());
627                 sort(longHomoPolymer.begin(), longHomoPolymer.end());
628                 
629                 int criteriaPercentile  = int(numSeqs * (criteria / (float) 100));
630                 
631                 for (int i = 0; i < optimize.size(); i++) {
632                         if (optimize[i] == "start") { startPos = startPosition[criteriaPercentile]; m->mothurOut("Optimizing start to " + toString(startPos) + "."); m->mothurOutEndLine(); }
633                         else if (optimize[i] == "end") { int endcriteriaPercentile = int(numSeqs * ((100 - criteria) / (float) 100));  endPos = endPosition[endcriteriaPercentile]; m->mothurOut("Optimizing end to " + toString(endPos) + "."); m->mothurOutEndLine();}
634                         else if (optimize[i] == "maxambig") { maxAmbig = ambigBases[criteriaPercentile]; m->mothurOut("Optimizing maxambig to " + toString(maxAmbig) + "."); m->mothurOutEndLine(); }
635                         else if (optimize[i] == "maxhomop") { maxHomoP = longHomoPolymer[criteriaPercentile]; m->mothurOut("Optimizing maxhomop to " + toString(maxHomoP) + "."); m->mothurOutEndLine(); }
636                         else if (optimize[i] == "minlength") { int mincriteriaPercentile = int(numSeqs * ((100 - criteria) / (float) 100)); minLength = seqLength[mincriteriaPercentile]; m->mothurOut("Optimizing minlength to " + toString(minLength) + "."); m->mothurOutEndLine(); }
637                         else if (optimize[i] == "maxlength") { maxLength = seqLength[criteriaPercentile]; m->mothurOut("Optimizing maxlength to " + toString(maxLength) + "."); m->mothurOutEndLine(); }
638                 }
639
640                 return 0;
641         }
642         catch(exception& e) {
643                 m->errorOut(e, "ScreenSeqsCommand", "getSummary");
644                 exit(1);
645         }
646 }
647 /**************************************************************************************/
648 int ScreenSeqsCommand::driverCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, linePair* filePos) {   
649         try {
650                 
651                 ifstream in;
652                 m->openInputFile(filename, in);
653                                 
654                 in.seekg(filePos->start);
655
656                 bool done = false;
657                 int count = 0;
658         
659                 while (!done) {
660                                 
661                         if (m->control_pressed) { in.close(); return 1; }
662                                         
663                         Sequence current(in); m->gobble(in);
664         
665                         if (current.getName() != "") {
666                                 startPosition.push_back(current.getStartPos());
667                                 endPosition.push_back(current.getEndPos());
668                                 seqLength.push_back(current.getNumBases());
669                                 ambigBases.push_back(current.getAmbigBases());
670                                 longHomoPolymer.push_back(current.getLongHomoPolymer());
671                                 count++;
672                         }
673                         
674                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
675                                 unsigned long int pos = in.tellg();
676                                 if ((pos == -1) || (pos >= filePos->end)) { break; }
677                         #else
678                                 if (in.eof()) { break; }
679                         #endif
680                         
681                 }
682                 
683                 in.close();
684                 
685                 return count;
686         }
687         catch(exception& e) {
688                 m->errorOut(e, "ScreenSeqsCommand", "driverCreateSummary");
689                 exit(1);
690         }
691 }
692 /**************************************************************************************************/
693 int ScreenSeqsCommand::createProcessesCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename) {
694         try {
695 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
696                 int process = 1;
697                 int num = 0;
698                 processIDS.clear();
699                 
700                 //loop through and create all the processes you want
701                 while (process != processors) {
702                         int pid = fork();
703                         
704                         if (pid > 0) {
705                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
706                                 process++;
707                         }else if (pid == 0){
708                                 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[process]);
709                                 
710                                 //pass numSeqs to parent
711                                 ofstream out;
712                                 string tempFile = fastafile + toString(getpid()) + ".num.temp";
713                                 m->openOutputFile(tempFile, out);
714                                 
715                                 out << num << endl;
716                                 for (int k = 0; k < startPosition.size(); k++)          {               out << startPosition[k] << '\t'; }  out << endl;
717                                 for (int k = 0; k < endPosition.size(); k++)            {               out << endPosition[k] << '\t'; }  out << endl;
718                                 for (int k = 0; k < seqLength.size(); k++)                      {               out << seqLength[k] << '\t'; }  out << endl;
719                                 for (int k = 0; k < ambigBases.size(); k++)                     {               out << ambigBases[k] << '\t'; }  out << endl;
720                                 for (int k = 0; k < longHomoPolymer.size(); k++)        {               out << longHomoPolymer[k] << '\t'; }  out << endl;
721                                 
722                                 out.close();
723                                 
724                                 exit(0);
725                         }else { 
726                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
727                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
728                                 exit(0);
729                         }
730                 }
731                 
732                 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
733                 
734                 //force parent to wait until all the processes are done
735                 for (int i=0;i<processIDS.size();i++) { 
736                         int temp = processIDS[i];
737                         wait(&temp);
738                 }
739                 
740                 //parent reads in and combine Filter info
741                 for (int i = 0; i < processIDS.size(); i++) {
742                         string tempFilename = fastafile + toString(processIDS[i]) + ".num.temp";
743                         ifstream in;
744                         m->openInputFile(tempFilename, in);
745                         
746                         int temp, tempNum;
747                         in >> tempNum; m->gobble(in); num += tempNum;
748                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; startPosition.push_back(temp);              }               m->gobble(in);
749                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; endPosition.push_back(temp);                }               m->gobble(in);
750                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; seqLength.push_back(temp);                  }               m->gobble(in);
751                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; ambigBases.push_back(temp);                 }               m->gobble(in);
752                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; longHomoPolymer.push_back(temp);    }               m->gobble(in);
753                                 
754                         in.close();
755                         remove(tempFilename.c_str());
756                 }
757                 
758                 return num;
759 #endif          
760         }
761         catch(exception& e) {
762                 m->errorOut(e, "ScreenSeqsCommand", "createProcessesCreateSummary");
763                 exit(1);
764         }
765 }
766
767 //***************************************************************************************************************
768
769 int ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
770         try {
771                 ifstream inputGroups;
772                 m->openInputFile(groupfile, inputGroups);
773                 string seqName, group;
774                 set<string>::iterator it;
775                 
776                 string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
777                 outputNames.push_back(goodGroupFile);  outputTypes["group"].push_back(goodGroupFile);
778                 ofstream goodGroupOut;  m->openOutputFile(goodGroupFile, goodGroupOut);
779                 
780                 while(!inputGroups.eof()){
781                         if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); return 0; }
782
783                         inputGroups >> seqName >> group;
784                         it = badSeqNames.find(seqName);
785                         
786                         if(it != badSeqNames.end()){
787                                 badSeqNames.erase(it);
788                         }
789                         else{
790                                 goodGroupOut << seqName << '\t' << group << endl;
791                         }
792                         m->gobble(inputGroups);
793                 }
794                 
795                 if (m->control_pressed) { goodGroupOut.close();  inputGroups.close(); remove(goodGroupFile.c_str());  return 0; }
796
797                 //we were unable to remove some of the bad sequences
798                 if (badSeqNames.size() != 0) {
799                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
800                                 m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct."); 
801                                 m->mothurOutEndLine();
802                         }
803                 }
804                 
805                 inputGroups.close();
806                 goodGroupOut.close();
807                 
808                 if (m->control_pressed) { remove(goodGroupFile.c_str());   }
809                 
810                 return 0;
811         
812         }
813         catch(exception& e) {
814                 m->errorOut(e, "ScreenSeqsCommand", "screenGroupFile");
815                 exit(1);
816         }
817 }
818
819 //***************************************************************************************************************
820
821 int ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
822         try {
823                 ifstream inputAlignReport;
824                 m->openInputFile(alignreport, inputAlignReport);
825                 string seqName, group;
826                 set<string>::iterator it;
827                 
828                 string goodAlignReportFile = outputDir + m->getRootName(m->getSimpleName(alignreport)) + "good" + m->getExtension(alignreport);
829                 outputNames.push_back(goodAlignReportFile);  outputTypes["alignreport"].push_back(goodAlignReportFile);
830                 ofstream goodAlignReportOut;    m->openOutputFile(goodAlignReportFile, goodAlignReportOut);
831
832                 while (!inputAlignReport.eof()) {               //      need to copy header
833                         char c = inputAlignReport.get();
834                         goodAlignReportOut << c;
835                         if (c == 10 || c == 13){        break;  }       
836                 }
837
838                 while(!inputAlignReport.eof()){
839                         if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); return 0; }
840
841                         inputAlignReport >> seqName;
842                         it = badSeqNames.find(seqName);
843                         string line;            
844                         while (!inputAlignReport.eof()) {               //      need to copy header
845                                 char c = inputAlignReport.get();
846                                 line += c;
847                                 if (c == 10 || c == 13){        break;  }       
848                         }
849                         
850                         if(it != badSeqNames.end()){
851                                 badSeqNames.erase(it);
852                         }
853                         else{
854                                 goodAlignReportOut << seqName << '\t' << line;
855                         }
856                         m->gobble(inputAlignReport);
857                 }
858                 
859                 if (m->control_pressed) { goodAlignReportOut.close();  inputAlignReport.close(); remove(goodAlignReportFile.c_str());  return 0; }
860
861                 //we were unable to remove some of the bad sequences
862                 if (badSeqNames.size() != 0) {
863                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
864                                 m->mothurOut("Your alignreport file does not include the sequence " + *it + " please correct."); 
865                                 m->mothurOutEndLine();
866                         }
867                 }
868
869                 inputAlignReport.close();
870                 goodAlignReportOut.close();
871                                 
872                 if (m->control_pressed) {  remove(goodAlignReportFile.c_str());  return 0; }
873                 
874                 return 0;
875         
876         }
877         catch(exception& e) {
878                 m->errorOut(e, "ScreenSeqsCommand", "screenAlignReport");
879                 exit(1);
880         }
881         
882 }
883 //**********************************************************************************************************************
884
885 int ScreenSeqsCommand::driver(linePair* filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
886         try {
887                 ofstream goodFile;
888                 m->openOutputFile(goodFName, goodFile);
889                 
890                 ofstream badAccnosFile;
891                 m->openOutputFile(badAccnosFName, badAccnosFile);
892                 
893                 ifstream inFASTA;
894                 m->openInputFile(filename, inFASTA);
895
896                 inFASTA.seekg(filePos->start);
897
898                 bool done = false;
899                 int count = 0;
900         
901                 while (!done) {
902                 
903                         if (m->control_pressed) {  return 0; }
904                         
905                         Sequence currSeq(inFASTA); m->gobble(inFASTA);
906                         if (currSeq.getName() != "") {
907                                 bool goodSeq = 1;               //      innocent until proven guilty
908                                 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos())                  {       goodSeq = 0;    }
909                                 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos())                                {       goodSeq = 0;    }
910                                 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases())                {       goodSeq = 0;    }
911                                 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer())   {       goodSeq = 0;    }
912                                 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases())                {       goodSeq = 0;    }
913                                 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases())                {       goodSeq = 0;    }
914                                 
915                                 if(goodSeq == 1){
916                                         currSeq.printSequence(goodFile);        
917                                 }
918                                 else{
919                                         badAccnosFile << currSeq.getName() << endl;
920                                         badSeqNames.insert(currSeq.getName());
921                                 }
922                         count++;
923                         }
924                         
925                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
926                                 unsigned long int pos = inFASTA.tellg();
927                                 if ((pos == -1) || (pos >= filePos->end)) { break; }
928                         #else
929                                 if (inFASTA.eof()) { break; }
930                         #endif
931                         
932                         //report progress
933                         if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine();         }
934                 }
935                 //report progress
936                 if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine();         }
937                 
938                         
939                 goodFile.close();
940                 inFASTA.close();
941                 badAccnosFile.close();
942                 
943                 return count;
944         }
945         catch(exception& e) {
946                 m->errorOut(e, "ScreenSeqsCommand", "driver");
947                 exit(1);
948         }
949 }
950 //**********************************************************************************************************************
951 #ifdef USE_MPI
952 int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badAccnosFile, vector<unsigned long int>& MPIPos, set<string>& badSeqNames){
953         try {
954                 string outputString = "";
955                 MPI_Status statusGood; 
956                 MPI_Status statusBadAccnos; 
957                 MPI_Status status; 
958                 int pid;
959                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
960
961                 for(int i=0;i<num;i++){
962                 
963                         if (m->control_pressed) {  return 0; }
964                         
965                         //read next sequence
966                         int length = MPIPos[start+i+1] - MPIPos[start+i];
967
968                         char* buf4 = new char[length];
969                         memcpy(buf4, outputString.c_str(), length);
970
971                         MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
972                         
973                         string tempBuf = buf4;  delete buf4;
974                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
975                         istringstream iss (tempBuf,istringstream::in);
976                         
977                         Sequence currSeq(iss);                  
978                         
979                         //process seq
980                         if (currSeq.getName() != "") {
981                                 bool goodSeq = 1;               //      innocent until proven guilty
982                                 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos())                  {       goodSeq = 0;    }
983                                 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos())                                {       goodSeq = 0;    }
984                                 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases())                {       goodSeq = 0;    }
985                                 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer())   {       goodSeq = 0;    }
986                                 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases())                {       goodSeq = 0;    }
987                                 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases())                {       goodSeq = 0;    }
988                                 
989                                 if(goodSeq == 1){
990                                         outputString =  ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
991                                 
992                                         //print good seq
993                                         length = outputString.length();
994                                         char* buf2 = new char[length];
995                                         memcpy(buf2, outputString.c_str(), length);
996                                         
997                                         MPI_File_write_shared(goodFile, buf2, length, MPI_CHAR, &statusGood);
998                                         delete buf2;
999                                 }
1000                                 else{
1001
1002                                         badSeqNames.insert(currSeq.getName());
1003                                         
1004                                         //write to bad accnos file
1005                                         outputString = currSeq.getName() + "\n";
1006                                 
1007                                         length = outputString.length();
1008                                         char* buf3 = new char[length];
1009                                         memcpy(buf3, outputString.c_str(), length);
1010                                         
1011                                         MPI_File_write_shared(badAccnosFile, buf3, length, MPI_CHAR, &statusBadAccnos);
1012                                         delete buf3;
1013                                 }
1014                         }
1015                         
1016                         //report progress
1017                         if((i) % 100 == 0){     m->mothurOut("Processing sequence: " + toString(i)); m->mothurOutEndLine();             }
1018                 }
1019                                 
1020                 return 1;
1021         }
1022         catch(exception& e) {
1023                 m->errorOut(e, "ScreenSeqsCommand", "driverMPI");
1024                 exit(1);
1025         }
1026 }
1027 #endif
1028 /**************************************************************************************************/
1029
1030 int ScreenSeqsCommand::createProcesses(string goodFileName, string badAccnos, string filename, set<string>& badSeqNames) {
1031         try {
1032 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1033                 int process = 0;
1034                 int num = 0;
1035                 
1036                 //loop through and create all the processes you want
1037                 while (process != processors) {
1038                         int pid = fork();
1039                         
1040                         if (pid > 0) {
1041                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1042                                 process++;
1043                         }else if (pid == 0){
1044                                 num = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
1045                                 
1046                                 //pass numSeqs to parent
1047                                 ofstream out;
1048                                 string tempFile = filename + toString(getpid()) + ".num.temp";
1049                                 m->openOutputFile(tempFile, out);
1050                                 out << num << endl;
1051                                 out.close();
1052                                 
1053                                 exit(0);
1054                         }else { 
1055                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
1056                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1057                                 exit(0);
1058                         }
1059                 }
1060                 
1061                 //force parent to wait until all the processes are done
1062                 for (int i=0;i<processors;i++) { 
1063                         int temp = processIDS[i];
1064                         wait(&temp);
1065                 }
1066                 
1067                 for (int i = 0; i < processIDS.size(); i++) {
1068                         ifstream in;
1069                         string tempFile =  filename + toString(processIDS[i]) + ".num.temp";
1070                         m->openInputFile(tempFile, in);
1071                         if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
1072                         in.close(); remove(tempFile.c_str());
1073                 }
1074                 
1075                 return num;
1076 #endif          
1077         }
1078         catch(exception& e) {
1079                 m->errorOut(e, "ScreenSeqsCommand", "createProcesses");
1080                 exit(1);
1081         }
1082 }
1083
1084 //***************************************************************************************************************
1085
1086