2 * screenseqscommand.cpp
5 * Created by Pat Schloss on 6/3/09.
6 * Copyright 2009 Patrick D. Schloss. All rights reserved.
10 #include "screenseqscommand.h"
13 //**********************************************************************************************************************
14 vector<string> ScreenSeqsCommand::setParameters(){
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);
34 vector<string> myArray;
35 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
39 m->errorOut(e, "ScreenSeqsCommand", "setParameters");
43 //**********************************************************************************************************************
44 string ScreenSeqsCommand::getHelpString(){
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";
69 m->errorOut(e, "ScreenSeqsCommand", "getHelpString");
73 //**********************************************************************************************************************
74 ScreenSeqsCommand::ScreenSeqsCommand(){
76 abort = true; calledHelp = true;
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;
88 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
92 //***************************************************************************************************************
94 ScreenSeqsCommand::ScreenSeqsCommand(string option) {
96 abort = false; calledHelp = false;
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;}
103 vector<string> myArray = setParameters();
105 OptionParser parser(option);
106 map<string,string> parameters = parser.getParameters();
108 ValidParameters validParameter("screen.seqs");
109 map<string,string>::iterator it;
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; }
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;
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 = ""; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
187 else if (fastafile == "not open") { abort = true; }
188 else { m->setFastaFile(fastafile); }
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); }
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); }
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); }
205 alignreport = validParameter.validFile(parameters, "alignreport", true);
206 if (alignreport == "not open") { abort = true; }
207 else if (alignreport == "not found") { alignreport = ""; }
209 taxonomy = validParameter.validFile(parameters, "taxonomy", true);
210 if (taxonomy == "not open") { abort = true; }
211 else if (taxonomy == "not found") { taxonomy = ""; }
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"){
216 outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it
219 //check for optional parameter and set defaults
220 // ...at some point should added some additional type checking...
222 temp = validParameter.validFile(parameters, "start", false); if (temp == "not found") { temp = "-1"; }
223 m->mothurConvert(temp, startPos);
225 temp = validParameter.validFile(parameters, "end", false); if (temp == "not found") { temp = "-1"; }
226 m->mothurConvert(temp, endPos);
228 temp = validParameter.validFile(parameters, "maxambig", false); if (temp == "not found") { temp = "-1"; }
229 m->mothurConvert(temp, maxAmbig);
231 temp = validParameter.validFile(parameters, "maxhomop", false); if (temp == "not found") { temp = "-1"; }
232 m->mothurConvert(temp, maxHomoP);
234 temp = validParameter.validFile(parameters, "minlength", false); if (temp == "not found") { temp = "-1"; }
235 m->mothurConvert(temp, minLength);
237 temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found") { temp = "-1"; }
238 m->mothurConvert(temp, maxLength);
240 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
241 m->setProcessors(temp);
242 m->mothurConvert(temp, processors);
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);
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);
259 if (optimize.size() == 1) { if (optimize[0] == "none") { optimize.clear(); } }
261 temp = validParameter.validFile(parameters, "criteria", false); if (temp == "not found"){ temp = "90"; }
262 m->mothurConvert(temp, criteria);
264 if (namefile == "") {
265 vector<string> files; files.push_back(fastafile);
266 parser.getNameFile(files);
271 catch(exception& e) {
272 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
276 //***************************************************************************************************************
278 int ScreenSeqsCommand::execute(){
281 if (abort == true) { if (calledHelp) { return 0; } return 2; }
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);
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)])); }
295 if(processors == 1){ lines.push_back(linePair(0, 1000)); }
297 int numFastaSeqs = 0;
298 positions = m->setFilePosFasta(fastafile, numFastaSeqs);
299 if (positions.size() < processors) { processors = positions.size(); }
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));
312 string goodSeqFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "good" + m->getExtension(fastafile);
313 string badAccnosFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "bad.accnos";
315 int numFastaSeqs = 0;
316 set<string> badSeqNames;
317 int start = time(NULL);
320 int pid, numSeqsPerProcessor;
322 vector<unsigned long long> MPIPos;
325 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
326 MPI_Comm_size(MPI_COMM_WORLD, &processors);
330 MPI_File outMPIBadAccnos;
332 int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
333 int inMode=MPI_MODE_RDONLY;
335 char outGoodFilename[1024];
336 strcpy(outGoodFilename, goodSeqFile.c_str());
338 char outBadAccnosFilename[1024];
339 strcpy(outBadAccnosFilename, badAccnosFile.c_str());
341 char inFileName[1024];
342 strcpy(inFileName, fastafile.c_str());
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);
348 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
350 if (pid == 0) { //you are the root process
352 MPIPos = m->setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
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);
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; }
366 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
368 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
370 for (int i = 1; i < processors; i++) {
373 MPI_Recv(&badSize, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
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);
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; }
386 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
388 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
391 int badSize = badSeqNames.size();
392 MPI_Send(&badSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
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
402 if(processors == 1){ numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames); }
403 else{ numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames); }
405 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
409 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
411 if (pid == 0) { //only one process should fix files
413 //read accnos file with all names in it, process 0 just has its names
414 MPI_File inMPIAccnos;
417 char inFileName[1024];
418 strcpy(inFileName, badAccnosFile.c_str());
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);
423 char* buffer = new char[size];
424 MPI_File_read(inMPIAccnos, buffer, size, MPI_CHAR, &status);
426 string tempBuf = buffer;
427 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size); }
428 istringstream iss (tempBuf,istringstream::in);
431 MPI_File_close(&inMPIAccnos);
436 iss >> tempName; m->gobble(iss);
437 badSeqNames.insert(tempName);
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
449 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
451 if(alignreport != "") { screenAlignReport(badSeqNames); }
452 if(qualfile != "") { screenQual(badSeqNames); }
453 if(taxonomy != "") { screenTaxonomy(badSeqNames); }
455 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
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();
469 //set fasta file as new current fastafile
471 itTypes = outputTypes.find("fasta");
472 if (itTypes != outputTypes.end()) {
473 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
476 itTypes = outputTypes.find("name");
477 if (itTypes != outputTypes.end()) {
478 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
481 itTypes = outputTypes.find("group");
482 if (itTypes != outputTypes.end()) {
483 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
486 itTypes = outputTypes.find("qfile");
487 if (itTypes != outputTypes.end()) {
488 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setQualFile(current); }
491 itTypes = outputTypes.find("taxonomy");
492 if (itTypes != outputTypes.end()) {
493 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
496 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
497 m->mothurOutEndLine();
501 catch(exception& e) {
502 m->errorOut(e, "ScreenSeqsCommand", "execute");
507 //***************************************************************************************************************
509 int ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
512 m->openInputFile(namefile, inputNames);
513 set<string> badSeqGroups;
514 string seqName, seqList, group;
515 set<string>::iterator it;
517 string goodNameFile = outputDir + m->getRootName(m->getSimpleName(namefile)) + "good" + m->getExtension(namefile);
518 outputNames.push_back(goodNameFile); outputTypes["name"].push_back(goodNameFile);
520 ofstream goodNameOut; m->openOutputFile(goodNameFile, goodNameOut);
522 while(!inputNames.eof()){
523 if (m->control_pressed) { goodNameOut.close(); inputNames.close(); m->mothurRemove(goodNameFile); return 0; }
525 inputNames >> seqName >> seqList;
526 it = badSeqNames.find(seqName);
528 if(it != badSeqNames.end()){
529 badSeqNames.erase(it);
533 for(int i=0;i<seqList.length();i++){
534 if(seqList[i] == ','){
535 badSeqGroups.insert(seqList.substr(start,i-start));
539 badSeqGroups.insert(seqList.substr(start,seqList.length()-start));
543 goodNameOut << seqName << '\t' << seqList << endl;
545 m->gobble(inputNames);
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();
560 ifstream inputGroups;
561 m->openInputFile(groupfile, inputGroups);
563 string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
564 outputNames.push_back(goodGroupFile); outputTypes["group"].push_back(goodGroupFile);
566 ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
568 while(!inputGroups.eof()){
569 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodNameFile); m->mothurRemove(goodGroupFile); return 0; }
571 inputGroups >> seqName >> group;
573 it = badSeqGroups.find(seqName);
575 if(it != badSeqGroups.end()){
576 badSeqGroups.erase(it);
579 goodGroupOut << seqName << '\t' << group << endl;
581 m->gobble(inputGroups);
584 goodGroupOut.close();
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();
599 catch(exception& e) {
600 m->errorOut(e, "ScreenSeqsCommand", "screenNameGroupFile");
604 //***************************************************************************************************************
605 int ScreenSeqsCommand::getSummary(vector<unsigned long long>& positions){
608 vector<int> startPosition;
609 vector<int> endPosition;
610 vector<int> seqLength;
611 vector<int> ambigBases;
612 vector<int> longHomoPolymer;
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)])); }
619 if(processors == 1){ lines.push_back(linePair(0, 1000)); }
621 int numFastaSeqs = 0;
622 positions = m->setFilePosFasta(fastafile, numFastaSeqs);
623 if (positions.size() < processors) { processors = positions.size(); }
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));
637 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
640 driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
643 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
645 numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
647 numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile);
650 if (m->control_pressed) { return 0; }
652 // numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
653 // if (m->control_pressed) { return 0; }
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());
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));
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(); }
678 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
679 MPI_Comm_size(MPI_COMM_WORLD, &processors);
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);
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);
699 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
703 catch(exception& e) {
704 m->errorOut(e, "ScreenSeqsCommand", "getSummary");
708 /**************************************************************************************/
709 int ScreenSeqsCommand::driverCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, linePair filePos) {
713 m->openInputFile(filename, in);
715 in.seekg(filePos.start);
722 if (m->control_pressed) { in.close(); return 1; }
724 Sequence current(in); m->gobble(in);
726 if (current.getName() != "") {
728 if (namefile != "") {
729 //make sure this sequence is in the namefile, else error
730 map<string, int>::iterator it = nameMap.find(current.getName());
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; }
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());
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; }
752 if (in.eof()) { break; }
761 catch(exception& e) {
762 m->errorOut(e, "ScreenSeqsCommand", "driverCreateSummary");
766 /**************************************************************************************************/
767 int ScreenSeqsCommand::createProcessesCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename) {
772 vector<int> processIDS;
774 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
776 //loop through and create all the processes you want
777 while (process != processors) {
781 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
784 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[process]);
786 //pass numSeqs to parent
788 string tempFile = fastafile + toString(getpid()) + ".num.temp";
789 m->openOutputFile(tempFile, out);
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;
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); }
809 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
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];
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";
821 m->openInputFile(tempFilename, in);
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);
833 m->mothurRemove(tempFilename);
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 //////////////////////////////////////////////////////////////////////////////////////////////////////
844 vector<sumData*> pDataArray;
845 DWORD dwThreadIdArray[processors-1];
846 HANDLE hThreadArray[processors-1];
848 //Create processor worker threads.
849 for( int i=0; i<processors-1; i++ ){
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);
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]);
861 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[processors-1]);
863 //Wait until all threads have terminated.
864 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
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];
881 catch(exception& e) {
882 m->errorOut(e, "ScreenSeqsCommand", "createProcessesCreateSummary");
887 //***************************************************************************************************************
889 int ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
891 ifstream inputGroups;
892 m->openInputFile(groupfile, inputGroups);
893 string seqName, group;
894 set<string>::iterator it;
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);
900 while(!inputGroups.eof()){
901 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodGroupFile); return 0; }
903 inputGroups >> seqName >> group;
904 it = badSeqNames.find(seqName);
906 if(it != badSeqNames.end()){
907 badSeqNames.erase(it);
910 goodGroupOut << seqName << '\t' << group << endl;
912 m->gobble(inputGroups);
915 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodGroupFile); return 0; }
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();
926 goodGroupOut.close();
928 if (m->control_pressed) { m->mothurRemove(goodGroupFile); }
933 catch(exception& e) {
934 m->errorOut(e, "ScreenSeqsCommand", "screenGroupFile");
939 //***************************************************************************************************************
941 int ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
943 ifstream inputAlignReport;
944 m->openInputFile(alignreport, inputAlignReport);
945 string seqName, group;
946 set<string>::iterator it;
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);
952 while (!inputAlignReport.eof()) { // need to copy header
953 char c = inputAlignReport.get();
954 goodAlignReportOut << c;
955 if (c == 10 || c == 13){ break; }
958 while(!inputAlignReport.eof()){
959 if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); m->mothurRemove(goodAlignReportFile); return 0; }
961 inputAlignReport >> seqName;
962 it = badSeqNames.find(seqName);
964 while (!inputAlignReport.eof()) { // need to copy header
965 char c = inputAlignReport.get();
967 if (c == 10 || c == 13){ break; }
970 if(it != badSeqNames.end()){
971 badSeqNames.erase(it);
974 goodAlignReportOut << seqName << '\t' << line;
976 m->gobble(inputAlignReport);
979 if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); m->mothurRemove(goodAlignReportFile); return 0; }
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();
989 inputAlignReport.close();
990 goodAlignReportOut.close();
992 if (m->control_pressed) { m->mothurRemove(goodAlignReportFile); return 0; }
997 catch(exception& e) {
998 m->errorOut(e, "ScreenSeqsCommand", "screenAlignReport");
1003 //***************************************************************************************************************
1005 int ScreenSeqsCommand::screenTaxonomy(set<string> badSeqNames){
1008 m->openInputFile(taxonomy, input);
1009 string seqName, tax;
1010 set<string>::iterator it;
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);
1016 while(!input.eof()){
1017 if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); return 0; }
1019 input >> seqName >> tax;
1020 it = badSeqNames.find(seqName);
1022 if(it != badSeqNames.end()){ badSeqNames.erase(it); }
1024 goodTaxOut << seqName << '\t' << tax << endl;
1029 if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); return 0; }
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();
1042 if (m->control_pressed) { m->mothurRemove(goodTaxFile); return 0; }
1047 catch(exception& e) {
1048 m->errorOut(e, "ScreenSeqsCommand", "screenTaxonomy");
1053 //***************************************************************************************************************
1055 int ScreenSeqsCommand::screenQual(set<string> badSeqNames){
1058 m->openInputFile(qualfile, in);
1059 set<string>::iterator it;
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);
1067 if (m->control_pressed) { goodQual.close(); in.close(); m->mothurRemove(goodQualFile); return 0; }
1069 string saveName = "";
1075 if (name.length() != 0) {
1076 saveName = name.substr(1);
1079 if (c == 10 || c == 13){ break; }
1086 char letter= in.get();
1087 if(letter == '>'){ in.putback(letter); break; }
1088 else{ scores += letter; }
1093 it = badSeqNames.find(saveName);
1095 if(it != badSeqNames.end()){
1096 badSeqNames.erase(it);
1098 goodQual << name << endl << scores;
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();
1115 if (m->control_pressed) { m->mothurRemove(goodQualFile); return 0; }
1120 catch(exception& e) {
1121 m->errorOut(e, "ScreenSeqsCommand", "screenQual");
1126 //**********************************************************************************************************************
1128 int ScreenSeqsCommand::driver(linePair filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
1131 m->openOutputFile(goodFName, goodFile);
1133 ofstream badAccnosFile;
1134 m->openOutputFile(badAccnosFName, badAccnosFile);
1137 m->openInputFile(filename, inFASTA);
1139 inFASTA.seekg(filePos.start);
1146 if (m->control_pressed) { return 0; }
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; }
1159 currSeq.printSequence(goodFile);
1162 badAccnosFile << currSeq.getName() << endl;
1163 badSeqNames.insert(currSeq.getName());
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; }
1172 if (inFASTA.eof()) { break; }
1176 if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
1179 if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
1184 badAccnosFile.close();
1188 catch(exception& e) {
1189 m->errorOut(e, "ScreenSeqsCommand", "driver");
1193 //**********************************************************************************************************************
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){
1197 string outputString = "";
1198 MPI_Status statusGood;
1199 MPI_Status statusBadAccnos;
1202 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
1204 for(int i=0;i<num;i++){
1206 if (m->control_pressed) { return 0; }
1208 //read next sequence
1209 int length = MPIPos[start+i+1] - MPIPos[start+i];
1211 char* buf4 = new char[length];
1212 memcpy(buf4, outputString.c_str(), length);
1214 MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
1216 string tempBuf = buf4; delete buf4;
1217 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
1218 istringstream iss (tempBuf,istringstream::in);
1220 Sequence currSeq(iss);
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; }
1233 outputString = ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
1236 length = outputString.length();
1237 char* buf2 = new char[length];
1238 memcpy(buf2, outputString.c_str(), length);
1240 MPI_File_write_shared(goodFile, buf2, length, MPI_CHAR, &statusGood);
1245 badSeqNames.insert(currSeq.getName());
1247 //write to bad accnos file
1248 outputString = currSeq.getName() + "\n";
1250 length = outputString.length();
1251 char* buf3 = new char[length];
1252 memcpy(buf3, outputString.c_str(), length);
1254 MPI_File_write_shared(badAccnosFile, buf3, length, MPI_CHAR, &statusBadAccnos);
1260 if((i) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(i)); m->mothurOutEndLine(); }
1265 catch(exception& e) {
1266 m->errorOut(e, "ScreenSeqsCommand", "driverMPI");
1271 /**************************************************************************************************/
1273 int ScreenSeqsCommand::createProcesses(string goodFileName, string badAccnos, string filename, set<string>& badSeqNames) {
1276 vector<int> processIDS;
1280 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1282 //loop through and create all the processes you want
1283 while (process != processors) {
1287 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1289 }else if (pid == 0){
1290 num = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
1292 //pass numSeqs to parent
1294 string tempFile = filename + toString(getpid()) + ".num.temp";
1295 m->openOutputFile(tempFile, out);
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); }
1307 num = driver(lines[0], goodFileName, badAccnos, filename, badSeqNames);
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];
1315 for (int i = 0; i < processIDS.size(); i++) {
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);
1322 m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
1323 m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
1325 m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
1326 m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
1329 //read badSeqs in because root process doesnt know what other "bad" seqs the children found
1331 int ableToOpen = m->openInputFile(badAccnos, inBad, "no error");
1333 if (ableToOpen == 0) {
1334 badSeqNames.clear();
1336 while (!inBad.eof()) {
1337 inBad >> tempName; m->gobble(inBad);
1338 badSeqNames.insert(tempName);
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 //////////////////////////////////////////////////////////////////////////////////////////////////////
1350 vector<sumScreenData*> pDataArray;
1351 DWORD dwThreadIdArray[processors-1];
1352 HANDLE hThreadArray[processors-1];
1354 //Create processor worker threads.
1355 for( int i=0; i<processors-1; i++ ){
1357 string extension = "";
1358 if (i!=0) {extension += toString(i) + ".temp"; processIDS.push_back(i); }
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);
1362 pDataArray.push_back(tempSum);
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]);
1369 num = driver(lines[processors-1], (goodFileName+toString(processors-1)+".temp"), (badAccnos+toString(processors-1)+".temp"), filename, badSeqNames);
1370 processIDS.push_back(processors-1);
1372 //Wait until all threads have terminated.
1373 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
1375 //Close all thread handles and free memory allocations.
1376 for(int i=0; i < pDataArray.size(); i++){
1377 num += pDataArray[i]->count;
1378 for (set<string>::iterator it = pDataArray[i]->badSeqNames.begin(); it != pDataArray[i]->badSeqNames.end(); it++) { badSeqNames.insert(*it); }
1379 CloseHandle(hThreadArray[i]);
1380 delete pDataArray[i];
1383 for (int i = 0; i < processIDS.size(); i++) {
1384 m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
1385 m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
1387 m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
1388 m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
1396 catch(exception& e) {
1397 m->errorOut(e, "ScreenSeqsCommand", "createProcesses");
1402 //***************************************************************************************************************