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 //**********************************************************************************************************************
14 vector<string> ScreenSeqsCommand::getValidParameters(){
16 string Array[] = {"fasta", "start", "end", "maxambig", "maxhomop","optimize","criteria", "minlength", "maxlength",
17 "name", "group", "alignreport","processors","outputdir","inputdir"};
18 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
22 m->errorOut(e, "ScreenSeqsCommand", "getValidParameters");
26 //**********************************************************************************************************************
27 ScreenSeqsCommand::ScreenSeqsCommand(){
29 abort = true; calledHelp = true;
30 vector<string> tempOutNames;
31 outputTypes["fasta"] = tempOutNames;
32 outputTypes["name"] = tempOutNames;
33 outputTypes["group"] = tempOutNames;
34 outputTypes["alignreport"] = tempOutNames;
35 outputTypes["accnos"] = tempOutNames;
38 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
42 //**********************************************************************************************************************
43 vector<string> ScreenSeqsCommand::getRequiredParameters(){
45 string Array[] = {"fasta"};
46 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
50 m->errorOut(e, "ScreenSeqsCommand", "getRequiredParameters");
54 //**********************************************************************************************************************
55 vector<string> ScreenSeqsCommand::getRequiredFiles(){
57 vector<string> myArray;
61 m->errorOut(e, "ScreenSeqsCommand", "getRequiredFiles");
65 //***************************************************************************************************************
67 ScreenSeqsCommand::ScreenSeqsCommand(string option) {
69 abort = false; calledHelp = false;
71 //allow user to run help
72 if(option == "help") { help(); abort = true; calledHelp = true; }
75 //valid paramters for this command
76 string AlignArray[] = {"fasta", "start", "end", "maxambig", "maxhomop","optimize","criteria", "minlength", "maxlength",
77 "name", "group", "alignreport","processors","outputdir","inputdir"};
78 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
80 OptionParser parser(option);
81 map<string,string> parameters = parser.getParameters();
83 ValidParameters validParameter("screen.seqs");
84 map<string,string>::iterator it;
86 //check to make sure all parameters are valid for command
87 for (it = parameters.begin(); it != parameters.end(); it++) {
88 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
91 //initialize outputTypes
92 vector<string> tempOutNames;
93 outputTypes["fasta"] = tempOutNames;
94 outputTypes["name"] = tempOutNames;
95 outputTypes["group"] = tempOutNames;
96 outputTypes["alignreport"] = tempOutNames;
97 outputTypes["accnos"] = tempOutNames;
99 //if the user changes the input directory command factory will send this info to us in the output parameter
100 string inputDir = validParameter.validFile(parameters, "inputdir", false);
101 if (inputDir == "not found"){ inputDir = ""; }
104 it = parameters.find("fasta");
105 //user has given a template file
106 if(it != parameters.end()){
107 path = m->hasPath(it->second);
108 //if the user has not given a path then, add inputdir. else leave path alone.
109 if (path == "") { parameters["fasta"] = inputDir + it->second; }
112 it = parameters.find("group");
113 //user has given a template file
114 if(it != parameters.end()){
115 path = m->hasPath(it->second);
116 //if the user has not given a path then, add inputdir. else leave path alone.
117 if (path == "") { parameters["group"] = inputDir + it->second; }
120 it = parameters.find("name");
121 //user has given a template file
122 if(it != parameters.end()){
123 path = m->hasPath(it->second);
124 //if the user has not given a path then, add inputdir. else leave path alone.
125 if (path == "") { parameters["name"] = inputDir + it->second; }
128 it = parameters.find("alignreport");
129 //user has given a template file
130 if(it != parameters.end()){
131 path = m->hasPath(it->second);
132 //if the user has not given a path then, add inputdir. else leave path alone.
133 if (path == "") { parameters["alignreport"] = inputDir + it->second; }
137 //check for required parameters
138 fastafile = validParameter.validFile(parameters, "fasta", true);
139 if (fastafile == "not found") { m->mothurOut("fasta is a required parameter for the screen.seqs command."); m->mothurOutEndLine(); abort = true; }
140 else if (fastafile == "not open") { abort = true; }
142 groupfile = validParameter.validFile(parameters, "group", true);
143 if (groupfile == "not open") { abort = true; }
144 else if (groupfile == "not found") { groupfile = ""; }
146 namefile = validParameter.validFile(parameters, "name", true);
147 if (namefile == "not open") { namefile = ""; abort = true; }
148 else if (namefile == "not found") { namefile = ""; }
150 alignreport = validParameter.validFile(parameters, "alignreport", true);
151 if (alignreport == "not open") { abort = true; }
152 else if (alignreport == "not found") { alignreport = ""; }
154 //if the user changes the output directory command factory will send this info to us in the output parameter
155 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
157 outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it
160 //check for optional parameter and set defaults
161 // ...at some point should added some additional type checking...
163 temp = validParameter.validFile(parameters, "start", false); if (temp == "not found") { temp = "-1"; }
164 convert(temp, startPos);
166 temp = validParameter.validFile(parameters, "end", false); if (temp == "not found") { temp = "-1"; }
167 convert(temp, endPos);
169 temp = validParameter.validFile(parameters, "maxambig", false); if (temp == "not found") { temp = "-1"; }
170 convert(temp, maxAmbig);
172 temp = validParameter.validFile(parameters, "maxhomop", false); if (temp == "not found") { temp = "-1"; }
173 convert(temp, maxHomoP);
175 temp = validParameter.validFile(parameters, "minlength", false); if (temp == "not found") { temp = "-1"; }
176 convert(temp, minLength);
178 temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found") { temp = "-1"; }
179 convert(temp, maxLength);
181 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
182 convert(temp, processors);
184 temp = validParameter.validFile(parameters, "optimize", false); //optimizing trumps the optimized values original value
185 if (temp == "not found"){ temp = ""; }
186 else { m->splitAtDash(temp, optimize); }
188 //check for invalid optimize options
189 set<string> validOptimizers;
190 validOptimizers.insert("start"); validOptimizers.insert("end"); validOptimizers.insert("maxambig"); validOptimizers.insert("maxhomop"); validOptimizers.insert("minlength"); validOptimizers.insert("maxlength");
192 for (int i = 0; i < optimize.size(); i++) {
193 if (validOptimizers.count(optimize[i]) == 0) {
194 m->mothurOut(optimize[i] + " is not a valid optimizer. Valid options are start, end, maxambig, maxhomop, minlength and maxlength."); m->mothurOutEndLine();
195 optimize.erase(optimize.begin()+i);
200 temp = validParameter.validFile(parameters, "criteria", false); if (temp == "not found"){ temp = "90"; }
201 convert(temp, criteria);
205 catch(exception& e) {
206 m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
210 //**********************************************************************************************************************
212 void ScreenSeqsCommand::help(){
214 m->mothurOut("The screen.seqs command reads a fastafile and creates .....\n");
215 m->mothurOut("The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group, optimize, criteria and processors.\n");
216 m->mothurOut("The fasta parameter is required.\n");
217 m->mothurOut("The start parameter .... The default is -1.\n");
218 m->mothurOut("The end parameter .... The default is -1.\n");
219 m->mothurOut("The maxambig parameter .... The default is -1.\n");
220 m->mothurOut("The maxhomop parameter .... The default is -1.\n");
221 m->mothurOut("The minlength parameter .... The default is -1.\n");
222 m->mothurOut("The maxlength parameter .... The default is -1.\n");
223 m->mothurOut("The processors parameter allows you to specify the number of processors to use while running the command. The default is 1.\n");
224 m->mothurOut("The optimize and criteria parameters allow you set the start, end, maxabig, maxhomop, minlength and maxlength parameters relative to your set of sequences .\n");
225 m->mothurOut("For example optimize=start-end, criteria=90, would set the start and end values to the position 90% of your sequences started and ended.\n");
226 m->mothurOut("The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n");
227 m->mothurOut("The screen.seqs command should be in the following format: \n");
228 m->mothurOut("screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig, \n");
229 m->mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n");
230 m->mothurOut("Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n");
231 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
234 catch(exception& e) {
235 m->errorOut(e, "ScreenSeqsCommand", "help");
240 //***************************************************************************************************************
242 ScreenSeqsCommand::~ScreenSeqsCommand(){ /* do nothing */ }
244 //***************************************************************************************************************
246 int ScreenSeqsCommand::execute(){
249 if (abort == true) { if (calledHelp) { return 0; } return 2; }
251 //if the user want to optimize we need to know the 90% mark
252 vector<unsigned long int> positions;
253 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
254 //use the namefile to optimize correctly
255 if (namefile != "") { nameMap = m->readNames(namefile); }
256 getSummary(positions);
259 positions = m->divideFile(fastafile, processors);
260 for (int i = 0; i < (positions.size()-1); i++) {
261 lines.push_back(new linePair(positions[i], positions[(i+1)]));
265 string goodSeqFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "good" + m->getExtension(fastafile);
266 string badAccnosFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "bad.accnos";
268 int numFastaSeqs = 0;
269 set<string> badSeqNames;
270 int start = time(NULL);
273 int pid, numSeqsPerProcessor;
275 vector<unsigned long int> MPIPos;
278 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
279 MPI_Comm_size(MPI_COMM_WORLD, &processors);
283 MPI_File outMPIBadAccnos;
285 int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
286 int inMode=MPI_MODE_RDONLY;
288 char outGoodFilename[1024];
289 strcpy(outGoodFilename, goodSeqFile.c_str());
291 char outBadAccnosFilename[1024];
292 strcpy(outBadAccnosFilename, badAccnosFile.c_str());
294 char inFileName[1024];
295 strcpy(inFileName, fastafile.c_str());
297 MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI); //comm, filename, mode, info, filepointer
298 MPI_File_open(MPI_COMM_WORLD, outGoodFilename, outMode, MPI_INFO_NULL, &outMPIGood);
299 MPI_File_open(MPI_COMM_WORLD, outBadAccnosFilename, outMode, MPI_INFO_NULL, &outMPIBadAccnos);
301 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
303 if (pid == 0) { //you are the root process
305 MPIPos = m->setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
307 //send file positions to all processes
308 for(int i = 1; i < processors; i++) {
309 MPI_Send(&numFastaSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
310 MPI_Send(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
313 //figure out how many sequences you have to align
314 numSeqsPerProcessor = numFastaSeqs / processors;
315 int startIndex = pid * numSeqsPerProcessor;
316 if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
317 // cout << pid << '\t' << numSeqsPerProcessor << '\t' << startIndex << endl;
319 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
320 //cout << pid << " done" << endl;
321 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
323 for (int i = 1; i < processors; i++) {
327 MPI_Recv(&badSize, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
328 /*for (int j = 0; j < badSize; j++) {
330 MPI_Recv(&length, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status); //recv the length of the name
331 char* buf2 = new char[length]; //make space to recieve it
332 MPI_Recv(buf2, length, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); //get name
334 string tempBuf = buf2;
335 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
338 badSeqNames.insert(tempBuf);
341 }else{ //you are a child process
342 MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
343 MPIPos.resize(numFastaSeqs+1);
344 MPI_Recv(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
346 //figure out how many sequences you have to align
347 numSeqsPerProcessor = numFastaSeqs / processors;
348 int startIndex = pid * numSeqsPerProcessor;
349 if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
350 //cout << pid << '\t' << numSeqsPerProcessor << '\t' << startIndex << endl;
352 driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
353 //cout << pid << " done" << endl;
354 if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
357 int badSize = badSeqNames.size();
358 MPI_Send(&badSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
361 set<string>::iterator it;
362 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
364 int length = name.length();
365 char* buf2 = new char[length];
366 memcpy(buf2, name.c_str(), length);
368 MPI_Send(&length, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
369 MPI_Send(buf2, length, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
374 MPI_File_close(&inMPI);
375 MPI_File_close(&outMPIGood);
376 MPI_File_close(&outMPIBadAccnos);
377 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
381 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
383 numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
385 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
388 processIDS.resize(0);
390 numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames);
392 rename((goodSeqFile + toString(processIDS[0]) + ".temp").c_str(), goodSeqFile.c_str());
393 rename((badAccnosFile + toString(processIDS[0]) + ".temp").c_str(), badAccnosFile.c_str());
395 //append alignment and report files
396 for(int i=1;i<processors;i++){
397 m->appendFiles((goodSeqFile + toString(processIDS[i]) + ".temp"), goodSeqFile);
398 remove((goodSeqFile + toString(processIDS[i]) + ".temp").c_str());
400 m->appendFiles((badAccnosFile + toString(processIDS[i]) + ".temp"), badAccnosFile);
401 remove((badAccnosFile + toString(processIDS[i]) + ".temp").c_str());
404 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
406 //read badSeqs in because root process doesnt know what other "bad" seqs the children found
408 int ableToOpen = m->openInputFile(badAccnosFile, inBad, "no error");
410 if (ableToOpen == 0) {
413 while (!inBad.eof()) {
414 inBad >> tempName; m->gobble(inBad);
415 badSeqNames.insert(tempName);
421 numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
423 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
430 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
432 if (pid == 0) { //only one process should fix files
434 //read accnos file with all names in it, process 0 just has its names
435 MPI_File inMPIAccnos;
438 char inFileName[1024];
439 strcpy(inFileName, badAccnosFile.c_str());
441 MPI_File_open(MPI_COMM_SELF, inFileName, inMode, MPI_INFO_NULL, &inMPIAccnos); //comm, filename, mode, info, filepointer
442 MPI_File_get_size(inMPIAccnos, &size);
444 char* buffer = new char[size];
445 MPI_File_read(inMPIAccnos, buffer, size, MPI_CHAR, &status);
447 string tempBuf = buffer;
448 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size); }
449 istringstream iss (tempBuf,istringstream::in);
452 MPI_File_close(&inMPIAccnos);
457 iss >> tempName; m->gobble(iss);
458 badSeqNames.insert(tempName);
462 if(namefile != "" && groupfile != "") {
463 screenNameGroupFile(badSeqNames);
464 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
465 }else if(namefile != "") {
466 screenNameGroupFile(badSeqNames);
467 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
468 }else if(groupfile != "") { screenGroupFile(badSeqNames); } // this screens just the group
470 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
472 if(alignreport != "") { screenAlignReport(badSeqNames); }
474 if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
480 m->mothurOutEndLine();
481 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
482 m->mothurOut(goodSeqFile); m->mothurOutEndLine(); outputTypes["fasta"].push_back(goodSeqFile);
483 m->mothurOut(badAccnosFile); m->mothurOutEndLine(); outputTypes["accnos"].push_back(badAccnosFile);
484 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
485 m->mothurOutEndLine();
486 m->mothurOutEndLine();
488 //set fasta file as new current fastafile
490 itTypes = outputTypes.find("fasta");
491 if (itTypes != outputTypes.end()) {
492 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
495 itTypes = outputTypes.find("name");
496 if (itTypes != outputTypes.end()) {
497 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
500 itTypes = outputTypes.find("group");
501 if (itTypes != outputTypes.end()) {
502 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
505 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
506 m->mothurOutEndLine();
510 catch(exception& e) {
511 m->errorOut(e, "ScreenSeqsCommand", "execute");
516 //***************************************************************************************************************
518 int ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
521 m->openInputFile(namefile, inputNames);
522 set<string> badSeqGroups;
523 string seqName, seqList, group;
524 set<string>::iterator it;
526 string goodNameFile = outputDir + m->getRootName(m->getSimpleName(namefile)) + "good" + m->getExtension(namefile);
527 outputNames.push_back(goodNameFile); outputTypes["name"].push_back(goodNameFile);
529 ofstream goodNameOut; m->openOutputFile(goodNameFile, goodNameOut);
531 while(!inputNames.eof()){
532 if (m->control_pressed) { goodNameOut.close(); inputNames.close(); remove(goodNameFile.c_str()); return 0; }
534 inputNames >> seqName >> seqList;
535 it = badSeqNames.find(seqName);
537 if(it != badSeqNames.end()){
538 badSeqNames.erase(it);
542 for(int i=0;i<seqList.length();i++){
543 if(seqList[i] == ','){
544 badSeqGroups.insert(seqList.substr(start,i-start));
548 badSeqGroups.insert(seqList.substr(start,seqList.length()-start));
552 goodNameOut << seqName << '\t' << seqList << endl;
554 m->gobble(inputNames);
559 //we were unable to remove some of the bad sequences
560 if (badSeqNames.size() != 0) {
561 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
562 m->mothurOut("Your namefile does not include the sequence " + *it + " please correct.");
563 m->mothurOutEndLine();
569 ifstream inputGroups;
570 m->openInputFile(groupfile, inputGroups);
572 string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
573 outputNames.push_back(goodGroupFile); outputTypes["group"].push_back(goodGroupFile);
575 ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
577 while(!inputGroups.eof()){
578 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodNameFile.c_str()); remove(goodGroupFile.c_str()); return 0; }
580 inputGroups >> seqName >> group;
582 it = badSeqGroups.find(seqName);
584 if(it != badSeqGroups.end()){
585 badSeqGroups.erase(it);
588 goodGroupOut << seqName << '\t' << group << endl;
590 m->gobble(inputGroups);
593 goodGroupOut.close();
595 //we were unable to remove some of the bad sequences
596 if (badSeqGroups.size() != 0) {
597 for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {
598 m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct.");
599 m->mothurOutEndLine();
608 catch(exception& e) {
609 m->errorOut(e, "ScreenSeqsCommand", "screenNameGroupFile");
613 //***************************************************************************************************************
614 int ScreenSeqsCommand::getSummary(vector<unsigned long int>& positions){
617 vector<int> startPosition;
618 vector<int> endPosition;
619 vector<int> seqLength;
620 vector<int> ambigBases;
621 vector<int> longHomoPolymer;
623 vector<unsigned long int> positions = m->divideFile(fastafile, processors);
625 for (int i = 0; i < (positions.size()-1); i++) {
626 lines.push_back(new linePair(positions[i], positions[(i+1)]));
630 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
632 numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
634 numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile);
637 if (m->control_pressed) { return 0; }
639 numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
640 if (m->control_pressed) { return 0; }
643 sort(startPosition.begin(), startPosition.end());
644 sort(endPosition.begin(), endPosition.end());
645 sort(seqLength.begin(), seqLength.end());
646 sort(ambigBases.begin(), ambigBases.end());
647 sort(longHomoPolymer.begin(), longHomoPolymer.end());
649 //numSeqs is the number of unique seqs, startPosition.size() is the total number of seqs, we want to optimize using all seqs
650 int criteriaPercentile = int(startPosition.size() * (criteria / (float) 100));
652 for (int i = 0; i < optimize.size(); i++) {
653 if (optimize[i] == "start") { startPos = startPosition[criteriaPercentile]; m->mothurOut("Optimizing start to " + toString(startPos) + "."); m->mothurOutEndLine(); }
654 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();}
655 else if (optimize[i] == "maxambig") { maxAmbig = ambigBases[criteriaPercentile]; m->mothurOut("Optimizing maxambig to " + toString(maxAmbig) + "."); m->mothurOutEndLine(); }
656 else if (optimize[i] == "maxhomop") { maxHomoP = longHomoPolymer[criteriaPercentile]; m->mothurOut("Optimizing maxhomop to " + toString(maxHomoP) + "."); m->mothurOutEndLine(); }
657 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(); }
658 else if (optimize[i] == "maxlength") { maxLength = seqLength[criteriaPercentile]; m->mothurOut("Optimizing maxlength to " + toString(maxLength) + "."); m->mothurOutEndLine(); }
663 catch(exception& e) {
664 m->errorOut(e, "ScreenSeqsCommand", "getSummary");
668 /**************************************************************************************/
669 int ScreenSeqsCommand::driverCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, linePair* filePos) {
673 m->openInputFile(filename, in);
675 in.seekg(filePos->start);
682 if (m->control_pressed) { in.close(); return 1; }
684 Sequence current(in); m->gobble(in);
686 if (current.getName() != "") {
688 if (namefile != "") {
689 //make sure this sequence is in the namefile, else error
690 map<string, int>::iterator it = nameMap.find(current.getName());
692 if (it == nameMap.end()) { m->mothurOut("[ERROR]: " + current.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); m->control_pressed = true; }
693 else { num = it->second; }
696 //for each sequence this sequence represents
697 for (int i = 0; i < num; i++) {
698 startPosition.push_back(current.getStartPos());
699 endPosition.push_back(current.getEndPos());
700 seqLength.push_back(current.getNumBases());
701 ambigBases.push_back(current.getAmbigBases());
702 longHomoPolymer.push_back(current.getLongHomoPolymer());
708 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
709 unsigned long int pos = in.tellg();
710 if ((pos == -1) || (pos >= filePos->end)) { break; }
712 if (in.eof()) { break; }
721 catch(exception& e) {
722 m->errorOut(e, "ScreenSeqsCommand", "driverCreateSummary");
726 /**************************************************************************************************/
727 int ScreenSeqsCommand::createProcessesCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename) {
729 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
734 //loop through and create all the processes you want
735 while (process != processors) {
739 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
742 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[process]);
744 //pass numSeqs to parent
746 string tempFile = fastafile + toString(getpid()) + ".num.temp";
747 m->openOutputFile(tempFile, out);
750 out << startPosition.size() << endl;
751 for (int k = 0; k < startPosition.size(); k++) { out << startPosition[k] << '\t'; } out << endl;
752 for (int k = 0; k < endPosition.size(); k++) { out << endPosition[k] << '\t'; } out << endl;
753 for (int k = 0; k < seqLength.size(); k++) { out << seqLength[k] << '\t'; } out << endl;
754 for (int k = 0; k < ambigBases.size(); k++) { out << ambigBases[k] << '\t'; } out << endl;
755 for (int k = 0; k < longHomoPolymer.size(); k++) { out << longHomoPolymer[k] << '\t'; } out << endl;
761 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
762 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
767 num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
769 //force parent to wait until all the processes are done
770 for (int i=0;i<processIDS.size();i++) {
771 int temp = processIDS[i];
775 //parent reads in and combine Filter info
776 for (int i = 0; i < processIDS.size(); i++) {
777 string tempFilename = fastafile + toString(processIDS[i]) + ".num.temp";
779 m->openInputFile(tempFilename, in);
782 in >> tempNum; m->gobble(in); num += tempNum;
783 in >> tempNum; m->gobble(in);
784 for (int k = 0; k < tempNum; k++) { in >> temp; startPosition.push_back(temp); } m->gobble(in);
785 for (int k = 0; k < tempNum; k++) { in >> temp; endPosition.push_back(temp); } m->gobble(in);
786 for (int k = 0; k < tempNum; k++) { in >> temp; seqLength.push_back(temp); } m->gobble(in);
787 for (int k = 0; k < tempNum; k++) { in >> temp; ambigBases.push_back(temp); } m->gobble(in);
788 for (int k = 0; k < tempNum; k++) { in >> temp; longHomoPolymer.push_back(temp); } m->gobble(in);
791 remove(tempFilename.c_str());
797 catch(exception& e) {
798 m->errorOut(e, "ScreenSeqsCommand", "createProcessesCreateSummary");
803 //***************************************************************************************************************
805 int ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
807 ifstream inputGroups;
808 m->openInputFile(groupfile, inputGroups);
809 string seqName, group;
810 set<string>::iterator it;
812 string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
813 outputNames.push_back(goodGroupFile); outputTypes["group"].push_back(goodGroupFile);
814 ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
816 while(!inputGroups.eof()){
817 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); return 0; }
819 inputGroups >> seqName >> group;
820 it = badSeqNames.find(seqName);
822 if(it != badSeqNames.end()){
823 badSeqNames.erase(it);
826 goodGroupOut << seqName << '\t' << group << endl;
828 m->gobble(inputGroups);
831 if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); return 0; }
833 //we were unable to remove some of the bad sequences
834 if (badSeqNames.size() != 0) {
835 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
836 m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct.");
837 m->mothurOutEndLine();
842 goodGroupOut.close();
844 if (m->control_pressed) { remove(goodGroupFile.c_str()); }
849 catch(exception& e) {
850 m->errorOut(e, "ScreenSeqsCommand", "screenGroupFile");
855 //***************************************************************************************************************
857 int ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
859 ifstream inputAlignReport;
860 m->openInputFile(alignreport, inputAlignReport);
861 string seqName, group;
862 set<string>::iterator it;
864 string goodAlignReportFile = outputDir + m->getRootName(m->getSimpleName(alignreport)) + "good" + m->getExtension(alignreport);
865 outputNames.push_back(goodAlignReportFile); outputTypes["alignreport"].push_back(goodAlignReportFile);
866 ofstream goodAlignReportOut; m->openOutputFile(goodAlignReportFile, goodAlignReportOut);
868 while (!inputAlignReport.eof()) { // need to copy header
869 char c = inputAlignReport.get();
870 goodAlignReportOut << c;
871 if (c == 10 || c == 13){ break; }
874 while(!inputAlignReport.eof()){
875 if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); return 0; }
877 inputAlignReport >> seqName;
878 it = badSeqNames.find(seqName);
880 while (!inputAlignReport.eof()) { // need to copy header
881 char c = inputAlignReport.get();
883 if (c == 10 || c == 13){ break; }
886 if(it != badSeqNames.end()){
887 badSeqNames.erase(it);
890 goodAlignReportOut << seqName << '\t' << line;
892 m->gobble(inputAlignReport);
895 if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); return 0; }
897 //we were unable to remove some of the bad sequences
898 if (badSeqNames.size() != 0) {
899 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
900 m->mothurOut("Your alignreport file does not include the sequence " + *it + " please correct.");
901 m->mothurOutEndLine();
905 inputAlignReport.close();
906 goodAlignReportOut.close();
908 if (m->control_pressed) { remove(goodAlignReportFile.c_str()); return 0; }
913 catch(exception& e) {
914 m->errorOut(e, "ScreenSeqsCommand", "screenAlignReport");
919 //**********************************************************************************************************************
921 int ScreenSeqsCommand::driver(linePair* filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
924 m->openOutputFile(goodFName, goodFile);
926 ofstream badAccnosFile;
927 m->openOutputFile(badAccnosFName, badAccnosFile);
930 m->openInputFile(filename, inFASTA);
932 inFASTA.seekg(filePos->start);
939 if (m->control_pressed) { return 0; }
941 Sequence currSeq(inFASTA); m->gobble(inFASTA);
942 if (currSeq.getName() != "") {
943 bool goodSeq = 1; // innocent until proven guilty
944 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos()) { goodSeq = 0; }
945 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos()) { goodSeq = 0; }
946 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases()) { goodSeq = 0; }
947 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer()) { goodSeq = 0; }
948 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases()) { goodSeq = 0; }
949 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases()) { goodSeq = 0; }
952 currSeq.printSequence(goodFile);
955 badAccnosFile << currSeq.getName() << endl;
956 badSeqNames.insert(currSeq.getName());
961 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
962 unsigned long int pos = inFASTA.tellg();
963 if ((pos == -1) || (pos >= filePos->end)) { break; }
965 if (inFASTA.eof()) { break; }
969 if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
972 if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
977 badAccnosFile.close();
981 catch(exception& e) {
982 m->errorOut(e, "ScreenSeqsCommand", "driver");
986 //**********************************************************************************************************************
988 int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badAccnosFile, vector<unsigned long int>& MPIPos, set<string>& badSeqNames){
990 string outputString = "";
991 MPI_Status statusGood;
992 MPI_Status statusBadAccnos;
995 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
997 for(int i=0;i<num;i++){
999 if (m->control_pressed) { return 0; }
1001 //read next sequence
1002 int length = MPIPos[start+i+1] - MPIPos[start+i];
1004 char* buf4 = new char[length];
1005 memcpy(buf4, outputString.c_str(), length);
1007 MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
1009 string tempBuf = buf4; delete buf4;
1010 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
1011 istringstream iss (tempBuf,istringstream::in);
1013 Sequence currSeq(iss);
1016 if (currSeq.getName() != "") {
1017 bool goodSeq = 1; // innocent until proven guilty
1018 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos()) { goodSeq = 0; }
1019 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos()) { goodSeq = 0; }
1020 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases()) { goodSeq = 0; }
1021 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer()) { goodSeq = 0; }
1022 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases()) { goodSeq = 0; }
1023 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases()) { goodSeq = 0; }
1026 outputString = ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
1029 length = outputString.length();
1030 char* buf2 = new char[length];
1031 memcpy(buf2, outputString.c_str(), length);
1033 MPI_File_write_shared(goodFile, buf2, length, MPI_CHAR, &statusGood);
1038 badSeqNames.insert(currSeq.getName());
1040 //write to bad accnos file
1041 outputString = currSeq.getName() + "\n";
1043 length = outputString.length();
1044 char* buf3 = new char[length];
1045 memcpy(buf3, outputString.c_str(), length);
1047 MPI_File_write_shared(badAccnosFile, buf3, length, MPI_CHAR, &statusBadAccnos);
1053 if((i) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(i)); m->mothurOutEndLine(); }
1058 catch(exception& e) {
1059 m->errorOut(e, "ScreenSeqsCommand", "driverMPI");
1064 /**************************************************************************************************/
1066 int ScreenSeqsCommand::createProcesses(string goodFileName, string badAccnos, string filename, set<string>& badSeqNames) {
1068 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1072 //loop through and create all the processes you want
1073 while (process != processors) {
1077 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1079 }else if (pid == 0){
1080 num = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
1082 //pass numSeqs to parent
1084 string tempFile = filename + toString(getpid()) + ".num.temp";
1085 m->openOutputFile(tempFile, out);
1091 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1092 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1097 //force parent to wait until all the processes are done
1098 for (int i=0;i<processors;i++) {
1099 int temp = processIDS[i];
1103 for (int i = 0; i < processIDS.size(); i++) {
1105 string tempFile = filename + toString(processIDS[i]) + ".num.temp";
1106 m->openInputFile(tempFile, in);
1107 if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
1108 in.close(); remove(tempFile.c_str());
1114 catch(exception& e) {
1115 m->errorOut(e, "ScreenSeqsCommand", "createProcesses");
1120 //***************************************************************************************************************