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 "sequence.hpp"
13 //***************************************************************************************************************
15 ScreenSeqsCommand::ScreenSeqsCommand(string option) {
19 //allow user to run help
20 if(option == "help") { help(); abort = true; }
23 //valid paramters for this command
24 string AlignArray[] = {"fasta", "start", "end", "maxambig", "maxhomop", "minlength", "maxlength",
25 "name", "group", "alignreport","processors","outputdir","inputdir"};
26 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
28 OptionParser parser(option);
29 map<string,string> parameters = parser.getParameters();
31 ValidParameters validParameter;
32 map<string,string>::iterator it;
34 //check to make sure all parameters are valid for command
35 for (it = parameters.begin(); it != parameters.end(); it++) {
36 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
39 //if the user changes the input directory command factory will send this info to us in the output parameter
40 string inputDir = validParameter.validFile(parameters, "inputdir", false);
41 if (inputDir == "not found"){ inputDir = ""; }
44 it = parameters.find("fasta");
45 //user has given a template file
46 if(it != parameters.end()){
47 path = hasPath(it->second);
48 //if the user has not given a path then, add inputdir. else leave path alone.
49 if (path == "") { parameters["fasta"] = inputDir + it->second; }
52 it = parameters.find("group");
53 //user has given a template file
54 if(it != parameters.end()){
55 path = hasPath(it->second);
56 //if the user has not given a path then, add inputdir. else leave path alone.
57 if (path == "") { parameters["group"] = inputDir + it->second; }
60 it = parameters.find("name");
61 //user has given a template file
62 if(it != parameters.end()){
63 path = hasPath(it->second);
64 //if the user has not given a path then, add inputdir. else leave path alone.
65 if (path == "") { parameters["name"] = inputDir + it->second; }
68 it = parameters.find("alignreport");
69 //user has given a template file
70 if(it != parameters.end()){
71 path = hasPath(it->second);
72 //if the user has not given a path then, add inputdir. else leave path alone.
73 if (path == "") { parameters["alignreport"] = inputDir + it->second; }
77 //check for required parameters
78 fastafile = validParameter.validFile(parameters, "fasta", true);
79 if (fastafile == "not found") { m->mothurOut("fasta is a required parameter for the screen.seqs command."); m->mothurOutEndLine(); abort = true; }
80 else if (fastafile == "not open") { abort = true; }
82 groupfile = validParameter.validFile(parameters, "group", true);
83 if (groupfile == "not open") { abort = true; }
84 else if (groupfile == "not found") { groupfile = ""; }
86 namefile = validParameter.validFile(parameters, "name", true);
87 if (namefile == "not open") { abort = true; }
88 else if (namefile == "not found") { namefile = ""; }
90 alignreport = validParameter.validFile(parameters, "alignreport", true);
91 if (alignreport == "not open") { abort = true; }
92 else if (alignreport == "not found") { alignreport = ""; }
94 //if the user changes the output directory command factory will send this info to us in the output parameter
95 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
97 outputDir += hasPath(fastafile); //if user entered a file with a path then preserve it
100 //check for optional parameter and set defaults
101 // ...at some point should added some additional type checking...
103 temp = validParameter.validFile(parameters, "start", false); if (temp == "not found") { temp = "-1"; }
104 convert(temp, startPos);
106 temp = validParameter.validFile(parameters, "end", false); if (temp == "not found") { temp = "-1"; }
107 convert(temp, endPos);
109 temp = validParameter.validFile(parameters, "maxambig", false); if (temp == "not found") { temp = "-1"; }
110 convert(temp, maxAmbig);
112 temp = validParameter.validFile(parameters, "maxhomop", false); if (temp == "not found") { temp = "-1"; }
113 convert(temp, maxHomoP);
115 temp = validParameter.validFile(parameters, "minlength", false); if (temp == "not found") { temp = "-1"; }
116 convert(temp, minLength);
118 temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found") { temp = "-1"; }
119 convert(temp, maxLength);
121 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
122 convert(temp, processors);
127 catch(exception& e) {
128 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
132 //**********************************************************************************************************************
134 void ScreenSeqsCommand::help(){
136 m->mothurOut("The screen.seqs command reads a fastafile and creates .....\n");
137 m->mothurOut("The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group and processors.\n");
138 m->mothurOut("The fasta parameter is required.\n");
139 m->mothurOut("The start parameter .... The default is -1.\n");
140 m->mothurOut("The end parameter .... The default is -1.\n");
141 m->mothurOut("The maxambig parameter .... The default is -1.\n");
142 m->mothurOut("The maxhomop parameter .... The default is -1.\n");
143 m->mothurOut("The minlength parameter .... The default is -1.\n");
144 m->mothurOut("The maxlength parameter .... The default is -1.\n");
145 m->mothurOut("The processors parameter allows you to specify the number of processors to use while running the command. The default is 1.\n");
146 m->mothurOut("The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n");
147 m->mothurOut("The screen.seqs command should be in the following format: \n");
148 m->mothurOut("screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig, \n");
149 m->mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n");
150 m->mothurOut("Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n");
151 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
154 catch(exception& e) {
155 m->errorOut(e, "ScreenSeqsCommand", "help");
160 //***************************************************************************************************************
162 ScreenSeqsCommand::~ScreenSeqsCommand(){ /* do nothing */ }
164 //***************************************************************************************************************
166 int ScreenSeqsCommand::execute(){
169 if (abort == true) { return 0; }
171 string goodSeqFile = outputDir + getRootName(getSimpleName(fastafile)) + "good" + getExtension(fastafile);
172 string badSeqFile = outputDir + getRootName(getSimpleName(fastafile)) + "bad" + getExtension(fastafile);
173 string badAccnosFile = outputDir + getRootName(getSimpleName(fastafile)) + "bad.accnos";
175 int numFastaSeqs = 0;
176 set<string> badSeqNames;
177 int start = time(NULL);
180 int pid, end, numSeqsPerProcessor;
185 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
186 MPI_Comm_size(MPI_COMM_WORLD, &processors);
191 MPI_File outMPIBadAccnos;
193 int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
194 int inMode=MPI_MODE_RDONLY;
196 char outGoodFilename[1024];
197 strcpy(outGoodFilename, goodSeqFile.c_str());
199 char outBadFilename[1024];
200 strcpy(outBadFilename, badSeqFile.c_str());
202 char outBadAccnosFilename[1024];
203 strcpy(outBadAccnosFilename, badAccnosFile.c_str());
205 char inFileName[1024];
206 strcpy(inFileName, fastafile.c_str());
208 MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI); //comm, filename, mode, info, filepointer
209 MPI_File_open(MPI_COMM_WORLD, outGoodFilename, outMode, MPI_INFO_NULL, &outMPIGood);
210 MPI_File_open(MPI_COMM_WORLD, outBadFilename, outMode, MPI_INFO_NULL, &outMPIBad);
211 MPI_File_open(MPI_COMM_WORLD, outBadAccnosFilename, outMode, MPI_INFO_NULL, &outMPIBadAccnos);
213 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBad); MPI_File_close(&outMPIBadAccnos); return 0; }
215 if (pid == 0) { //you are the root process
217 MPIPos = setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
219 //send file positions to all processes
220 for(int i = 1; i < processors; i++) {
221 MPI_Send(&numFastaSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
222 MPI_Send(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
225 //figure out how many sequences you have to align
226 numSeqsPerProcessor = numFastaSeqs / processors;
227 int startIndex = pid * numSeqsPerProcessor;
228 if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
231 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBad, outMPIBadAccnos, MPIPos, badSeqNames);
233 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); MPI_File_close(&outMPIBad); return 0; }
235 for (int i = 1; i < processors; i++) {
239 MPI_Recv(&badSize, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
240 /*for (int j = 0; j < badSize; j++) {
242 MPI_Recv(&length, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status); //recv the length of the name
243 char* buf2 = new char[length]; //make space to recieve it
244 MPI_Recv(buf2, length, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); //get name
246 string tempBuf = buf2;
247 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
250 badSeqNames.insert(tempBuf);
253 }else{ //you are a child process
254 MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
255 MPIPos.resize(numFastaSeqs+1);
256 MPI_Recv(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
258 //figure out how many sequences you have to align
259 numSeqsPerProcessor = numFastaSeqs / processors;
260 int startIndex = pid * numSeqsPerProcessor;
261 if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
264 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBad, outMPIBadAccnos, MPIPos, badSeqNames);
266 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBad); MPI_File_close(&outMPIBadAccnos); return 0; }
269 int badSize = badSeqNames.size();
270 MPI_Send(&badSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
273 set<string>::iterator it;
274 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
276 int length = name.length();
277 char* buf2 = new char[length];
278 memcpy(buf2, name.c_str(), length);
280 MPI_Send(&length, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
281 MPI_Send(buf2, length, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
286 MPI_File_close(&inMPI);
287 MPI_File_close(&outMPIGood);
288 MPI_File_close(&outMPIBad);
289 MPI_File_close(&outMPIBadAccnos);
290 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
294 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
297 openInputFile(fastafile, inFASTA);
298 numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
301 lines.push_back(new linePair(0, numFastaSeqs));
303 driver(lines[0], goodSeqFile, badSeqFile, badAccnosFile, fastafile, badSeqNames);
305 if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
308 vector<int> positions;
309 processIDS.resize(0);
312 openInputFile(fastafile, inFASTA);
315 while(!inFASTA.eof()){
316 input = getline(inFASTA);
317 if (input.length() != 0) {
318 if(input[0] == '>'){ long int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1); }
323 numFastaSeqs = positions.size();
325 int numSeqsPerProcessor = numFastaSeqs / processors;
327 for (int i = 0; i < processors; i++) {
328 long int startPos = positions[ i * numSeqsPerProcessor ];
329 if(i == processors - 1){
330 numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
332 lines.push_back(new linePair(startPos, numSeqsPerProcessor));
335 createProcesses(goodSeqFile, badSeqFile, badAccnosFile, fastafile, badSeqNames);
337 rename((goodSeqFile + toString(processIDS[0]) + ".temp").c_str(), goodSeqFile.c_str());
338 rename((badSeqFile + toString(processIDS[0]) + ".temp").c_str(), badSeqFile.c_str());
339 rename((badAccnosFile + toString(processIDS[0]) + ".temp").c_str(), badAccnosFile.c_str());
341 //append alignment and report files
342 for(int i=1;i<processors;i++){
343 appendFiles((goodSeqFile + toString(processIDS[i]) + ".temp"), goodSeqFile);
344 remove((goodSeqFile + toString(processIDS[i]) + ".temp").c_str());
346 appendFiles((badSeqFile + toString(processIDS[i]) + ".temp"), badSeqFile);
347 remove((badSeqFile + toString(processIDS[i]) + ".temp").c_str());
349 appendFiles((badAccnosFile + toString(processIDS[i]) + ".temp"), badAccnosFile);
350 remove((badAccnosFile + toString(processIDS[i]) + ".temp").c_str());
353 if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
355 //read badSeqs in because root process doesnt know what other "bad" seqs the children found
357 int ableToOpen = openInputFile(badAccnosFile, inBad, "no error");
359 if (ableToOpen == 0) {
362 while (!inBad.eof()) {
363 inBad >> tempName; gobble(inBad);
364 badSeqNames.insert(tempName);
371 openInputFile(fastafile, inFASTA);
372 numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
375 lines.push_back(new linePair(0, numFastaSeqs));
377 driver(lines[0], goodSeqFile, badSeqFile, badAccnosFile, fastafile, badSeqNames);
379 if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
386 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
388 if (pid == 0) { //only one process should fix files
390 //read accnos file with all names in it, process 0 just has its names
391 MPI_File inMPIAccnos;
394 char inFileName[1024];
395 strcpy(inFileName, badAccnosFile.c_str());
397 MPI_File_open(MPI_COMM_SELF, inFileName, inMode, MPI_INFO_NULL, &inMPIAccnos); //comm, filename, mode, info, filepointer
398 MPI_File_get_size(inMPIAccnos, &size);
400 char* buffer = new char[size];
401 MPI_File_read(inMPIAccnos, buffer, size, MPI_CHAR, &status);
403 string tempBuf = buffer;
404 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size); }
405 istringstream iss (tempBuf,istringstream::in);
408 MPI_File_close(&inMPIAccnos);
413 iss >> tempName; gobble(iss);
414 badSeqNames.insert(tempName);
418 if(namefile != "" && groupfile != "") {
419 screenNameGroupFile(badSeqNames);
420 if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
421 }else if(namefile != "") {
422 screenNameGroupFile(badSeqNames);
423 if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
424 }else if(groupfile != "") { screenGroupFile(badSeqNames); } // this screens just the group
426 if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
428 if(alignreport != "") { screenAlignReport(badSeqNames); }
430 if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
436 m->mothurOutEndLine();
437 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
438 m->mothurOut(goodSeqFile); m->mothurOutEndLine();
439 m->mothurOut(badSeqFile); m->mothurOutEndLine();
440 m->mothurOut(badAccnosFile); m->mothurOutEndLine();
441 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
442 m->mothurOutEndLine();
443 m->mothurOutEndLine();
445 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
446 m->mothurOutEndLine();
450 catch(exception& e) {
451 m->errorOut(e, "ScreenSeqsCommand", "execute");
456 //***************************************************************************************************************
458 int ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
461 openInputFile(namefile, inputNames);
462 set<string> badSeqGroups;
463 string seqName, seqList, group;
464 set<string>::iterator it;
466 string goodNameFile = outputDir + getRootName(getSimpleName(namefile)) + "good" + getExtension(namefile);
467 string badNameFile = outputDir + getRootName(getSimpleName(namefile)) + "bad" + getExtension(namefile);
469 outputNames.push_back(goodNameFile); outputNames.push_back(badNameFile);
471 ofstream goodNameOut; openOutputFile(goodNameFile, goodNameOut);
472 ofstream badNameOut; openOutputFile(badNameFile, badNameOut);
474 while(!inputNames.eof()){
475 if (m->control_pressed) { goodNameOut.close(); badNameOut.close(); inputNames.close(); remove(goodNameFile.c_str()); remove(badNameFile.c_str()); return 0; }
477 inputNames >> seqName >> seqList;
478 it = badSeqNames.find(seqName);
480 if(it != badSeqNames.end()){
481 badSeqNames.erase(it);
482 badNameOut << seqName << '\t' << seqList << endl;
485 for(int i=0;i<seqList.length();i++){
486 if(seqList[i] == ','){
487 badSeqGroups.insert(seqList.substr(start,i-start));
491 badSeqGroups.insert(seqList.substr(start,seqList.length()-start));
495 goodNameOut << seqName << '\t' << seqList << endl;
503 //we were unable to remove some of the bad sequences
504 if (badSeqNames.size() != 0) {
505 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
506 m->mothurOut("Your namefile does not include the sequence " + *it + " please correct.");
507 m->mothurOutEndLine();
513 ifstream inputGroups;
514 openInputFile(groupfile, inputGroups);
516 string goodGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "good" + getExtension(groupfile);
517 string badGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "bad" + getExtension(groupfile);
519 outputNames.push_back(goodGroupFile); outputNames.push_back(badGroupFile);
521 ofstream goodGroupOut; openOutputFile(goodGroupFile, goodGroupOut);
522 ofstream badGroupOut; openOutputFile(badGroupFile, badGroupOut);
524 while(!inputGroups.eof()){
525 if (m->control_pressed) { goodGroupOut.close(); badGroupOut.close(); inputGroups.close(); remove(goodNameFile.c_str()); remove(badNameFile.c_str()); remove(goodGroupFile.c_str()); remove(badGroupFile.c_str()); return 0; }
527 inputGroups >> seqName >> group;
529 it = badSeqGroups.find(seqName);
531 if(it != badSeqGroups.end()){
532 badSeqGroups.erase(it);
533 badGroupOut << seqName << '\t' << group << endl;
536 goodGroupOut << seqName << '\t' << group << endl;
541 goodGroupOut.close();
544 //we were unable to remove some of the bad sequences
545 if (badSeqGroups.size() != 0) {
546 for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {
547 m->mothurOut("Your namefile does not include the sequence " + *it + " please correct.");
548 m->mothurOutEndLine();
556 catch(exception& e) {
557 m->errorOut(e, "ScreenSeqsCommand", "screenNameGroupFile");
562 //***************************************************************************************************************
564 int ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
566 ifstream inputGroups;
567 openInputFile(groupfile, inputGroups);
568 string seqName, group;
569 set<string>::iterator it;
571 string goodGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "good" + getExtension(groupfile);
572 string badGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "bad" + getExtension(groupfile);
574 outputNames.push_back(goodGroupFile); outputNames.push_back(badGroupFile);
576 ofstream goodGroupOut; openOutputFile(goodGroupFile, goodGroupOut);
577 ofstream badGroupOut; openOutputFile(badGroupFile, badGroupOut);
579 while(!inputGroups.eof()){
580 if (m->control_pressed) { goodGroupOut.close(); badGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); remove(badGroupFile.c_str()); return 0; }
582 inputGroups >> seqName >> group;
583 it = badSeqNames.find(seqName);
585 if(it != badSeqNames.end()){
586 badSeqNames.erase(it);
587 badGroupOut << seqName << '\t' << group << endl;
590 goodGroupOut << seqName << '\t' << group << endl;
595 if (m->control_pressed) { goodGroupOut.close(); badGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); remove(badGroupFile.c_str()); return 0; }
597 //we were unable to remove some of the bad sequences
598 if (badSeqNames.size() != 0) {
599 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
600 m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct.");
601 m->mothurOutEndLine();
606 goodGroupOut.close();
609 if (m->control_pressed) { remove(goodGroupFile.c_str()); remove(badGroupFile.c_str()); }
615 catch(exception& e) {
616 m->errorOut(e, "ScreenSeqsCommand", "screenGroupFile");
621 //***************************************************************************************************************
623 int ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
625 ifstream inputAlignReport;
626 openInputFile(alignreport, inputAlignReport);
627 string seqName, group;
628 set<string>::iterator it;
630 string goodAlignReportFile = outputDir + getRootName(getSimpleName(alignreport)) + "good" + getExtension(alignreport);
631 string badAlignReportFile = outputDir + getRootName(getSimpleName(alignreport)) + "bad" + getExtension(alignreport);
633 outputNames.push_back(goodAlignReportFile); outputNames.push_back(badAlignReportFile);
635 ofstream goodAlignReportOut; openOutputFile(goodAlignReportFile, goodAlignReportOut);
636 ofstream badAlignReportOut; openOutputFile(badAlignReportFile, badAlignReportOut);
638 while (!inputAlignReport.eof()) { // need to copy header
639 char c = inputAlignReport.get();
640 goodAlignReportOut << c;
641 badAlignReportOut << c;
642 if (c == 10 || c == 13){ break; }
645 while(!inputAlignReport.eof()){
646 if (m->control_pressed) { goodAlignReportOut.close(); badAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); remove(badAlignReportFile.c_str()); return 0; }
648 inputAlignReport >> seqName;
649 it = badSeqNames.find(seqName);
651 while (!inputAlignReport.eof()) { // need to copy header
652 char c = inputAlignReport.get();
654 if (c == 10 || c == 13){ break; }
657 if(it != badSeqNames.end()){
658 badSeqNames.erase(it);
659 badAlignReportOut << seqName << '\t' << line;
662 goodAlignReportOut << seqName << '\t' << line;
664 gobble(inputAlignReport);
667 if (m->control_pressed) { goodAlignReportOut.close(); badAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); remove(badAlignReportFile.c_str()); return 0; }
669 //we were unable to remove some of the bad sequences
670 if (badSeqNames.size() != 0) {
671 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
672 m->mothurOut("Your file does not include the sequence " + *it + " please correct.");
673 m->mothurOutEndLine();
677 inputAlignReport.close();
678 goodAlignReportOut.close();
679 badAlignReportOut.close();
681 if (m->control_pressed) { remove(goodAlignReportFile.c_str()); remove(badAlignReportFile.c_str()); return 0; }
686 catch(exception& e) {
687 m->errorOut(e, "ScreenSeqsCommand", "screenAlignReport");
692 //**********************************************************************************************************************
694 int ScreenSeqsCommand::driver(linePair* line, string goodFName, string badFName, string badAccnosFName, string filename, set<string>& badSeqNames){
697 openOutputFile(goodFName, goodFile);
700 openOutputFile(badFName, badFile);
702 ofstream badAccnosFile;
703 openOutputFile(badAccnosFName, badAccnosFile);
706 openInputFile(filename, inFASTA);
708 inFASTA.seekg(line->start);
710 for(int i=0;i<line->numSeqs;i++){
712 if (m->control_pressed) { return 0; }
714 Sequence currSeq(inFASTA);
715 if (currSeq.getName() != "") {
716 bool goodSeq = 1; // innocent until proven guilty
717 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos()) { goodSeq = 0; }
718 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos()) { goodSeq = 0; }
719 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases()) { goodSeq = 0; }
720 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer()) { goodSeq = 0; }
721 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases()) { goodSeq = 0; }
722 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases()) { goodSeq = 0; }
725 currSeq.printSequence(goodFile);
728 currSeq.printSequence(badFile);
729 badAccnosFile << currSeq.getName() << endl;
730 badSeqNames.insert(currSeq.getName());
740 badAccnosFile.close();
744 catch(exception& e) {
745 m->errorOut(e, "ScreenSeqsCommand", "driver");
749 //**********************************************************************************************************************
751 int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badFile, MPI_File& badAccnosFile, vector<long>& MPIPos, set<string>& badSeqNames){
753 string outputString = "";
754 MPI_Status statusGood;
755 MPI_Status statusBad;
756 MPI_Status statusBadAccnos;
759 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
761 for(int i=0;i<num;i++){
763 if (m->control_pressed) { return 0; }
766 int length = MPIPos[start+i+1] - MPIPos[start+i];
768 char* buf4 = new char[length];
769 memcpy(buf4, outputString.c_str(), length);
771 MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
773 string tempBuf = buf4; delete buf4;
774 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
775 istringstream iss (tempBuf,istringstream::in);
777 Sequence currSeq(iss);
780 if (currSeq.getName() != "") {
781 bool goodSeq = 1; // innocent until proven guilty
782 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos()) { goodSeq = 0; }
783 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos()) { goodSeq = 0; }
784 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases()) { goodSeq = 0; }
785 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer()) { goodSeq = 0; }
786 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases()) { goodSeq = 0; }
787 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases()) { goodSeq = 0; }
790 outputString = ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
793 length = outputString.length();
794 char* buf2 = new char[length];
795 memcpy(buf2, outputString.c_str(), length);
797 MPI_File_write_shared(goodFile, buf2, length, MPI_CHAR, &statusGood);
801 outputString = ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
803 //print bad seq to fasta
804 length = outputString.length();
805 char* buf2 = new char[length];
806 memcpy(buf2, outputString.c_str(), length);
808 MPI_File_write_shared(badFile, buf2, length, MPI_CHAR, &statusBad);
811 badSeqNames.insert(currSeq.getName());
813 //write to bad accnos file
814 outputString = currSeq.getName() + "\n";
816 length = outputString.length();
817 char* buf3 = new char[length];
818 memcpy(buf3, outputString.c_str(), length);
820 MPI_File_write_shared(badAccnosFile, buf3, length, MPI_CHAR, &statusBadAccnos);
828 catch(exception& e) {
829 m->errorOut(e, "ScreenSeqsCommand", "driverMPI");
834 /**************************************************************************************************/
836 int ScreenSeqsCommand::createProcesses(string goodFileName, string badFileName, string badAccnos, string filename, set<string>& badSeqNames) {
838 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
842 //loop through and create all the processes you want
843 while (process != processors) {
847 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
850 exitCommand = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
852 }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
855 //force parent to wait until all the processes are done
856 for (int i=0;i<processors;i++) {
857 int temp = processIDS[i];
864 catch(exception& e) {
865 m->errorOut(e, "ScreenSeqsCommand", "createProcesses");
870 //***************************************************************************************************************