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"
11 #include "counttable.h"
13 //**********************************************************************************************************************
14 vector<string> ScreenSeqsCommand::setParameters(){
16 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","fasta",false,true,true); parameters.push_back(pfasta);
17 CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","name",false,false,true); parameters.push_back(pname);
18 CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","count",false,false,true); parameters.push_back(pcount);
19 CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","group",false,false,true); parameters.push_back(pgroup);
20 CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "none","qfile",false,false); parameters.push_back(pqfile);
21 CommandParameter palignreport("alignreport", "InputTypes", "", "", "none", "none", "none","alignreport",false,false); parameters.push_back(palignreport);
22 CommandParameter ptax("taxonomy", "InputTypes", "", "", "none", "none", "none","taxonomy",false,false); parameters.push_back(ptax);
23 CommandParameter pstart("start", "Number", "", "-1", "", "", "","",false,false,true); parameters.push_back(pstart);
24 CommandParameter pend("end", "Number", "", "-1", "", "", "","",false,false,true); parameters.push_back(pend);
25 CommandParameter pmaxambig("maxambig", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pmaxambig);
26 CommandParameter pmaxhomop("maxhomop", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pmaxhomop);
27 CommandParameter pminlength("minlength", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pminlength);
28 CommandParameter pmaxlength("maxlength", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pmaxlength);
29 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
30 CommandParameter pcriteria("criteria", "Number", "", "90", "", "", "","",false,false); parameters.push_back(pcriteria);
31 CommandParameter poptimize("optimize", "Multiple", "none-start-end-maxambig-maxhomop-minlength-maxlength", "none", "", "", "","",true,false); parameters.push_back(poptimize);
32 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
33 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
35 vector<string> myArray;
36 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
40 m->errorOut(e, "ScreenSeqsCommand", "setParameters");
44 //**********************************************************************************************************************
45 string ScreenSeqsCommand::getHelpString(){
47 string helpString = "";
48 helpString += "The screen.seqs command reads a fastafile and screens sequences.\n";
49 helpString += "The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group, count, qfile, alignreport, taxonomy, optimize, criteria and processors.\n";
50 helpString += "The fasta parameter is required.\n";
51 helpString += "The alignreport and taxonomy parameters allow you to remove bad seqs from taxonomy and alignreport files.\n";
52 helpString += "The start parameter is used to set a position the \"good\" sequences must start by. The default is -1.\n";
53 helpString += "The end parameter is used to set a position the \"good\" sequences must end after. The default is -1.\n";
54 helpString += "The maxambig parameter allows you to set the maximum number of ambigious bases allowed. The default is -1.\n";
55 helpString += "The maxhomop parameter allows you to set a maximum homopolymer length. \n";
56 helpString += "The minlength parameter allows you to set and minimum sequence length. \n";
57 helpString += "The maxlength parameter allows you to set and maximum sequence length. \n";
58 helpString += "The processors parameter allows you to specify the number of processors to use while running the command. The default is 1.\n";
59 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";
60 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";
61 helpString += "The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n";
62 helpString += "The screen.seqs command should be in the following format: \n";
63 helpString += "screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig, \n";
64 helpString += "maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n";
65 helpString += "Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n";
66 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
70 m->errorOut(e, "ScreenSeqsCommand", "getHelpString");
74 //**********************************************************************************************************************
75 string ScreenSeqsCommand::getOutputPattern(string type) {
79 if (type == "fasta") { pattern = "[filename],good,[extension]"; }
80 else if (type == "taxonomy") { pattern = "[filename],good,[extension]"; }
81 else if (type == "name") { pattern = "[filename],good,[extension]"; }
82 else if (type == "group") { pattern = "[filename],good,[extension]"; }
83 else if (type == "count") { pattern = "[filename],good,[extension]"; }
84 else if (type == "accnos") { pattern = "[filename],bad.accnos"; }
85 else if (type == "qfile") { pattern = "[filename],good,[extension]"; }
86 else if (type == "alignreport") { pattern = "[filename],good.align.report"; }
87 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
92 m->errorOut(e, "ScreenSeqsCommand", "getOutputPattern");
96 //**********************************************************************************************************************
97 ScreenSeqsCommand::ScreenSeqsCommand(){
99 abort = true; calledHelp = true;
101 vector<string> tempOutNames;
102 outputTypes["fasta"] = tempOutNames;
103 outputTypes["name"] = tempOutNames;
104 outputTypes["group"] = tempOutNames;
105 outputTypes["alignreport"] = tempOutNames;
106 outputTypes["accnos"] = tempOutNames;
107 outputTypes["qfile"] = tempOutNames;
108 outputTypes["taxonomy"] = tempOutNames;
109 outputTypes["count"] = tempOutNames;
111 catch(exception& e) {
112 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
116 //***************************************************************************************************************
118 ScreenSeqsCommand::ScreenSeqsCommand(string option) {
120 abort = false; calledHelp = false;
122 //allow user to run help
123 if(option == "help") { help(); abort = true; calledHelp = true; }
124 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
127 vector<string> myArray = setParameters();
129 OptionParser parser(option);
130 map<string,string> parameters = parser.getParameters();
132 ValidParameters validParameter("screen.seqs");
133 map<string,string>::iterator it;
135 //check to make sure all parameters are valid for command
136 for (it = parameters.begin(); it != parameters.end(); it++) {
137 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
140 //initialize outputTypes
141 vector<string> tempOutNames;
142 outputTypes["fasta"] = tempOutNames;
143 outputTypes["name"] = tempOutNames;
144 outputTypes["group"] = tempOutNames;
145 outputTypes["alignreport"] = tempOutNames;
146 outputTypes["accnos"] = tempOutNames;
147 outputTypes["qfile"] = tempOutNames;
148 outputTypes["taxonomy"] = tempOutNames;
149 outputTypes["count"] = tempOutNames;
151 //if the user changes the input directory command factory will send this info to us in the output parameter
152 string inputDir = validParameter.validFile(parameters, "inputdir", false);
153 if (inputDir == "not found"){ inputDir = ""; }
156 it = parameters.find("fasta");
157 //user has given a template file
158 if(it != parameters.end()){
159 path = m->hasPath(it->second);
160 //if the user has not given a path then, add inputdir. else leave path alone.
161 if (path == "") { parameters["fasta"] = inputDir + it->second; }
164 it = parameters.find("group");
165 //user has given a template file
166 if(it != parameters.end()){
167 path = m->hasPath(it->second);
168 //if the user has not given a path then, add inputdir. else leave path alone.
169 if (path == "") { parameters["group"] = inputDir + it->second; }
172 it = parameters.find("name");
173 //user has given a template file
174 if(it != parameters.end()){
175 path = m->hasPath(it->second);
176 //if the user has not given a path then, add inputdir. else leave path alone.
177 if (path == "") { parameters["name"] = inputDir + it->second; }
180 it = parameters.find("alignreport");
181 //user has given a template file
182 if(it != parameters.end()){
183 path = m->hasPath(it->second);
184 //if the user has not given a path then, add inputdir. else leave path alone.
185 if (path == "") { parameters["alignreport"] = inputDir + it->second; }
188 it = parameters.find("qfile");
189 //user has given a template file
190 if(it != parameters.end()){
191 path = m->hasPath(it->second);
192 //if the user has not given a path then, add inputdir. else leave path alone.
193 if (path == "") { parameters["qfile"] = inputDir + it->second; }
196 it = parameters.find("taxonomy");
197 //user has given a template file
198 if(it != parameters.end()){
199 path = m->hasPath(it->second);
200 //if the user has not given a path then, add inputdir. else leave path alone.
201 if (path == "") { parameters["taxonomy"] = inputDir + it->second; }
204 it = parameters.find("count");
205 //user has given a template file
206 if(it != parameters.end()){
207 path = m->hasPath(it->second);
208 //if the user has not given a path then, add inputdir. else leave path alone.
209 if (path == "") { parameters["count"] = inputDir + it->second; }
213 //check for required parameters
214 fastafile = validParameter.validFile(parameters, "fasta", true);
215 if (fastafile == "not found") {
216 fastafile = m->getFastaFile();
217 if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
218 else { m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
220 else if (fastafile == "not open") { abort = true; }
221 else { m->setFastaFile(fastafile); }
223 groupfile = validParameter.validFile(parameters, "group", true);
224 if (groupfile == "not open") { abort = true; }
225 else if (groupfile == "not found") { groupfile = ""; }
226 else { m->setGroupFile(groupfile); }
228 qualfile = validParameter.validFile(parameters, "qfile", true);
229 if (qualfile == "not open") { abort = true; }
230 else if (qualfile == "not found") { qualfile = ""; }
231 else { m->setQualFile(qualfile); }
233 namefile = validParameter.validFile(parameters, "name", true);
234 if (namefile == "not open") { namefile = ""; abort = true; }
235 else if (namefile == "not found") { namefile = ""; }
236 else { m->setNameFile(namefile); }
238 countfile = validParameter.validFile(parameters, "count", true);
239 if (countfile == "not open") { countfile = ""; abort = true; }
240 else if (countfile == "not found") { countfile = ""; }
241 else { m->setCountTableFile(countfile); }
243 if ((namefile != "") && (countfile != "")) {
244 m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
247 if ((groupfile != "") && (countfile != "")) {
248 m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true;
251 alignreport = validParameter.validFile(parameters, "alignreport", true);
252 if (alignreport == "not open") { abort = true; }
253 else if (alignreport == "not found") { alignreport = ""; }
255 taxonomy = validParameter.validFile(parameters, "taxonomy", true);
256 if (taxonomy == "not open") { abort = true; }
257 else if (taxonomy == "not found") { taxonomy = ""; }
259 //if the user changes the output directory command factory will send this info to us in the output parameter
260 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
262 outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it
265 //check for optional parameter and set defaults
266 // ...at some point should added some additional type checking...
268 temp = validParameter.validFile(parameters, "start", false); if (temp == "not found") { temp = "-1"; }
269 m->mothurConvert(temp, startPos);
271 temp = validParameter.validFile(parameters, "end", false); if (temp == "not found") { temp = "-1"; }
272 m->mothurConvert(temp, endPos);
274 temp = validParameter.validFile(parameters, "maxambig", false); if (temp == "not found") { temp = "-1"; }
275 m->mothurConvert(temp, maxAmbig);
277 temp = validParameter.validFile(parameters, "maxhomop", false); if (temp == "not found") { temp = "-1"; }
278 m->mothurConvert(temp, maxHomoP);
280 temp = validParameter.validFile(parameters, "minlength", false); if (temp == "not found") { temp = "-1"; }
281 m->mothurConvert(temp, minLength);
283 temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found") { temp = "-1"; }
284 m->mothurConvert(temp, maxLength);
286 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
287 m->setProcessors(temp);
288 m->mothurConvert(temp, processors);
290 temp = validParameter.validFile(parameters, "optimize", false); //optimizing trumps the optimized values original value
291 if (temp == "not found"){ temp = "none"; }
292 m->splitAtDash(temp, optimize);
294 //check for invalid optimize options
295 set<string> validOptimizers;
296 validOptimizers.insert("none"); validOptimizers.insert("start"); validOptimizers.insert("end"); validOptimizers.insert("maxambig"); validOptimizers.insert("maxhomop"); validOptimizers.insert("minlength"); validOptimizers.insert("maxlength");
297 for (int i = 0; i < optimize.size(); i++) {
298 if (validOptimizers.count(optimize[i]) == 0) {
299 m->mothurOut(optimize[i] + " is not a valid optimizer. Valid options are start, end, maxambig, maxhomop, minlength and maxlength."); m->mothurOutEndLine();
300 optimize.erase(optimize.begin()+i);
305 if (optimize.size() == 1) { if (optimize[0] == "none") { optimize.clear(); } }
307 temp = validParameter.validFile(parameters, "criteria", false); if (temp == "not found"){ temp = "90"; }
308 m->mothurConvert(temp, criteria);
310 if (countfile == "") {
311 if (namefile == "") {
312 vector<string> files; files.push_back(fastafile);
313 parser.getNameFile(files);
319 catch(exception& e) {
320 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
324 //***************************************************************************************************************
326 int ScreenSeqsCommand::execute(){
329 if (abort == true) { if (calledHelp) { return 0; } return 2; }
331 //if the user want to optimize we need to know the 90% mark
332 vector<unsigned long long> positions;
333 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
334 //use the namefile to optimize correctly
335 if (namefile != "") { nameMap = m->readNames(namefile); }
336 else if (countfile != "") {
338 ct.readTable(countfile);
339 nameMap = ct.getNameMap();
341 getSummary(positions);
344 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
345 positions = m->divideFile(fastafile, processors);
346 for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }
348 if(processors == 1){ lines.push_back(linePair(0, 1000)); }
350 int numFastaSeqs = 0;
351 positions = m->setFilePosFasta(fastafile, numFastaSeqs);
352 if (positions.size() < processors) { processors = positions.size(); }
354 //figure out how many sequences you have to process
355 int numSeqsPerProcessor = numFastaSeqs / processors;
356 for (int i = 0; i < processors; i++) {
357 int startIndex = i * numSeqsPerProcessor;
358 if(i == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor; }
359 lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
365 map<string, string> variables;
366 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(fastafile));
367 string badAccnosFile = getOutputFileName("accnos",variables);
368 variables["[extension]"] = m->getExtension(fastafile);
369 string goodSeqFile = getOutputFileName("fasta", variables);
372 int numFastaSeqs = 0;
373 set<string> badSeqNames;
374 int start = time(NULL);
377 int pid, numSeqsPerProcessor;
379 vector<unsigned long long> MPIPos;
382 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
383 MPI_Comm_size(MPI_COMM_WORLD, &processors);
387 MPI_File outMPIBadAccnos;
389 int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
390 int inMode=MPI_MODE_RDONLY;
392 char outGoodFilename[1024];
393 strcpy(outGoodFilename, goodSeqFile.c_str());
395 char outBadAccnosFilename[1024];
396 strcpy(outBadAccnosFilename, badAccnosFile.c_str());
398 char inFileName[1024];
399 strcpy(inFileName, fastafile.c_str());
401 MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI); //comm, filename, mode, info, filepointer
402 MPI_File_open(MPI_COMM_WORLD, outGoodFilename, outMode, MPI_INFO_NULL, &outMPIGood);
403 MPI_File_open(MPI_COMM_WORLD, outBadAccnosFilename, outMode, MPI_INFO_NULL, &outMPIBadAccnos);
405 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
407 if (pid == 0) { //you are the root process
409 MPIPos = m->setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
411 //send file positions to all processes
412 for(int i = 1; i < processors; i++) {
413 MPI_Send(&numFastaSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
414 MPI_Send(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
417 //figure out how many sequences you have to align
418 numSeqsPerProcessor = numFastaSeqs / processors;
419 int startIndex = pid * numSeqsPerProcessor;
420 if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
423 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
425 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
427 for (int i = 1; i < processors; i++) {
430 MPI_Recv(&badSize, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
432 }else{ //you are a child process
433 MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
434 MPIPos.resize(numFastaSeqs+1);
435 MPI_Recv(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
437 //figure out how many sequences you have to align
438 numSeqsPerProcessor = numFastaSeqs / processors;
439 int startIndex = pid * numSeqsPerProcessor;
440 if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
443 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
445 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
448 int badSize = badSeqNames.size();
449 MPI_Send(&badSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
453 MPI_File_close(&inMPI);
454 MPI_File_close(&outMPIGood);
455 MPI_File_close(&outMPIBadAccnos);
456 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
459 if(processors == 1){ numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames); }
460 else{ numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames); }
462 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
466 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
468 if (pid == 0) { //only one process should fix files
470 //read accnos file with all names in it, process 0 just has its names
471 MPI_File inMPIAccnos;
474 char inFileName[1024];
475 strcpy(inFileName, badAccnosFile.c_str());
477 MPI_File_open(MPI_COMM_SELF, inFileName, inMode, MPI_INFO_NULL, &inMPIAccnos); //comm, filename, mode, info, filepointer
478 MPI_File_get_size(inMPIAccnos, &size);
480 char* buffer = new char[size];
481 MPI_File_read(inMPIAccnos, buffer, size, MPI_CHAR, &status);
483 string tempBuf = buffer;
484 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size); }
485 istringstream iss (tempBuf,istringstream::in);
488 MPI_File_close(&inMPIAccnos);
493 iss >> tempName; m->gobble(iss);
494 badSeqNames.insert(tempName);
498 if(namefile != "" && groupfile != "") {
499 screenNameGroupFile(badSeqNames);
500 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
501 }else if(namefile != "") {
502 screenNameGroupFile(badSeqNames);
503 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
504 }else if(groupfile != "") { screenGroupFile(badSeqNames); } // this screens just the group
505 else if (countfile != "") { screenCountFile(badSeqNames); }
508 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
510 if(alignreport != "") { screenAlignReport(badSeqNames); }
511 if(qualfile != "") { screenQual(badSeqNames); }
512 if(taxonomy != "") { screenTaxonomy(badSeqNames); }
514 if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
520 m->mothurOutEndLine();
521 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
522 m->mothurOut(goodSeqFile); m->mothurOutEndLine(); outputTypes["fasta"].push_back(goodSeqFile);
523 m->mothurOut(badAccnosFile); m->mothurOutEndLine(); outputTypes["accnos"].push_back(badAccnosFile);
524 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
525 m->mothurOutEndLine();
526 m->mothurOutEndLine();
528 //set fasta file as new current fastafile
530 itTypes = outputTypes.find("fasta");
531 if (itTypes != outputTypes.end()) {
532 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
535 itTypes = outputTypes.find("name");
536 if (itTypes != outputTypes.end()) {
537 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
540 itTypes = outputTypes.find("group");
541 if (itTypes != outputTypes.end()) {
542 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
545 itTypes = outputTypes.find("qfile");
546 if (itTypes != outputTypes.end()) {
547 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setQualFile(current); }
550 itTypes = outputTypes.find("taxonomy");
551 if (itTypes != outputTypes.end()) {
552 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
555 itTypes = outputTypes.find("count");
556 if (itTypes != outputTypes.end()) {
557 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
560 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
561 m->mothurOutEndLine();
565 catch(exception& e) {
566 m->errorOut(e, "ScreenSeqsCommand", "execute");
571 //***************************************************************************************************************
573 int ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
576 m->openInputFile(namefile, inputNames);
577 set<string> badSeqGroups;
578 string seqName, seqList, group;
579 set<string>::iterator it;
580 map<string, string> variables;
581 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(namefile));
582 variables["[extension]"] = m->getExtension(namefile);
583 string goodNameFile = getOutputFileName("name", variables);
584 outputNames.push_back(goodNameFile); outputTypes["name"].push_back(goodNameFile);
586 ofstream goodNameOut; m->openOutputFile(goodNameFile, goodNameOut);
588 while(!inputNames.eof()){
589 if (m->control_pressed) { goodNameOut.close(); inputNames.close(); m->mothurRemove(goodNameFile); return 0; }
591 inputNames >> seqName >> seqList;
592 it = badSeqNames.find(seqName);
594 if(it != badSeqNames.end()){
595 badSeqNames.erase(it);
599 for(int i=0;i<seqList.length();i++){
600 if(seqList[i] == ','){
601 badSeqGroups.insert(seqList.substr(start,i-start));
605 badSeqGroups.insert(seqList.substr(start,seqList.length()-start));
609 goodNameOut << seqName << '\t' << seqList << endl;
611 m->gobble(inputNames);
616 //we were unable to remove some of the bad sequences
617 if (badSeqNames.size() != 0) {
618 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
619 m->mothurOut("Your namefile does not include the sequence " + *it + " please correct.");
620 m->mothurOutEndLine();
626 ifstream inputGroups;
627 m->openInputFile(groupfile, inputGroups);
628 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(groupfile));
629 variables["[extension]"] = m->getExtension(groupfile);
630 string goodGroupFile = getOutputFileName("group", variables);
632 outputNames.push_back(goodGroupFile); outputTypes["group"].push_back(goodGroupFile);
634 ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
636 while(!inputGroups.eof()){
637 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodNameFile); m->mothurRemove(goodGroupFile); return 0; }
639 inputGroups >> seqName >> group;
641 it = badSeqGroups.find(seqName);
643 if(it != badSeqGroups.end()){
644 badSeqGroups.erase(it);
647 goodGroupOut << seqName << '\t' << group << endl;
649 m->gobble(inputGroups);
652 goodGroupOut.close();
654 //we were unable to remove some of the bad sequences
655 if (badSeqGroups.size() != 0) {
656 for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {
657 m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct.");
658 m->mothurOutEndLine();
667 catch(exception& e) {
668 m->errorOut(e, "ScreenSeqsCommand", "screenNameGroupFile");
672 //***************************************************************************************************************
673 int ScreenSeqsCommand::getSummary(vector<unsigned long long>& positions){
676 vector<int> startPosition;
677 vector<int> endPosition;
678 vector<int> seqLength;
679 vector<int> ambigBases;
680 vector<int> longHomoPolymer;
682 vector<unsigned long long> positions;
683 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
684 positions = m->divideFile(fastafile, processors);
685 for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }
687 if(processors == 1){ lines.push_back(linePair(0, 1000)); }
689 int numFastaSeqs = 0;
690 positions = m->setFilePosFasta(fastafile, numFastaSeqs);
691 if (positions.size() < processors) { processors = positions.size(); }
693 //figure out how many sequences you have to process
694 int numSeqsPerProcessor = numFastaSeqs / processors;
695 for (int i = 0; i < processors; i++) {
696 int startIndex = i * numSeqsPerProcessor;
697 if(i == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor; }
698 lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
705 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
708 driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
711 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
713 numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
715 numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile);
718 if (m->control_pressed) { return 0; }
720 // numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
721 // if (m->control_pressed) { return 0; }
724 sort(startPosition.begin(), startPosition.end());
725 sort(endPosition.begin(), endPosition.end());
726 sort(seqLength.begin(), seqLength.end());
727 sort(ambigBases.begin(), ambigBases.end());
728 sort(longHomoPolymer.begin(), longHomoPolymer.end());
730 //numSeqs is the number of unique seqs, startPosition.size() is the total number of seqs, we want to optimize using all seqs
731 int criteriaPercentile = int(startPosition.size() * (criteria / (float) 100));
733 for (int i = 0; i < optimize.size(); i++) {
734 if (optimize[i] == "start") { startPos = startPosition[criteriaPercentile]; m->mothurOut("Optimizing start to " + toString(startPos) + "."); m->mothurOutEndLine(); }
735 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();}
736 else if (optimize[i] == "maxambig") { maxAmbig = ambigBases[criteriaPercentile]; m->mothurOut("Optimizing maxambig to " + toString(maxAmbig) + "."); m->mothurOutEndLine(); }
737 else if (optimize[i] == "maxhomop") { maxHomoP = longHomoPolymer[criteriaPercentile]; m->mothurOut("Optimizing maxhomop to " + toString(maxHomoP) + "."); m->mothurOutEndLine(); }
738 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(); }
739 else if (optimize[i] == "maxlength") { maxLength = seqLength[criteriaPercentile]; m->mothurOut("Optimizing maxlength to " + toString(maxLength) + "."); m->mothurOutEndLine(); }
746 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
747 MPI_Comm_size(MPI_COMM_WORLD, &processors);
750 //send file positions to all processes
751 for(int i = 1; i < processors; i++) {
752 MPI_Send(&startPos, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
753 MPI_Send(&endPos, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
754 MPI_Send(&maxAmbig, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
755 MPI_Send(&maxHomoP, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
756 MPI_Send(&minLength, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
757 MPI_Send(&maxLength, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
760 MPI_Recv(&startPos, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
761 MPI_Recv(&endPos, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
762 MPI_Recv(&maxAmbig, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
763 MPI_Recv(&maxHomoP, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
764 MPI_Recv(&minLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
765 MPI_Recv(&maxLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
767 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
771 catch(exception& e) {
772 m->errorOut(e, "ScreenSeqsCommand", "getSummary");
776 /**************************************************************************************/
777 int ScreenSeqsCommand::driverCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, linePair filePos) {
781 m->openInputFile(filename, in);
783 in.seekg(filePos.start);
790 if (m->control_pressed) { in.close(); return 1; }
792 Sequence current(in); m->gobble(in);
794 if (current.getName() != "") {
796 if (namefile != "") {
797 //make sure this sequence is in the namefile, else error
798 map<string, int>::iterator it = nameMap.find(current.getName());
800 if (it == nameMap.end()) { m->mothurOut("[ERROR]: " + current.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); m->control_pressed = true; }
801 else { num = it->second; }
804 //for each sequence this sequence represents
805 for (int i = 0; i < num; i++) {
806 startPosition.push_back(current.getStartPos());
807 endPosition.push_back(current.getEndPos());
808 seqLength.push_back(current.getNumBases());
809 ambigBases.push_back(current.getAmbigBases());
810 longHomoPolymer.push_back(current.getLongHomoPolymer());
815 //if((count) % 100 == 0){ m->mothurOut("Optimizing sequence: " + toString(count)); m->mothurOutEndLine(); }
816 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
817 unsigned long long pos = in.tellg();
818 if ((pos == -1) || (pos >= filePos.end)) { break; }
820 if (in.eof()) { break; }
829 catch(exception& e) {
830 m->errorOut(e, "ScreenSeqsCommand", "driverCreateSummary");
834 /**************************************************************************************************/
835 int ScreenSeqsCommand::createProcessesCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename) {
840 vector<int> processIDS;
842 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
844 //loop through and create all the processes you want
845 while (process != processors) {
849 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
852 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[process]);
854 //pass numSeqs to parent
856 string tempFile = fastafile + toString(getpid()) + ".num.temp";
857 m->openOutputFile(tempFile, out);
860 out << startPosition.size() << endl;
861 for (int k = 0; k < startPosition.size(); k++) { out << startPosition[k] << '\t'; } out << endl;
862 for (int k = 0; k < endPosition.size(); k++) { out << endPosition[k] << '\t'; } out << endl;
863 for (int k = 0; k < seqLength.size(); k++) { out << seqLength[k] << '\t'; } out << endl;
864 for (int k = 0; k < ambigBases.size(); k++) { out << ambigBases[k] << '\t'; } out << endl;
865 for (int k = 0; k < longHomoPolymer.size(); k++) { out << longHomoPolymer[k] << '\t'; } out << endl;
871 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
872 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
877 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
879 //force parent to wait until all the processes are done
880 for (int i=0;i<processIDS.size();i++) {
881 int temp = processIDS[i];
885 //parent reads in and combine Filter info
886 for (int i = 0; i < processIDS.size(); i++) {
887 string tempFilename = fastafile + toString(processIDS[i]) + ".num.temp";
889 m->openInputFile(tempFilename, in);
892 in >> tempNum; m->gobble(in); num += tempNum;
893 in >> tempNum; m->gobble(in);
894 for (int k = 0; k < tempNum; k++) { in >> temp; startPosition.push_back(temp); } m->gobble(in);
895 for (int k = 0; k < tempNum; k++) { in >> temp; endPosition.push_back(temp); } m->gobble(in);
896 for (int k = 0; k < tempNum; k++) { in >> temp; seqLength.push_back(temp); } m->gobble(in);
897 for (int k = 0; k < tempNum; k++) { in >> temp; ambigBases.push_back(temp); } m->gobble(in);
898 for (int k = 0; k < tempNum; k++) { in >> temp; longHomoPolymer.push_back(temp); } m->gobble(in);
901 m->mothurRemove(tempFilename);
906 //////////////////////////////////////////////////////////////////////////////////////////////////////
907 //Windows version shared memory, so be careful when passing variables through the seqSumData struct.
908 //Above fork() will clone, so memory is separate, but that's not the case with windows,
909 //Taking advantage of shared memory to allow both threads to add info to vectors.
910 //////////////////////////////////////////////////////////////////////////////////////////////////////
912 vector<sumData*> pDataArray;
913 DWORD dwThreadIdArray[processors-1];
914 HANDLE hThreadArray[processors-1];
916 //Create processor worker threads.
917 for( int i=0; i<processors-1; i++ ){
919 // Allocate memory for thread data.
920 sumData* tempSum = new sumData(filename, m, lines[i].start, lines[i].end, namefile, nameMap);
921 pDataArray.push_back(tempSum);
923 //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
924 //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
925 hThreadArray[i] = CreateThread(NULL, 0, MySumThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
929 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[processors-1]);
931 //Wait until all threads have terminated.
932 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
934 //Close all thread handles and free memory allocations.
935 for(int i=0; i < pDataArray.size(); i++){
936 num += pDataArray[i]->count;
937 for (int k = 0; k < pDataArray[i]->startPosition.size(); k++) { startPosition.push_back(pDataArray[i]->startPosition[k]); }
938 for (int k = 0; k < pDataArray[i]->endPosition.size(); k++) { endPosition.push_back(pDataArray[i]->endPosition[k]); }
939 for (int k = 0; k < pDataArray[i]->seqLength.size(); k++) { seqLength.push_back(pDataArray[i]->seqLength[k]); }
940 for (int k = 0; k < pDataArray[i]->ambigBases.size(); k++) { ambigBases.push_back(pDataArray[i]->ambigBases[k]); }
941 for (int k = 0; k < pDataArray[i]->longHomoPolymer.size(); k++) { longHomoPolymer.push_back(pDataArray[i]->longHomoPolymer[k]); }
942 CloseHandle(hThreadArray[i]);
943 delete pDataArray[i];
949 catch(exception& e) {
950 m->errorOut(e, "ScreenSeqsCommand", "createProcessesCreateSummary");
955 //***************************************************************************************************************
957 int ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
959 ifstream inputGroups;
960 m->openInputFile(groupfile, inputGroups);
961 string seqName, group;
962 set<string>::iterator it;
963 map<string, string> variables;
964 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(groupfile));
965 variables["[extension]"] = m->getExtension(groupfile);
966 string goodGroupFile = getOutputFileName("group", variables);
967 outputNames.push_back(goodGroupFile); outputTypes["group"].push_back(goodGroupFile);
968 ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
970 while(!inputGroups.eof()){
971 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodGroupFile); return 0; }
973 inputGroups >> seqName >> group;
974 it = badSeqNames.find(seqName);
976 if(it != badSeqNames.end()){
977 badSeqNames.erase(it);
980 goodGroupOut << seqName << '\t' << group << endl;
982 m->gobble(inputGroups);
985 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodGroupFile); return 0; }
987 //we were unable to remove some of the bad sequences
988 if (badSeqNames.size() != 0) {
989 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
990 m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct.");
991 m->mothurOutEndLine();
996 goodGroupOut.close();
998 if (m->control_pressed) { m->mothurRemove(goodGroupFile); }
1003 catch(exception& e) {
1004 m->errorOut(e, "ScreenSeqsCommand", "screenGroupFile");
1008 //***************************************************************************************************************
1009 int ScreenSeqsCommand::screenCountFile(set<string> badSeqNames){
1012 m->openInputFile(countfile, in);
1013 set<string>::iterator it;
1014 map<string, string> variables;
1015 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(countfile));
1016 variables["[extension]"] = m->getExtension(countfile);
1017 string goodCountFile = getOutputFileName("count", variables);
1019 outputNames.push_back(goodCountFile); outputTypes["count"].push_back(goodCountFile);
1020 ofstream goodCountOut; m->openOutputFile(goodCountFile, goodCountOut);
1022 string headers = m->getline(in); m->gobble(in);
1023 goodCountOut << headers << endl;
1025 string name, rest; int thisTotal;
1028 if (m->control_pressed) { goodCountOut.close(); in.close(); m->mothurRemove(goodCountFile); return 0; }
1030 in >> name; m->gobble(in);
1031 in >> thisTotal; m->gobble(in);
1032 rest = m->getline(in); m->gobble(in);
1034 it = badSeqNames.find(name);
1036 if(it != badSeqNames.end()){
1037 badSeqNames.erase(it);
1040 goodCountOut << name << '\t' << thisTotal << '\t' << rest << endl;
1044 if (m->control_pressed) { goodCountOut.close(); in.close(); m->mothurRemove(goodCountFile); return 0; }
1046 //we were unable to remove some of the bad sequences
1047 if (badSeqNames.size() != 0) {
1048 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
1049 m->mothurOut("Your count file does not include the sequence " + *it + " please correct.");
1050 m->mothurOutEndLine();
1055 goodCountOut.close();
1057 //check for groups that have been eliminated
1059 if (ct.testGroups(goodCountFile)) {
1060 ct.readTable(goodCountFile);
1061 ct.printTable(goodCountFile);
1064 if (m->control_pressed) { m->mothurRemove(goodCountFile); }
1069 catch(exception& e) {
1070 m->errorOut(e, "ScreenSeqsCommand", "screenCountFile");
1074 //***************************************************************************************************************
1076 int ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
1078 ifstream inputAlignReport;
1079 m->openInputFile(alignreport, inputAlignReport);
1080 string seqName, group;
1081 set<string>::iterator it;
1083 map<string, string> variables;
1084 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(alignreport));
1085 string goodAlignReportFile = getOutputFileName("alignreport", variables);
1087 outputNames.push_back(goodAlignReportFile); outputTypes["alignreport"].push_back(goodAlignReportFile);
1088 ofstream goodAlignReportOut; m->openOutputFile(goodAlignReportFile, goodAlignReportOut);
1090 while (!inputAlignReport.eof()) { // need to copy header
1091 char c = inputAlignReport.get();
1092 goodAlignReportOut << c;
1093 if (c == 10 || c == 13){ break; }
1096 while(!inputAlignReport.eof()){
1097 if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); m->mothurRemove(goodAlignReportFile); return 0; }
1099 inputAlignReport >> seqName;
1100 it = badSeqNames.find(seqName);
1102 while (!inputAlignReport.eof()) { // need to copy header
1103 char c = inputAlignReport.get();
1105 if (c == 10 || c == 13){ break; }
1108 if(it != badSeqNames.end()){
1109 badSeqNames.erase(it);
1112 goodAlignReportOut << seqName << '\t' << line;
1114 m->gobble(inputAlignReport);
1117 if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); m->mothurRemove(goodAlignReportFile); return 0; }
1119 //we were unable to remove some of the bad sequences
1120 if (badSeqNames.size() != 0) {
1121 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
1122 m->mothurOut("Your alignreport file does not include the sequence " + *it + " please correct.");
1123 m->mothurOutEndLine();
1127 inputAlignReport.close();
1128 goodAlignReportOut.close();
1130 if (m->control_pressed) { m->mothurRemove(goodAlignReportFile); return 0; }
1135 catch(exception& e) {
1136 m->errorOut(e, "ScreenSeqsCommand", "screenAlignReport");
1141 //***************************************************************************************************************
1143 int ScreenSeqsCommand::screenTaxonomy(set<string> badSeqNames){
1146 m->openInputFile(taxonomy, input);
1147 string seqName, tax;
1148 set<string>::iterator it;
1149 map<string, string> variables;
1150 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(taxonomy));
1151 variables["[extension]"] = m->getExtension(taxonomy);
1152 string goodTaxFile = getOutputFileName("taxonomy", variables);
1154 outputNames.push_back(goodTaxFile); outputTypes["taxonomy"].push_back(goodTaxFile);
1155 ofstream goodTaxOut; m->openOutputFile(goodTaxFile, goodTaxOut);
1157 while(!input.eof()){
1158 if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); return 0; }
1160 input >> seqName >> tax;
1161 it = badSeqNames.find(seqName);
1163 if(it != badSeqNames.end()){ badSeqNames.erase(it); }
1165 goodTaxOut << seqName << '\t' << tax << endl;
1170 if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); return 0; }
1172 //we were unable to remove some of the bad sequences
1173 if (badSeqNames.size() != 0) {
1174 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
1175 m->mothurOut("Your taxonomy file does not include the sequence " + *it + " please correct.");
1176 m->mothurOutEndLine();
1183 if (m->control_pressed) { m->mothurRemove(goodTaxFile); return 0; }
1188 catch(exception& e) {
1189 m->errorOut(e, "ScreenSeqsCommand", "screenTaxonomy");
1194 //***************************************************************************************************************
1196 int ScreenSeqsCommand::screenQual(set<string> badSeqNames){
1199 m->openInputFile(qualfile, in);
1200 set<string>::iterator it;
1201 map<string, string> variables;
1202 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(qualfile));
1203 variables["[extension]"] = m->getExtension(qualfile);
1204 string goodQualFile = getOutputFileName("qfile", variables);
1206 outputNames.push_back(goodQualFile); outputTypes["qfile"].push_back(goodQualFile);
1207 ofstream goodQual; m->openOutputFile(goodQualFile, goodQual);
1211 if (m->control_pressed) { goodQual.close(); in.close(); m->mothurRemove(goodQualFile); return 0; }
1213 string saveName = "";
1219 if (name.length() != 0) {
1220 saveName = name.substr(1);
1223 if (c == 10 || c == 13){ break; }
1230 char letter= in.get();
1231 if(letter == '>'){ in.putback(letter); break; }
1232 else{ scores += letter; }
1237 it = badSeqNames.find(saveName);
1239 if(it != badSeqNames.end()){
1240 badSeqNames.erase(it);
1242 goodQual << name << endl << scores;
1251 //we were unable to remove some of the bad sequences
1252 if (badSeqNames.size() != 0) {
1253 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
1254 m->mothurOut("Your qual file does not include the sequence " + *it + " please correct.");
1255 m->mothurOutEndLine();
1259 if (m->control_pressed) { m->mothurRemove(goodQualFile); return 0; }
1264 catch(exception& e) {
1265 m->errorOut(e, "ScreenSeqsCommand", "screenQual");
1270 //**********************************************************************************************************************
1272 int ScreenSeqsCommand::driver(linePair filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
1275 m->openOutputFile(goodFName, goodFile);
1277 ofstream badAccnosFile;
1278 m->openOutputFile(badAccnosFName, badAccnosFile);
1281 m->openInputFile(filename, inFASTA);
1283 inFASTA.seekg(filePos.start);
1290 if (m->control_pressed) { return 0; }
1292 Sequence currSeq(inFASTA); m->gobble(inFASTA);
1293 if (currSeq.getName() != "") {
1294 bool goodSeq = 1; // innocent until proven guilty
1295 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos()) { goodSeq = 0; }
1296 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos()) { goodSeq = 0; }
1297 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases()) { goodSeq = 0; }
1298 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer()) { goodSeq = 0; }
1299 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases()) { goodSeq = 0; }
1300 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases()) { goodSeq = 0; }
1303 currSeq.printSequence(goodFile);
1306 badAccnosFile << currSeq.getName() << endl;
1307 badSeqNames.insert(currSeq.getName());
1312 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1313 unsigned long long pos = inFASTA.tellg();
1314 if ((pos == -1) || (pos >= filePos.end)) { break; }
1316 if (inFASTA.eof()) { break; }
1320 if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
1323 if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
1328 badAccnosFile.close();
1332 catch(exception& e) {
1333 m->errorOut(e, "ScreenSeqsCommand", "driver");
1337 //**********************************************************************************************************************
1339 int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badAccnosFile, vector<unsigned long long>& MPIPos, set<string>& badSeqNames){
1341 string outputString = "";
1342 MPI_Status statusGood;
1343 MPI_Status statusBadAccnos;
1346 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
1348 for(int i=0;i<num;i++){
1350 if (m->control_pressed) { return 0; }
1352 //read next sequence
1353 int length = MPIPos[start+i+1] - MPIPos[start+i];
1355 char* buf4 = new char[length];
1357 MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
1359 string tempBuf = buf4; delete buf4;
1360 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
1361 istringstream iss (tempBuf,istringstream::in);
1363 Sequence currSeq(iss);
1366 if (currSeq.getName() != "") {
1367 bool goodSeq = 1; // innocent until proven guilty
1368 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos()) { goodSeq = 0; }
1369 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos()) { goodSeq = 0; }
1370 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases()) { goodSeq = 0; }
1371 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer()) { goodSeq = 0; }
1372 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases()) { goodSeq = 0; }
1373 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases()) { goodSeq = 0; }
1376 outputString = ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
1379 length = outputString.length();
1380 char* buf2 = new char[length];
1381 memcpy(buf2, outputString.c_str(), length);
1383 MPI_File_write_shared(goodFile, buf2, length, MPI_CHAR, &statusGood);
1388 badSeqNames.insert(currSeq.getName());
1390 //write to bad accnos file
1391 outputString = currSeq.getName() + "\n";
1393 length = outputString.length();
1394 char* buf3 = new char[length];
1395 memcpy(buf3, outputString.c_str(), length);
1397 MPI_File_write_shared(badAccnosFile, buf3, length, MPI_CHAR, &statusBadAccnos);
1403 if((i) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(i)); m->mothurOutEndLine(); }
1408 catch(exception& e) {
1409 m->errorOut(e, "ScreenSeqsCommand", "driverMPI");
1414 /**************************************************************************************************/
1416 int ScreenSeqsCommand::createProcesses(string goodFileName, string badAccnos, string filename, set<string>& badSeqNames) {
1419 vector<int> processIDS;
1423 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1425 //loop through and create all the processes you want
1426 while (process != processors) {
1430 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1432 }else if (pid == 0){
1433 num = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
1435 //pass numSeqs to parent
1437 string tempFile = filename + toString(getpid()) + ".num.temp";
1438 m->openOutputFile(tempFile, out);
1444 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1445 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1450 num = driver(lines[0], goodFileName, badAccnos, filename, badSeqNames);
1452 //force parent to wait until all the processes are done
1453 for (int i=0;i<processIDS.size();i++) {
1454 int temp = processIDS[i];
1458 for (int i = 0; i < processIDS.size(); i++) {
1460 string tempFile = filename + toString(processIDS[i]) + ".num.temp";
1461 m->openInputFile(tempFile, in);
1462 if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
1463 in.close(); m->mothurRemove(tempFile);
1465 m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
1466 m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
1468 m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
1469 m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
1472 //read badSeqs in because root process doesnt know what other "bad" seqs the children found
1474 int ableToOpen = m->openInputFile(badAccnos, inBad, "no error");
1476 if (ableToOpen == 0) {
1477 badSeqNames.clear();
1479 while (!inBad.eof()) {
1480 inBad >> tempName; m->gobble(inBad);
1481 badSeqNames.insert(tempName);
1487 //////////////////////////////////////////////////////////////////////////////////////////////////////
1488 //Windows version shared memory, so be careful when passing variables through the sumScreenData struct.
1489 //Above fork() will clone, so memory is separate, but that's not the case with windows,
1490 //Taking advantage of shared memory to allow both threads to add info to badSeqNames.
1491 //////////////////////////////////////////////////////////////////////////////////////////////////////
1493 vector<sumScreenData*> pDataArray;
1494 DWORD dwThreadIdArray[processors-1];
1495 HANDLE hThreadArray[processors-1];
1497 //Create processor worker threads.
1498 for( int i=0; i<processors-1; i++ ){
1500 string extension = "";
1501 if (i!=0) {extension += toString(i) + ".temp"; processIDS.push_back(i); }
1503 // Allocate memory for thread data.
1504 sumScreenData* tempSum = new sumScreenData(startPos, endPos, maxAmbig, maxHomoP, minLength, maxLength, filename, m, lines[i].start, lines[i].end,goodFileName+extension, badAccnos+extension);
1505 pDataArray.push_back(tempSum);
1507 //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
1508 hThreadArray[i] = CreateThread(NULL, 0, MySumScreenThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
1512 num = driver(lines[processors-1], (goodFileName+toString(processors-1)+".temp"), (badAccnos+toString(processors-1)+".temp"), filename, badSeqNames);
1513 processIDS.push_back(processors-1);
1515 //Wait until all threads have terminated.
1516 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
1518 //Close all thread handles and free memory allocations.
1519 for(int i=0; i < pDataArray.size(); i++){
1520 num += pDataArray[i]->count;
1521 for (set<string>::iterator it = pDataArray[i]->badSeqNames.begin(); it != pDataArray[i]->badSeqNames.end(); it++) { badSeqNames.insert(*it); }
1522 CloseHandle(hThreadArray[i]);
1523 delete pDataArray[i];
1526 for (int i = 0; i < processIDS.size(); i++) {
1527 m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
1528 m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
1530 m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
1531 m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
1539 catch(exception& e) {
1540 m->errorOut(e, "ScreenSeqsCommand", "createProcesses");
1545 //***************************************************************************************************************