]> git.donarmstrong.com Git - mothur.git/blob - screenseqscommand.cpp
added unix to ifdefs. minor changes while testing 1.24.0.
[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
12
13 //**********************************************************************************************************************
14 vector<string> ScreenSeqsCommand::setParameters(){      
15         try {
16                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
17                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
18                 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
19                 CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pqfile);
20                 CommandParameter palignreport("alignreport", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(palignreport);
21                 CommandParameter ptax("taxonomy", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(ptax);
22                 CommandParameter pstart("start", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pstart);
23                 CommandParameter pend("end", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pend);
24                 CommandParameter pmaxambig("maxambig", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pmaxambig);
25                 CommandParameter pmaxhomop("maxhomop", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pmaxhomop);
26                 CommandParameter pminlength("minlength", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pminlength);
27                 CommandParameter pmaxlength("maxlength", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pmaxlength);
28                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
29                 CommandParameter pcriteria("criteria", "Number", "", "90", "", "", "",false,false); parameters.push_back(pcriteria);
30                 CommandParameter poptimize("optimize", "Multiple", "none-start-end-maxambig-maxhomop-minlength-maxlength", "none", "", "", "",true,false); parameters.push_back(poptimize);
31                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
32                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
33                 
34                 vector<string> myArray;
35                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
36                 return myArray;
37         }
38         catch(exception& e) {
39                 m->errorOut(e, "ScreenSeqsCommand", "setParameters");
40                 exit(1);
41         }
42 }
43 //**********************************************************************************************************************
44 string ScreenSeqsCommand::getHelpString(){      
45         try {
46                 string helpString = "";
47                 helpString += "The screen.seqs command reads a fastafile and creates .....\n";
48                 helpString += "The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group, qfile, alignreport, taxonomy, optimize, criteria and processors.\n";
49                 helpString += "The fasta parameter is required.\n";
50                 helpString += "The alignreport and taxonomy parameters allow you to remove bad seqs from taxonomy and alignreport files.\n";
51                 helpString += "The start parameter .... The default is -1.\n";
52                 helpString += "The end parameter .... The default is -1.\n";
53                 helpString += "The maxambig parameter allows you to set the maximum number of ambigious bases allowed. The default is -1.\n";
54                 helpString += "The maxhomop parameter allows you to set a maximum homopolymer length. \n";
55                 helpString += "The minlength parameter allows you to set and minimum sequence length. \n";
56                 helpString += "The maxlength parameter allows you to set and maximum sequence length. \n";
57                 helpString += "The processors parameter allows you to specify the number of processors to use while running the command. The default is 1.\n";
58                 helpString += "The optimize and criteria parameters allow you set the start, end, maxabig, maxhomop, minlength and maxlength parameters relative to your set of sequences .\n";
59                 helpString += "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";
60                 helpString += "The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n";
61                 helpString += "The screen.seqs command should be in the following format: \n";
62                 helpString += "screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig,  \n";
63                 helpString += "maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength)  \n";    
64                 helpString += "Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n";
65                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
66                 return helpString;
67         }
68         catch(exception& e) {
69                 m->errorOut(e, "ScreenSeqsCommand", "getHelpString");
70                 exit(1);
71         }
72 }
73 //**********************************************************************************************************************
74 ScreenSeqsCommand::ScreenSeqsCommand(){ 
75         try {
76                 abort = true; calledHelp = true; 
77                 setParameters();
78                 vector<string> tempOutNames;
79                 outputTypes["fasta"] = tempOutNames;
80                 outputTypes["name"] = tempOutNames;
81                 outputTypes["group"] = tempOutNames;
82                 outputTypes["alignreport"] = tempOutNames;
83                 outputTypes["accnos"] = tempOutNames;
84                 outputTypes["qfile"] = tempOutNames;
85                 outputTypes["taxonomy"] = tempOutNames;
86         }
87         catch(exception& e) {
88                 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
89                 exit(1);
90         }
91 }
92 //***************************************************************************************************************
93
94 ScreenSeqsCommand::ScreenSeqsCommand(string option)  {
95         try {
96                 abort = false; calledHelp = false;   
97                 
98                 //allow user to run help
99                 if(option == "help") { help(); abort = true; calledHelp = true; }
100                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
101                 
102                 else {
103                         vector<string> myArray = setParameters();
104                         
105                         OptionParser parser(option);
106                         map<string,string> parameters = parser.getParameters();
107                         
108                         ValidParameters validParameter("screen.seqs");
109                         map<string,string>::iterator it;
110                         
111                         //check to make sure all parameters are valid for command
112                         for (it = parameters.begin(); it != parameters.end(); it++) { 
113                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
114                         }
115                         
116                         //initialize outputTypes
117                         vector<string> tempOutNames;
118                         outputTypes["fasta"] = tempOutNames;
119                         outputTypes["name"] = tempOutNames;
120                         outputTypes["group"] = tempOutNames;
121                         outputTypes["alignreport"] = tempOutNames;
122                         outputTypes["accnos"] = tempOutNames;
123                         outputTypes["qfile"] = tempOutNames;
124                         outputTypes["taxonomy"] = tempOutNames;
125                         
126                         //if the user changes the input directory command factory will send this info to us in the output parameter 
127                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
128                         if (inputDir == "not found"){   inputDir = "";          }
129                         else {
130                                 string path;
131                                 it = parameters.find("fasta");
132                                 //user has given a template file
133                                 if(it != parameters.end()){ 
134                                         path = m->hasPath(it->second);
135                                         //if the user has not given a path then, add inputdir. else leave path alone.
136                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
137                                 }
138                                 
139                                 it = parameters.find("group");
140                                 //user has given a template file
141                                 if(it != parameters.end()){ 
142                                         path = m->hasPath(it->second);
143                                         //if the user has not given a path then, add inputdir. else leave path alone.
144                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
145                                 }
146                                 
147                                 it = parameters.find("name");
148                                 //user has given a template file
149                                 if(it != parameters.end()){ 
150                                         path = m->hasPath(it->second);
151                                         //if the user has not given a path then, add inputdir. else leave path alone.
152                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
153                                 }
154                                 
155                                 it = parameters.find("alignreport");
156                                 //user has given a template file
157                                 if(it != parameters.end()){ 
158                                         path = m->hasPath(it->second);
159                                         //if the user has not given a path then, add inputdir. else leave path alone.
160                                         if (path == "") {       parameters["alignreport"] = inputDir + it->second;              }
161                                 }
162                                 
163                                 it = parameters.find("qfile");
164                                 //user has given a template file
165                                 if(it != parameters.end()){ 
166                                         path = m->hasPath(it->second);
167                                         //if the user has not given a path then, add inputdir. else leave path alone.
168                                         if (path == "") {       parameters["qfile"] = inputDir + it->second;            }
169                                 }
170                                 
171                                 it = parameters.find("taxonomy");
172                                 //user has given a template file
173                                 if(it != parameters.end()){ 
174                                         path = m->hasPath(it->second);
175                                         //if the user has not given a path then, add inputdir. else leave path alone.
176                                         if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
177                                 }
178                         }
179
180                         //check for required parameters
181                         fastafile = validParameter.validFile(parameters, "fasta", true);
182                         if (fastafile == "not found") {                         
183                                 fastafile = m->getFastaFile(); 
184                                 if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
185                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
186                         }
187                         else if (fastafile == "not open") { abort = true; }
188                         else { m->setFastaFile(fastafile); }
189         
190                         groupfile = validParameter.validFile(parameters, "group", true);
191                         if (groupfile == "not open") { abort = true; }  
192                         else if (groupfile == "not found") { groupfile = ""; }
193                         else { m->setGroupFile(groupfile); }
194                         
195                         qualfile = validParameter.validFile(parameters, "qfile", true);
196                         if (qualfile == "not open") { abort = true; }   
197                         else if (qualfile == "not found") { qualfile = ""; }
198                         else { m->setQualFile(qualfile); }
199                         
200                         namefile = validParameter.validFile(parameters, "name", true);
201                         if (namefile == "not open") { namefile = ""; abort = true; }
202                         else if (namefile == "not found") { namefile = ""; }    
203                         else { m->setNameFile(namefile); }
204                         
205                         alignreport = validParameter.validFile(parameters, "alignreport", true);
206                         if (alignreport == "not open") { abort = true; }
207                         else if (alignreport == "not found") { alignreport = ""; }
208                         
209                         taxonomy = validParameter.validFile(parameters, "taxonomy", true);
210                         if (taxonomy == "not open") { abort = true; }
211                         else if (taxonomy == "not found") { taxonomy = ""; }    
212                         
213                         //if the user changes the output directory command factory will send this info to us in the output parameter 
214                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  
215                                 outputDir = ""; 
216                                 outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it       
217                         }
218
219                         //check for optional parameter and set defaults
220                         // ...at some point should added some additional type checking...
221                         string temp;
222                         temp = validParameter.validFile(parameters, "start", false);            if (temp == "not found") { temp = "-1"; }
223                         m->mothurConvert(temp, startPos); 
224                 
225                         temp = validParameter.validFile(parameters, "end", false);                      if (temp == "not found") { temp = "-1"; }
226                         m->mothurConvert(temp, endPos);  
227
228                         temp = validParameter.validFile(parameters, "maxambig", false);         if (temp == "not found") { temp = "-1"; }
229                         m->mothurConvert(temp, maxAmbig);  
230
231                         temp = validParameter.validFile(parameters, "maxhomop", false);         if (temp == "not found") { temp = "-1"; }
232                         m->mothurConvert(temp, maxHomoP);  
233
234                         temp = validParameter.validFile(parameters, "minlength", false);        if (temp == "not found") { temp = "-1"; }
235                         m->mothurConvert(temp, minLength); 
236                         
237                         temp = validParameter.validFile(parameters, "maxlength", false);        if (temp == "not found") { temp = "-1"; }
238                         m->mothurConvert(temp, maxLength); 
239                         
240                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
241                         m->setProcessors(temp);
242                         m->mothurConvert(temp, processors);
243                         
244                         temp = validParameter.validFile(parameters, "optimize", false); //optimizing trumps the optimized values original value
245                         if (temp == "not found"){       temp = "none";          }
246                         m->splitAtDash(temp, optimize);         
247                         
248                         //check for invalid optimize options
249                         set<string> validOptimizers;
250                         validOptimizers.insert("none"); validOptimizers.insert("start"); validOptimizers.insert("end"); validOptimizers.insert("maxambig"); validOptimizers.insert("maxhomop"); validOptimizers.insert("minlength"); validOptimizers.insert("maxlength");
251                         for (int i = 0; i < optimize.size(); i++) { 
252                                 if (validOptimizers.count(optimize[i]) == 0) { 
253                                         m->mothurOut(optimize[i] + " is not a valid optimizer. Valid options are start, end, maxambig, maxhomop, minlength and maxlength."); m->mothurOutEndLine();
254                                         optimize.erase(optimize.begin()+i);
255                                         i--;
256                                 }
257                         }
258                         
259                         if (optimize.size() == 1) { if (optimize[0] == "none") { optimize.clear(); } }
260                         
261                         temp = validParameter.validFile(parameters, "criteria", false); if (temp == "not found"){       temp = "90";                            }
262                         m->mothurConvert(temp, criteria); 
263                         
264                         if (namefile == "") {
265                                 vector<string> files; files.push_back(fastafile);
266                                 parser.getNameFile(files);
267                         }
268                 }
269
270         }
271         catch(exception& e) {
272                 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
273                 exit(1);
274         }
275 }
276 //***************************************************************************************************************
277
278 int ScreenSeqsCommand::execute(){
279         try{
280                 
281                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
282                 
283                 //if the user want to optimize we need to know the 90% mark
284                 vector<unsigned long long> positions;
285                 if (optimize.size() != 0) {  //get summary is paralellized so we need to divideFile, no need to do this step twice so I moved it here
286                         //use the namefile to optimize correctly
287                         if (namefile != "") { nameMap = m->readNames(namefile); }
288                         getSummary(positions); 
289                 } 
290                 else { 
291                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
292                 positions = m->divideFile(fastafile, processors);
293                 for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }
294                         #else 
295                 if(processors == 1){ lines.push_back(linePair(0, 1000));  }
296                 else {
297                     int numFastaSeqs = 0;
298                     positions = m->setFilePosFasta(fastafile, numFastaSeqs); 
299                     if (positions.size() < processors) { processors = positions.size(); }
300                 
301                     //figure out how many sequences you have to process
302                     int numSeqsPerProcessor = numFastaSeqs / processors;
303                     for (int i = 0; i < processors; i++) {
304                         int startIndex =  i * numSeqsPerProcessor;
305                         if(i == (processors - 1)){      numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;   }
306                         lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
307                     }
308                 }
309                         #endif
310                 }
311                                         
312                 string goodSeqFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "good" + m->getExtension(fastafile);
313                 string badAccnosFile =  outputDir + m->getRootName(m->getSimpleName(fastafile)) + "bad.accnos";
314                 
315                 int numFastaSeqs = 0;
316                 set<string> badSeqNames;
317                 int start = time(NULL);
318         
319 #ifdef USE_MPI  
320                         int pid, numSeqsPerProcessor; 
321                         int tag = 2001;
322                         vector<unsigned long long> MPIPos;
323                         
324                         MPI_Status status; 
325                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
326                         MPI_Comm_size(MPI_COMM_WORLD, &processors); 
327         
328                         MPI_File inMPI;
329                         MPI_File outMPIGood;
330                         MPI_File outMPIBadAccnos;
331                         
332                         int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
333                         int inMode=MPI_MODE_RDONLY; 
334                         
335                         char outGoodFilename[1024];
336                         strcpy(outGoodFilename, goodSeqFile.c_str());
337
338                         char outBadAccnosFilename[1024];
339                         strcpy(outBadAccnosFilename, badAccnosFile.c_str());
340
341                         char inFileName[1024];
342                         strcpy(inFileName, fastafile.c_str());
343                         
344                         MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
345                         MPI_File_open(MPI_COMM_WORLD, outGoodFilename, outMode, MPI_INFO_NULL, &outMPIGood);
346                         MPI_File_open(MPI_COMM_WORLD, outBadAccnosFilename, outMode, MPI_INFO_NULL, &outMPIBadAccnos);
347                         
348                         if (m->control_pressed) { MPI_File_close(&inMPI);  MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
349                         
350                         if (pid == 0) { //you are the root process 
351                                 
352                                 MPIPos = m->setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
353                                 
354                                 //send file positions to all processes
355                                 for(int i = 1; i < processors; i++) { 
356                                         MPI_Send(&numFastaSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
357                                         MPI_Send(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
358                                 }
359                                 
360                                 //figure out how many sequences you have to align
361                                 numSeqsPerProcessor = numFastaSeqs / processors;
362                                 int startIndex =  pid * numSeqsPerProcessor;
363                                 if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
364
365                                 //align your part
366                                 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
367
368                                 if (m->control_pressed) { MPI_File_close(&inMPI);  MPI_File_close(&outMPIGood);  MPI_File_close(&outMPIBadAccnos);  return 0; }
369
370                                 for (int i = 1; i < processors; i++) {
371                                         //get bad lists
372                                         int badSize;
373                                         MPI_Recv(&badSize, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
374                                 }
375                         }else{ //you are a child process
376                                 MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
377                                 MPIPos.resize(numFastaSeqs+1);
378                                 MPI_Recv(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
379
380                                 //figure out how many sequences you have to align
381                                 numSeqsPerProcessor = numFastaSeqs / processors;
382                                 int startIndex =  pid * numSeqsPerProcessor;
383                                 if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
384
385                                 //align your part
386                                 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
387
388                                 if (m->control_pressed) { MPI_File_close(&inMPI);  MPI_File_close(&outMPIGood);  MPI_File_close(&outMPIBadAccnos); return 0; }
389                                 
390                                 //send bad list 
391                                 int badSize = badSeqNames.size();
392                                 MPI_Send(&badSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
393                         }
394                         
395                         //close files 
396                         MPI_File_close(&inMPI);
397                         MPI_File_close(&outMPIGood);
398                         MPI_File_close(&outMPIBadAccnos);
399                         MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
400                                         
401 #else
402         if(processors == 1){ numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);       }       
403         else{ numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames); }
404         
405         if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
406 #endif          
407
408                 #ifdef USE_MPI
409                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
410                                         
411                         if (pid == 0) { //only one process should fix files
412                         
413                                 //read accnos file with all names in it, process 0 just has its names
414                                 MPI_File inMPIAccnos;
415                                 MPI_Offset size;
416                         
417                                 char inFileName[1024];
418                                 strcpy(inFileName, badAccnosFile.c_str());
419                         
420                                 MPI_File_open(MPI_COMM_SELF, inFileName, inMode, MPI_INFO_NULL, &inMPIAccnos);  //comm, filename, mode, info, filepointer
421                                 MPI_File_get_size(inMPIAccnos, &size);
422                         
423                                 char* buffer = new char[size];
424                                 MPI_File_read(inMPIAccnos, buffer, size, MPI_CHAR, &status);
425                         
426                                 string tempBuf = buffer;
427                                 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size);  }
428                                 istringstream iss (tempBuf,istringstream::in);
429
430                                 delete buffer;
431                                 MPI_File_close(&inMPIAccnos);
432                                 
433                                 badSeqNames.clear();
434                                 string tempName;
435                                 while (!iss.eof()) {
436                                         iss >> tempName; m->gobble(iss);
437                                         badSeqNames.insert(tempName);
438                                 }
439                 #endif
440                                                                                                                                                                         
441                 if(namefile != "" && groupfile != "")   {       
442                         screenNameGroupFile(badSeqNames);       
443                         if (m->control_pressed) {  m->mothurRemove(goodSeqFile); return 0; }
444                 }else if(namefile != "")        {       
445                         screenNameGroupFile(badSeqNames);
446                         if (m->control_pressed) {  m->mothurRemove(goodSeqFile);  return 0; }   
447                 }else if(groupfile != "")                               {       screenGroupFile(badSeqNames);           }       // this screens just the group
448                 
449                 if (m->control_pressed) { m->mothurRemove(goodSeqFile);  return 0; }
450
451                 if(alignreport != "")                                   {       screenAlignReport(badSeqNames);         }
452                 if(qualfile != "")                                              {       screenQual(badSeqNames);                        }
453                 if(taxonomy != "")                                              {       screenTaxonomy(badSeqNames);            }
454                 
455                 if (m->control_pressed) { m->mothurRemove(goodSeqFile);  return 0; }
456                 
457                 #ifdef USE_MPI
458                         }
459                 #endif
460
461                 m->mothurOutEndLine();
462                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
463                 m->mothurOut(goodSeqFile); m->mothurOutEndLine();       outputTypes["fasta"].push_back(goodSeqFile);
464                 m->mothurOut(badAccnosFile); m->mothurOutEndLine();      outputTypes["accnos"].push_back(badAccnosFile);
465                 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
466                 m->mothurOutEndLine();
467                 m->mothurOutEndLine();
468                 
469                 //set fasta file as new current fastafile
470                 string current = "";
471                 itTypes = outputTypes.find("fasta");
472                 if (itTypes != outputTypes.end()) {
473                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
474                 }
475                 
476                 itTypes = outputTypes.find("name");
477                 if (itTypes != outputTypes.end()) {
478                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
479                 }
480                 
481                 itTypes = outputTypes.find("group");
482                 if (itTypes != outputTypes.end()) {
483                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
484                 }
485                 
486                 itTypes = outputTypes.find("qfile");
487                 if (itTypes != outputTypes.end()) {
488                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setQualFile(current); }
489                 }
490                 
491                 itTypes = outputTypes.find("taxonomy");
492                 if (itTypes != outputTypes.end()) {
493                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
494                 }
495
496                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
497                 m->mothurOutEndLine();
498
499                 return 0;
500         }
501         catch(exception& e) {
502                 m->errorOut(e, "ScreenSeqsCommand", "execute");
503                 exit(1);
504         }
505 }
506
507 //***************************************************************************************************************
508
509 int ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
510         try {
511                 ifstream inputNames;
512                 m->openInputFile(namefile, inputNames);
513                 set<string> badSeqGroups;
514                 string seqName, seqList, group;
515                 set<string>::iterator it;
516
517                 string goodNameFile = outputDir + m->getRootName(m->getSimpleName(namefile)) + "good" + m->getExtension(namefile);
518                 outputNames.push_back(goodNameFile);  outputTypes["name"].push_back(goodNameFile);
519                 
520                 ofstream goodNameOut;   m->openOutputFile(goodNameFile, goodNameOut);
521                 
522                 while(!inputNames.eof()){
523                         if (m->control_pressed) { goodNameOut.close();  inputNames.close(); m->mothurRemove(goodNameFile);  return 0; }
524
525                         inputNames >> seqName >> seqList;
526                         it = badSeqNames.find(seqName);
527                                 
528                         if(it != badSeqNames.end()){
529                                 badSeqNames.erase(it);
530                                 
531                                 if(namefile != ""){
532                                         int start = 0;
533                                         for(int i=0;i<seqList.length();i++){
534                                                 if(seqList[i] == ','){
535                                                         badSeqGroups.insert(seqList.substr(start,i-start));
536                                                         start = i+1;
537                                                 }                                       
538                                         }
539                                         badSeqGroups.insert(seqList.substr(start,seqList.length()-start));
540                                 }
541                         }
542                         else{
543                                 goodNameOut << seqName << '\t' << seqList << endl;
544                         }
545                         m->gobble(inputNames);
546                 }
547                 inputNames.close();
548                 goodNameOut.close();
549         
550                 //we were unable to remove some of the bad sequences
551                 if (badSeqNames.size() != 0) {
552                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
553                                 m->mothurOut("Your namefile does not include the sequence " + *it + " please correct."); 
554                                 m->mothurOutEndLine();
555                         }
556                 }
557
558                 if(groupfile != ""){
559                         
560                         ifstream inputGroups;
561                         m->openInputFile(groupfile, inputGroups);
562
563                         string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
564                         outputNames.push_back(goodGroupFile);   outputTypes["group"].push_back(goodGroupFile);
565                         
566                         ofstream goodGroupOut;  m->openOutputFile(goodGroupFile, goodGroupOut);
567                         
568                         while(!inputGroups.eof()){
569                                 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodNameFile);  m->mothurRemove(goodGroupFile); return 0; }
570
571                                 inputGroups >> seqName >> group;
572                                 
573                                 it = badSeqGroups.find(seqName);
574                                 
575                                 if(it != badSeqGroups.end()){
576                                         badSeqGroups.erase(it);
577                                 }
578                                 else{
579                                         goodGroupOut << seqName << '\t' << group << endl;
580                                 }
581                                 m->gobble(inputGroups);
582                         }
583                         inputGroups.close();
584                         goodGroupOut.close();
585                         
586                         //we were unable to remove some of the bad sequences
587                         if (badSeqGroups.size() != 0) {
588                                 for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {  
589                                         m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct."); 
590                                         m->mothurOutEndLine();
591                                 }
592                         }
593                 }
594                 
595                 
596                 return 0;
597         
598         }
599         catch(exception& e) {
600                 m->errorOut(e, "ScreenSeqsCommand", "screenNameGroupFile");
601                 exit(1);
602         }
603 }
604 //***************************************************************************************************************
605 int ScreenSeqsCommand::getSummary(vector<unsigned long long>& positions){
606         try {
607                 
608                 vector<int> startPosition;
609                 vector<int> endPosition;
610                 vector<int> seqLength;
611                 vector<int> ambigBases;
612                 vector<int> longHomoPolymer;
613                 
614         vector<unsigned long long> positions;
615 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
616                 positions = m->divideFile(fastafile, processors);
617                 for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }   
618 #else
619                 if(processors == 1){ lines.push_back(linePair(0, 1000));  }
620         else {
621             int numFastaSeqs = 0;
622             positions = m->setFilePosFasta(fastafile, numFastaSeqs); 
623             if (positions.size() < processors) { processors = positions.size(); }
624             
625             //figure out how many sequences you have to process
626             int numSeqsPerProcessor = numFastaSeqs / processors;
627             for (int i = 0; i < processors; i++) {
628                 int startIndex =  i * numSeqsPerProcessor;
629                 if(i == (processors - 1)){      numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;   }
630                 lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
631             }
632         }
633 #endif
634                 
635 #ifdef USE_MPI
636                 int pid;
637                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
638                 
639                 if (pid == 0) { 
640                         driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
641 #else
642                 int numSeqs = 0;
643                 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
644                         if(processors == 1){
645                                 numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
646                         }else{
647                                 numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile); 
648                         }
649                                 
650                         if (m->control_pressed) {  return 0; }
651                 //#else
652                 //      numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
653                 //      if (m->control_pressed) {  return 0; }
654                 //#endif
655 #endif
656                 sort(startPosition.begin(), startPosition.end());
657                 sort(endPosition.begin(), endPosition.end());
658                 sort(seqLength.begin(), seqLength.end());
659                 sort(ambigBases.begin(), ambigBases.end());
660                 sort(longHomoPolymer.begin(), longHomoPolymer.end());
661                 
662                 //numSeqs is the number of unique seqs, startPosition.size() is the total number of seqs, we want to optimize using all seqs
663                 int criteriaPercentile  = int(startPosition.size() * (criteria / (float) 100));
664                 
665                 for (int i = 0; i < optimize.size(); i++) {
666                         if (optimize[i] == "start") { startPos = startPosition[criteriaPercentile]; m->mothurOut("Optimizing start to " + toString(startPos) + "."); m->mothurOutEndLine(); }
667                         else if (optimize[i] == "end") { int endcriteriaPercentile = int(endPosition.size() * ((100 - criteria) / (float) 100));  endPos = endPosition[endcriteriaPercentile]; m->mothurOut("Optimizing end to " + toString(endPos) + "."); m->mothurOutEndLine();}
668                         else if (optimize[i] == "maxambig") { maxAmbig = ambigBases[criteriaPercentile]; m->mothurOut("Optimizing maxambig to " + toString(maxAmbig) + "."); m->mothurOutEndLine(); }
669                         else if (optimize[i] == "maxhomop") { maxHomoP = longHomoPolymer[criteriaPercentile]; m->mothurOut("Optimizing maxhomop to " + toString(maxHomoP) + "."); m->mothurOutEndLine(); }
670                         else if (optimize[i] == "minlength") { int mincriteriaPercentile = int(seqLength.size() * ((100 - criteria) / (float) 100)); minLength = seqLength[mincriteriaPercentile]; m->mothurOut("Optimizing minlength to " + toString(minLength) + "."); m->mothurOutEndLine(); }
671                         else if (optimize[i] == "maxlength") { maxLength = seqLength[criteriaPercentile]; m->mothurOut("Optimizing maxlength to " + toString(maxLength) + "."); m->mothurOutEndLine(); }
672                 }
673
674 #ifdef USE_MPI
675                 }
676                 
677                 MPI_Status status; 
678                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
679                 MPI_Comm_size(MPI_COMM_WORLD, &processors); 
680                         
681                 if (pid == 0) { 
682                         //send file positions to all processes
683                         for(int i = 1; i < processors; i++) { 
684                                 MPI_Send(&startPos, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
685                                 MPI_Send(&endPos, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
686                                 MPI_Send(&maxAmbig, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
687                                 MPI_Send(&maxHomoP, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
688                                 MPI_Send(&minLength, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
689                                 MPI_Send(&maxLength, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
690                         }
691                 }else {
692                         MPI_Recv(&startPos, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
693                         MPI_Recv(&endPos, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
694                         MPI_Recv(&maxAmbig, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
695                         MPI_Recv(&maxHomoP, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
696                         MPI_Recv(&minLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
697                         MPI_Recv(&maxLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
698                 }
699                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
700 #endif
701                 return 0;
702         }
703         catch(exception& e) {
704                 m->errorOut(e, "ScreenSeqsCommand", "getSummary");
705                 exit(1);
706         }
707 }
708 /**************************************************************************************/
709 int ScreenSeqsCommand::driverCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, linePair filePos) {    
710         try {
711                 
712                 ifstream in;
713                 m->openInputFile(filename, in);
714                                 
715                 in.seekg(filePos.start);
716
717                 bool done = false;
718                 int count = 0;
719         
720                 while (!done) {
721                                 
722                         if (m->control_pressed) { in.close(); return 1; }
723                                         
724                         Sequence current(in); m->gobble(in);
725         
726                         if (current.getName() != "") {
727                                 int num = 1;
728                                 if (namefile != "") {
729                                         //make sure this sequence is in the namefile, else error 
730                                         map<string, int>::iterator it = nameMap.find(current.getName());
731                                         
732                                         if (it == nameMap.end()) { m->mothurOut("[ERROR]: " + current.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); m->control_pressed = true; }
733                                         else { num = it->second; }
734                                 }
735                                 
736                                 //for each sequence this sequence represents
737                                 for (int i = 0; i < num; i++) {
738                                         startPosition.push_back(current.getStartPos());
739                                         endPosition.push_back(current.getEndPos());
740                                         seqLength.push_back(current.getNumBases());
741                                         ambigBases.push_back(current.getAmbigBases());
742                                         longHomoPolymer.push_back(current.getLongHomoPolymer());
743                                 }
744                                 
745                                 count++;
746                         }
747                         //if((count) % 100 == 0){       m->mothurOut("Optimizing sequence: " + toString(count)); m->mothurOutEndLine();         }
748                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
749                                 unsigned long long pos = in.tellg();
750                                 if ((pos == -1) || (pos >= filePos.end)) { break; }
751                         #else
752                                 if (in.eof()) { break; }
753                         #endif
754                         
755                 }
756                 
757                 in.close();
758                 
759                 return count;
760         }
761         catch(exception& e) {
762                 m->errorOut(e, "ScreenSeqsCommand", "driverCreateSummary");
763                 exit(1);
764         }
765 }
766 /**************************************************************************************************/
767 int ScreenSeqsCommand::createProcessesCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename) {
768         try {
769         
770         int process = 1;
771                 int num = 0;
772                 vector<int> processIDS;
773
774 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
775                                 
776                 //loop through and create all the processes you want
777                 while (process != processors) {
778                         int pid = fork();
779                         
780                         if (pid > 0) {
781                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
782                                 process++;
783                         }else if (pid == 0){
784                                 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[process]);
785                                 
786                                 //pass numSeqs to parent
787                                 ofstream out;
788                                 string tempFile = fastafile + toString(getpid()) + ".num.temp";
789                                 m->openOutputFile(tempFile, out);
790                                 
791                                 out << num << endl;
792                                 out << startPosition.size() << endl;
793                                 for (int k = 0; k < startPosition.size(); k++)          {               out << startPosition[k] << '\t'; }  out << endl;
794                                 for (int k = 0; k < endPosition.size(); k++)            {               out << endPosition[k] << '\t'; }  out << endl;
795                                 for (int k = 0; k < seqLength.size(); k++)                      {               out << seqLength[k] << '\t'; }  out << endl;
796                                 for (int k = 0; k < ambigBases.size(); k++)                     {               out << ambigBases[k] << '\t'; }  out << endl;
797                                 for (int k = 0; k < longHomoPolymer.size(); k++)        {               out << longHomoPolymer[k] << '\t'; }  out << endl;
798                                 
799                                 out.close();
800                                 
801                                 exit(0);
802                         }else { 
803                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
804                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
805                                 exit(0);
806                         }
807                 }
808                 
809                 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
810                 
811                 //force parent to wait until all the processes are done
812                 for (int i=0;i<processIDS.size();i++) { 
813                         int temp = processIDS[i];
814                         wait(&temp);
815                 }
816                 
817                 //parent reads in and combine Filter info
818                 for (int i = 0; i < processIDS.size(); i++) {
819                         string tempFilename = fastafile + toString(processIDS[i]) + ".num.temp";
820                         ifstream in;
821                         m->openInputFile(tempFilename, in);
822                         
823                         int temp, tempNum;
824                         in >> tempNum; m->gobble(in); num += tempNum;
825                         in >> tempNum; m->gobble(in);
826                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; startPosition.push_back(temp);              }               m->gobble(in);
827                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; endPosition.push_back(temp);                }               m->gobble(in);
828                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; seqLength.push_back(temp);                  }               m->gobble(in);
829                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; ambigBases.push_back(temp);                 }               m->gobble(in);
830                         for (int k = 0; k < tempNum; k++)                       {               in >> temp; longHomoPolymer.push_back(temp);    }               m->gobble(in);
831                                 
832                         in.close();
833                         m->mothurRemove(tempFilename);
834                 }
835                 
836                 
837 #else 
838         //////////////////////////////////////////////////////////////////////////////////////////////////////
839                 //Windows version shared memory, so be careful when passing variables through the seqSumData struct. 
840                 //Above fork() will clone, so memory is separate, but that's not the case with windows, 
841                 //Taking advantage of shared memory to allow both threads to add info to vectors.
842                 //////////////////////////////////////////////////////////////////////////////////////////////////////
843                 
844                 vector<sumData*> pDataArray; 
845                 DWORD   dwThreadIdArray[processors-1];
846                 HANDLE  hThreadArray[processors-1]; 
847                 
848                 //Create processor worker threads.
849                 for( int i=0; i<processors-1; i++ ){
850             
851                         // Allocate memory for thread data.
852                         sumData* tempSum = new sumData(filename, m, lines[i].start, lines[i].end, namefile, nameMap);
853                         pDataArray.push_back(tempSum);
854                         
855                         //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
856                         //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
857                         hThreadArray[i] = CreateThread(NULL, 0, MySumThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);   
858                 }
859                 
860         //do your part
861                 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[processors-1]);
862          
863                 //Wait until all threads have terminated.
864                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
865                 
866                 //Close all thread handles and free memory allocations.
867                 for(int i=0; i < pDataArray.size(); i++){
868                         num += pDataArray[i]->count;
869             for (int k = 0; k < pDataArray[i]->startPosition.size(); k++) {     startPosition.push_back(pDataArray[i]->startPosition[k]);       }
870                         for (int k = 0; k < pDataArray[i]->endPosition.size(); k++) {   endPosition.push_back(pDataArray[i]->endPosition[k]);       }
871             for (int k = 0; k < pDataArray[i]->seqLength.size(); k++) { seqLength.push_back(pDataArray[i]->seqLength[k]);       }
872             for (int k = 0; k < pDataArray[i]->ambigBases.size(); k++) {        ambigBases.push_back(pDataArray[i]->ambigBases[k]);       }
873             for (int k = 0; k < pDataArray[i]->longHomoPolymer.size(); k++) {   longHomoPolymer.push_back(pDataArray[i]->longHomoPolymer[k]);       }
874                         CloseHandle(hThreadArray[i]);
875                         delete pDataArray[i];
876                 }
877
878 #endif          
879         return num;
880         }
881         catch(exception& e) {
882                 m->errorOut(e, "ScreenSeqsCommand", "createProcessesCreateSummary");
883                 exit(1);
884         }
885 }
886
887 //***************************************************************************************************************
888
889 int ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
890         try {
891                 ifstream inputGroups;
892                 m->openInputFile(groupfile, inputGroups);
893                 string seqName, group;
894                 set<string>::iterator it;
895                 
896                 string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
897                 outputNames.push_back(goodGroupFile);  outputTypes["group"].push_back(goodGroupFile);
898                 ofstream goodGroupOut;  m->openOutputFile(goodGroupFile, goodGroupOut);
899                 
900                 while(!inputGroups.eof()){
901                         if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodGroupFile); return 0; }
902
903                         inputGroups >> seqName >> group;
904                         it = badSeqNames.find(seqName);
905                         
906                         if(it != badSeqNames.end()){
907                                 badSeqNames.erase(it);
908                         }
909                         else{
910                                 goodGroupOut << seqName << '\t' << group << endl;
911                         }
912                         m->gobble(inputGroups);
913                 }
914                 
915                 if (m->control_pressed) { goodGroupOut.close();  inputGroups.close(); m->mothurRemove(goodGroupFile);  return 0; }
916
917                 //we were unable to remove some of the bad sequences
918                 if (badSeqNames.size() != 0) {
919                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
920                                 m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct."); 
921                                 m->mothurOutEndLine();
922                         }
923                 }
924                 
925                 inputGroups.close();
926                 goodGroupOut.close();
927                 
928                 if (m->control_pressed) { m->mothurRemove(goodGroupFile);   }
929                 
930                 return 0;
931         
932         }
933         catch(exception& e) {
934                 m->errorOut(e, "ScreenSeqsCommand", "screenGroupFile");
935                 exit(1);
936         }
937 }
938
939 //***************************************************************************************************************
940
941 int ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
942         try {
943                 ifstream inputAlignReport;
944                 m->openInputFile(alignreport, inputAlignReport);
945                 string seqName, group;
946                 set<string>::iterator it;
947                 
948                 string goodAlignReportFile = outputDir + m->getRootName(m->getSimpleName(alignreport)) + "good" + m->getExtension(alignreport);
949                 outputNames.push_back(goodAlignReportFile);  outputTypes["alignreport"].push_back(goodAlignReportFile);
950                 ofstream goodAlignReportOut;    m->openOutputFile(goodAlignReportFile, goodAlignReportOut);
951
952                 while (!inputAlignReport.eof()) {               //      need to copy header
953                         char c = inputAlignReport.get();
954                         goodAlignReportOut << c;
955                         if (c == 10 || c == 13){        break;  }       
956                 }
957
958                 while(!inputAlignReport.eof()){
959                         if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); m->mothurRemove(goodAlignReportFile); return 0; }
960
961                         inputAlignReport >> seqName;
962                         it = badSeqNames.find(seqName);
963                         string line;            
964                         while (!inputAlignReport.eof()) {               //      need to copy header
965                                 char c = inputAlignReport.get();
966                                 line += c;
967                                 if (c == 10 || c == 13){        break;  }       
968                         }
969                         
970                         if(it != badSeqNames.end()){
971                                 badSeqNames.erase(it);
972                         }
973                         else{
974                                 goodAlignReportOut << seqName << '\t' << line;
975                         }
976                         m->gobble(inputAlignReport);
977                 }
978                 
979                 if (m->control_pressed) { goodAlignReportOut.close();  inputAlignReport.close(); m->mothurRemove(goodAlignReportFile);  return 0; }
980
981                 //we were unable to remove some of the bad sequences
982                 if (badSeqNames.size() != 0) {
983                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
984                                 m->mothurOut("Your alignreport file does not include the sequence " + *it + " please correct."); 
985                                 m->mothurOutEndLine();
986                         }
987                 }
988
989                 inputAlignReport.close();
990                 goodAlignReportOut.close();
991                                 
992                 if (m->control_pressed) {  m->mothurRemove(goodAlignReportFile);  return 0; }
993                 
994                 return 0;
995         
996         }
997         catch(exception& e) {
998                 m->errorOut(e, "ScreenSeqsCommand", "screenAlignReport");
999                 exit(1);
1000         }
1001         
1002 }
1003 //***************************************************************************************************************
1004
1005 int ScreenSeqsCommand::screenTaxonomy(set<string> badSeqNames){
1006         try {
1007                 ifstream input;
1008                 m->openInputFile(taxonomy, input);
1009                 string seqName, tax;
1010                 set<string>::iterator it;
1011                 
1012                 string goodTaxFile = outputDir + m->getRootName(m->getSimpleName(taxonomy)) + "good" + m->getExtension(taxonomy);
1013                 outputNames.push_back(goodTaxFile);  outputTypes["taxonomy"].push_back(goodTaxFile);
1014                 ofstream goodTaxOut;    m->openOutputFile(goodTaxFile, goodTaxOut);
1015                                 
1016                 while(!input.eof()){
1017                         if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); return 0; }
1018                         
1019                         input >> seqName >> tax;
1020                         it = badSeqNames.find(seqName);
1021                         
1022                         if(it != badSeqNames.end()){ badSeqNames.erase(it); }
1023                         else{
1024                                 goodTaxOut << seqName << '\t' << tax << endl;
1025                         }
1026                         m->gobble(input);
1027                 }
1028                 
1029                 if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); return 0; }
1030                 
1031                 //we were unable to remove some of the bad sequences
1032                 if (badSeqNames.size() != 0) {
1033                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
1034                                 m->mothurOut("Your taxonomy file does not include the sequence " + *it + " please correct."); 
1035                                 m->mothurOutEndLine();
1036                         }
1037                 }
1038                 
1039                 input.close();
1040                 goodTaxOut.close();
1041                 
1042                 if (m->control_pressed) {  m->mothurRemove(goodTaxFile);  return 0; }
1043                 
1044                 return 0;
1045                 
1046         }
1047         catch(exception& e) {
1048                 m->errorOut(e, "ScreenSeqsCommand", "screenTaxonomy");
1049                 exit(1);
1050         }
1051         
1052 }
1053 //***************************************************************************************************************
1054
1055 int ScreenSeqsCommand::screenQual(set<string> badSeqNames){
1056         try {
1057                 ifstream in;
1058                 m->openInputFile(qualfile, in);
1059                 set<string>::iterator it;
1060                 
1061                 string goodQualFile = outputDir + m->getRootName(m->getSimpleName(qualfile)) + "good" + m->getExtension(qualfile);
1062                 outputNames.push_back(goodQualFile);  outputTypes["qfile"].push_back(goodQualFile);
1063                 ofstream goodQual;      m->openOutputFile(goodQualFile, goodQual);
1064                 
1065                 while(!in.eof()){       
1066                         
1067                         if (m->control_pressed) { goodQual.close(); in.close(); m->mothurRemove(goodQualFile); return 0; }
1068
1069                         string saveName = "";
1070                         string name = "";
1071                         string scores = "";
1072                         
1073                         in >> name; 
1074                         
1075                         if (name.length() != 0) { 
1076                                 saveName = name.substr(1);
1077                                 while (!in.eof())       {       
1078                                         char c = in.get(); 
1079                                         if (c == 10 || c == 13){        break;  }
1080                                         else { name += c; }     
1081                                 } 
1082                                 m->gobble(in);
1083                         }
1084                         
1085                         while(in){
1086                                 char letter= in.get();
1087                                 if(letter == '>'){      in.putback(letter);     break;  }
1088                                 else{ scores += letter; }
1089                         }
1090                         
1091                         m->gobble(in);
1092                         
1093                         it = badSeqNames.find(saveName);
1094                         
1095                         if(it != badSeqNames.end()){
1096                                 badSeqNames.erase(it);
1097                         }else{                          
1098                                 goodQual << name << endl << scores;
1099                         }
1100                         
1101                         m->gobble(in);
1102                 }
1103                 
1104                 in.close();
1105                 goodQual.close();
1106                 
1107                 //we were unable to remove some of the bad sequences
1108                 if (badSeqNames.size() != 0) {
1109                         for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
1110                                 m->mothurOut("Your qual file does not include the sequence " + *it + " please correct."); 
1111                                 m->mothurOutEndLine();
1112                         }
1113                 }
1114                 
1115                 if (m->control_pressed) {  m->mothurRemove(goodQualFile);  return 0; }
1116                 
1117                 return 0;
1118                 
1119         }
1120         catch(exception& e) {
1121                 m->errorOut(e, "ScreenSeqsCommand", "screenQual");
1122                 exit(1);
1123         }
1124         
1125 }
1126 //**********************************************************************************************************************
1127
1128 int ScreenSeqsCommand::driver(linePair filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
1129         try {
1130                 ofstream goodFile;
1131                 m->openOutputFile(goodFName, goodFile);
1132                 
1133                 ofstream badAccnosFile;
1134                 m->openOutputFile(badAccnosFName, badAccnosFile);
1135                 
1136                 ifstream inFASTA;
1137                 m->openInputFile(filename, inFASTA);
1138
1139                 inFASTA.seekg(filePos.start);
1140
1141                 bool done = false;
1142                 int count = 0;
1143         
1144                 while (!done) {
1145                 
1146                         if (m->control_pressed) {  return 0; }
1147                         
1148                         Sequence currSeq(inFASTA); m->gobble(inFASTA);
1149                         if (currSeq.getName() != "") {
1150                                 bool goodSeq = 1;               //      innocent until proven guilty
1151                                 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos())                  {       goodSeq = 0;    }
1152                                 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos())                                {       goodSeq = 0;    }
1153                                 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases())                {       goodSeq = 0;    }
1154                                 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer())   {       goodSeq = 0;    }
1155                                 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases())                {       goodSeq = 0;    }
1156                                 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases())                {       goodSeq = 0;    }
1157                                 
1158                                 if(goodSeq == 1){
1159                                         currSeq.printSequence(goodFile);        
1160                                 }
1161                                 else{
1162                                         badAccnosFile << currSeq.getName() << endl;
1163                                         badSeqNames.insert(currSeq.getName());
1164                                 }
1165                         count++;
1166                         }
1167                         
1168                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1169                                 unsigned long long pos = inFASTA.tellg();
1170                                 if ((pos == -1) || (pos >= filePos.end)) { break; }
1171                         #else
1172                                 if (inFASTA.eof()) { break; }
1173                         #endif
1174                         
1175                         //report progress
1176                         if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine();         }
1177                 }
1178                 //report progress
1179                 if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine();         }
1180                 
1181                         
1182                 goodFile.close();
1183                 inFASTA.close();
1184                 badAccnosFile.close();
1185                 
1186                 return count;
1187         }
1188         catch(exception& e) {
1189                 m->errorOut(e, "ScreenSeqsCommand", "driver");
1190                 exit(1);
1191         }
1192 }
1193 //**********************************************************************************************************************
1194 #ifdef USE_MPI
1195 int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badAccnosFile, vector<unsigned long long>& MPIPos, set<string>& badSeqNames){
1196         try {
1197                 string outputString = "";
1198                 MPI_Status statusGood; 
1199                 MPI_Status statusBadAccnos; 
1200                 MPI_Status status; 
1201                 int pid;
1202                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
1203
1204                 for(int i=0;i<num;i++){
1205                 
1206                         if (m->control_pressed) {  return 0; }
1207                         
1208                         //read next sequence
1209                         int length = MPIPos[start+i+1] - MPIPos[start+i];
1210
1211                         char* buf4 = new char[length];
1212                         memcpy(buf4, outputString.c_str(), length);
1213
1214                         MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
1215                         
1216                         string tempBuf = buf4;  delete buf4;
1217                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
1218                         istringstream iss (tempBuf,istringstream::in);
1219                         
1220                         Sequence currSeq(iss);                  
1221                         
1222                         //process seq
1223                         if (currSeq.getName() != "") {
1224                                 bool goodSeq = 1;               //      innocent until proven guilty
1225                                 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos())                  {       goodSeq = 0;    }
1226                                 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos())                                {       goodSeq = 0;    }
1227                                 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases())                {       goodSeq = 0;    }
1228                                 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer())   {       goodSeq = 0;    }
1229                                 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases())                {       goodSeq = 0;    }
1230                                 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases())                {       goodSeq = 0;    }
1231                                 
1232                                 if(goodSeq == 1){
1233                                         outputString =  ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
1234                                 
1235                                         //print good seq
1236                                         length = outputString.length();
1237                                         char* buf2 = new char[length];
1238                                         memcpy(buf2, outputString.c_str(), length);
1239                                         
1240                                         MPI_File_write_shared(goodFile, buf2, length, MPI_CHAR, &statusGood);
1241                                         delete buf2;
1242                                 }
1243                                 else{
1244
1245                                         badSeqNames.insert(currSeq.getName());
1246                                         
1247                                         //write to bad accnos file
1248                                         outputString = currSeq.getName() + "\n";
1249                                 
1250                                         length = outputString.length();
1251                                         char* buf3 = new char[length];
1252                                         memcpy(buf3, outputString.c_str(), length);
1253                                         
1254                                         MPI_File_write_shared(badAccnosFile, buf3, length, MPI_CHAR, &statusBadAccnos);
1255                                         delete buf3;
1256                                 }
1257                         }
1258                         
1259                         //report progress
1260                         if((i) % 100 == 0){     m->mothurOut("Processing sequence: " + toString(i)); m->mothurOutEndLine();             }
1261                 }
1262                                 
1263                 return 1;
1264         }
1265         catch(exception& e) {
1266                 m->errorOut(e, "ScreenSeqsCommand", "driverMPI");
1267                 exit(1);
1268         }
1269 }
1270 #endif
1271 /**************************************************************************************************/
1272
1273 int ScreenSeqsCommand::createProcesses(string goodFileName, string badAccnos, string filename, set<string>& badSeqNames) {
1274         try {
1275         
1276         vector<int> processIDS;   
1277         int process = 1;
1278                 int num = 0;
1279
1280 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1281                                 
1282                 //loop through and create all the processes you want
1283                 while (process != processors) {
1284                         int pid = fork();
1285                         
1286                         if (pid > 0) {
1287                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1288                                 process++;
1289                         }else if (pid == 0){
1290                                 num = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
1291                                 
1292                                 //pass numSeqs to parent
1293                                 ofstream out;
1294                                 string tempFile = filename + toString(getpid()) + ".num.temp";
1295                                 m->openOutputFile(tempFile, out);
1296                                 out << num << endl;
1297                                 out.close();
1298                                 
1299                                 exit(0);
1300                         }else { 
1301                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
1302                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1303                                 exit(0);
1304                         }
1305                 }
1306                 
1307         num = driver(lines[0], goodFileName, badAccnos, filename, badSeqNames);
1308         
1309                 //force parent to wait until all the processes are done
1310                 for (int i=0;i<processIDS.size();i++) { 
1311                         int temp = processIDS[i];
1312                         wait(&temp);
1313                 }
1314                 
1315                 for (int i = 0; i < processIDS.size(); i++) {
1316                         ifstream in;
1317                         string tempFile =  filename + toString(processIDS[i]) + ".num.temp";
1318                         m->openInputFile(tempFile, in);
1319                         if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
1320                         in.close(); m->mothurRemove(tempFile);
1321             
1322             m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
1323             m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
1324                         
1325             m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
1326             m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
1327                 }
1328                 
1329         //read badSeqs in because root process doesnt know what other "bad" seqs the children found
1330         ifstream inBad;
1331         int ableToOpen = m->openInputFile(badAccnos, inBad, "no error");
1332         
1333         if (ableToOpen == 0) {
1334             badSeqNames.clear();
1335             string tempName;
1336             while (!inBad.eof()) {
1337                 inBad >> tempName; m->gobble(inBad);
1338                 badSeqNames.insert(tempName);
1339             }
1340             inBad.close();
1341         }
1342 #else
1343         
1344         //////////////////////////////////////////////////////////////////////////////////////////////////////
1345                 //Windows version shared memory, so be careful when passing variables through the sumScreenData struct. 
1346                 //Above fork() will clone, so memory is separate, but that's not the case with windows, 
1347                 //Taking advantage of shared memory to allow both threads to add info to badSeqNames.
1348                 //////////////////////////////////////////////////////////////////////////////////////////////////////
1349                 
1350                 vector<sumScreenData*> pDataArray; 
1351                 DWORD   dwThreadIdArray[processors-1];
1352                 HANDLE  hThreadArray[processors-1]; 
1353                 
1354                 //Create processor worker threads.
1355                 for( int i=0; i<processors-1; i++ ){
1356             
1357             string extension = "";
1358             if (i!=0) {extension += toString(i) + ".temp"; processIDS.push_back(i); }
1359             
1360                         // Allocate memory for thread data.
1361                         sumScreenData* tempSum = new sumScreenData(startPos, endPos, maxAmbig, maxHomoP, minLength, maxLength, filename, m, lines[i].start, lines[i].end,goodFileName+extension, badAccnos+extension, &badSeqNames);
1362                         pDataArray.push_back(tempSum);
1363                         
1364                         //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
1365                         hThreadArray[i] = CreateThread(NULL, 0, MySumScreenThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);   
1366                 }
1367                 
1368         //do your part
1369         num = driver(lines[processors-1], (goodFileName+toString(processors-1)+".temp"), (badAccnos+toString(processors-1)+".temp"), filename, badSeqNames);
1370         processIDS.push_back(processors-1);
1371         
1372                 //Wait until all threads have terminated.
1373                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
1374                 
1375                 //Close all thread handles and free memory allocations.
1376                 for(int i=0; i < pDataArray.size(); i++){
1377                         num += pDataArray[i]->count;
1378                         CloseHandle(hThreadArray[i]);
1379                         delete pDataArray[i];
1380                 }
1381         
1382         for (int i = 0; i < processIDS.size(); i++) {
1383             m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
1384             m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
1385                         
1386             m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
1387             m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
1388                 }
1389
1390 #endif  
1391         
1392         return num;
1393         
1394         }
1395         catch(exception& e) {
1396                 m->errorOut(e, "ScreenSeqsCommand", "createProcesses");
1397                 exit(1);
1398         }
1399 }
1400
1401 //***************************************************************************************************************
1402
1403