#include "screenseqscommand.h"
#include "sequence.hpp"
+//**********************************************************************************************************************
+vector<string> ScreenSeqsCommand::getValidParameters(){
+ try {
+ string Array[] = {"fasta", "start", "end", "maxambig", "maxhomop","optimize","criteria", "minlength", "maxlength",
+ "name", "group", "alignreport","processors","outputdir","inputdir"};
+ vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ return myArray;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "getValidParameters");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+ScreenSeqsCommand::ScreenSeqsCommand(){
+ try {
+ abort = true;
+ //initialize outputTypes
+ vector<string> tempOutNames;
+ outputTypes["fasta"] = tempOutNames;
+ outputTypes["name"] = tempOutNames;
+ outputTypes["group"] = tempOutNames;
+ outputTypes["alignreport"] = tempOutNames;
+ outputTypes["accnos"] = tempOutNames;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+vector<string> ScreenSeqsCommand::getRequiredParameters(){
+ try {
+ string Array[] = {"fasta"};
+ vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ return myArray;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "getRequiredParameters");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+vector<string> ScreenSeqsCommand::getRequiredFiles(){
+ try {
+ vector<string> myArray;
+ return myArray;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "getRequiredFiles");
+ exit(1);
+ }
+}
//***************************************************************************************************************
ScreenSeqsCommand::ScreenSeqsCommand(string option) {
else {
//valid paramters for this command
- string AlignArray[] = {"fasta", "start", "end", "maxambig", "maxhomop", "minlength", "maxlength",
+ string AlignArray[] = {"fasta", "start", "end", "maxambig", "maxhomop","optimize","criteria", "minlength", "maxlength",
"name", "group", "alignreport","processors","outputdir","inputdir"};
vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
OptionParser parser(option);
map<string,string> parameters = parser.getParameters();
- ValidParameters validParameter;
+ ValidParameters validParameter("screen.seqs");
map<string,string>::iterator it;
//check to make sure all parameters are valid for command
if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
}
+ //initialize outputTypes
+ vector<string> tempOutNames;
+ outputTypes["fasta"] = tempOutNames;
+ outputTypes["name"] = tempOutNames;
+ outputTypes["group"] = tempOutNames;
+ outputTypes["alignreport"] = tempOutNames;
+ outputTypes["accnos"] = tempOutNames;
+
//if the user changes the input directory command factory will send this info to us in the output parameter
string inputDir = validParameter.validFile(parameters, "inputdir", false);
if (inputDir == "not found"){ inputDir = ""; }
it = parameters.find("fasta");
//user has given a template file
if(it != parameters.end()){
- path = hasPath(it->second);
+ path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["fasta"] = inputDir + it->second; }
}
it = parameters.find("group");
//user has given a template file
if(it != parameters.end()){
- path = hasPath(it->second);
+ path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["group"] = inputDir + it->second; }
}
it = parameters.find("name");
//user has given a template file
if(it != parameters.end()){
- path = hasPath(it->second);
+ path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["name"] = inputDir + it->second; }
}
it = parameters.find("alignreport");
//user has given a template file
if(it != parameters.end()){
- path = hasPath(it->second);
+ path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["alignreport"] = inputDir + it->second; }
}
//if the user changes the output directory command factory will send this info to us in the output parameter
outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
outputDir = "";
- outputDir += hasPath(fastafile); //if user entered a file with a path then preserve it
+ outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it
}
//check for optional parameter and set defaults
temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
convert(temp, processors);
-
+
+ temp = validParameter.validFile(parameters, "optimize", false); //optimizing trumps the optimized values original value
+ if (temp == "not found"){ temp = ""; }
+ else { m->splitAtDash(temp, optimize); }
+
+ //check for invalid optimize options
+ set<string> validOptimizers;
+ validOptimizers.insert("start"); validOptimizers.insert("end"); validOptimizers.insert("maxambig"); validOptimizers.insert("maxhomop"); validOptimizers.insert("minlength"); validOptimizers.insert("maxlength");
+
+ for (int i = 0; i < optimize.size(); i++) {
+ if (validOptimizers.count(optimize[i]) == 0) {
+ m->mothurOut(optimize[i] + " is not a valid optimizer. Valid options are start, end, maxambig, maxhomop, minlength and maxlength."); m->mothurOutEndLine();
+ optimize.erase(optimize.begin()+i);
+ i--;
+ }
+ }
+
+ temp = validParameter.validFile(parameters, "criteria", false); if (temp == "not found"){ temp = "90"; }
+ convert(temp, criteria);
}
}
void ScreenSeqsCommand::help(){
try {
m->mothurOut("The screen.seqs command reads a fastafile and creates .....\n");
- m->mothurOut("The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, and group.\n");
+ m->mothurOut("The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group, optimize, criteria and processors.\n");
m->mothurOut("The fasta parameter is required.\n");
m->mothurOut("The start parameter .... The default is -1.\n");
m->mothurOut("The end parameter .... The default is -1.\n");
m->mothurOut("The maxhomop parameter .... The default is -1.\n");
m->mothurOut("The minlength parameter .... The default is -1.\n");
m->mothurOut("The maxlength parameter .... The default is -1.\n");
+ m->mothurOut("The processors parameter allows you to specify the number of processors to use while running the command. The default is 1.\n");
+ 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");
+ 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");
m->mothurOut("The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n");
m->mothurOut("The screen.seqs command should be in the following format: \n");
m->mothurOut("screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig, \n");
try{
if (abort == true) { return 0; }
+
+ //if the user want to optimize we need to no the 90% mark
+ vector<unsigned long int> positions;
+ if (optimize.size() != 0) { getSummary(positions); } //get summary is paralellized so we need to divideFile, no need to do this step twice so I moved it here
+ else {
+ positions = m->divideFile(fastafile, processors);
+ for (int i = 0; i < (positions.size()-1); i++) {
+ lines.push_back(new linePair(positions[i], positions[(i+1)]));
+ }
+ }
- string goodSeqFile = outputDir + getRootName(getSimpleName(fastafile)) + "good" + getExtension(fastafile);
- string badSeqFile = outputDir + getRootName(getSimpleName(fastafile)) + "bad" + getExtension(fastafile);
- string badAccnosFile = outputDir + getRootName(getSimpleName(fastafile)) + "bad.accnos";
+ string goodSeqFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "good" + m->getExtension(fastafile);
+ string badAccnosFile = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "bad.accnos";
int numFastaSeqs = 0;
set<string> badSeqNames;
#ifdef USE_MPI
int pid, end, numSeqsPerProcessor;
int tag = 2001;
- vector<long> MPIPos;
+ vector<unsigned long int> MPIPos;
MPI_Status status;
MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
MPI_File inMPI;
MPI_File outMPIGood;
- MPI_File outMPIBad;
MPI_File outMPIBadAccnos;
int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
char outGoodFilename[1024];
strcpy(outGoodFilename, goodSeqFile.c_str());
- char outBadFilename[1024];
- strcpy(outBadFilename, badSeqFile.c_str());
-
char outBadAccnosFilename[1024];
strcpy(outBadAccnosFilename, badAccnosFile.c_str());
MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI); //comm, filename, mode, info, filepointer
MPI_File_open(MPI_COMM_WORLD, outGoodFilename, outMode, MPI_INFO_NULL, &outMPIGood);
- MPI_File_open(MPI_COMM_WORLD, outBadFilename, outMode, MPI_INFO_NULL, &outMPIBad);
MPI_File_open(MPI_COMM_WORLD, outBadAccnosFilename, outMode, MPI_INFO_NULL, &outMPIBadAccnos);
- if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBad); MPI_File_close(&outMPIBadAccnos); return 0; }
+ if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
if (pid == 0) { //you are the root process
- MPIPos = setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
+ MPIPos = m->setFilePosFasta(fastafile, numFastaSeqs); //fills MPIPos, returns numSeqs
//send file positions to all processes
- MPI_Bcast(&numFastaSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD); //send numSeqs
- MPI_Bcast(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos
+ for(int i = 1; i < processors; i++) {
+ MPI_Send(&numFastaSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
+ MPI_Send(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
+ }
//figure out how many sequences you have to align
numSeqsPerProcessor = numFastaSeqs / processors;
int startIndex = pid * numSeqsPerProcessor;
if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
-
+ // cout << pid << '\t' << numSeqsPerProcessor << '\t' << startIndex << endl;
//align your part
- driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBad, outMPIBadAccnos, MPIPos, badSeqNames);
-
- if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); MPI_File_close(&outMPIBad); return 0; }
+ driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
+ //cout << pid << " done" << endl;
+ if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
for (int i = 1; i < processors; i++) {
}*/
}
}else{ //you are a child process
- MPI_Bcast(&numFastaSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs
+ MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
MPIPos.resize(numFastaSeqs+1);
- MPI_Bcast(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions
-
+ MPI_Recv(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
+
//figure out how many sequences you have to align
numSeqsPerProcessor = numFastaSeqs / processors;
int startIndex = pid * numSeqsPerProcessor;
if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; }
-
+ //cout << pid << '\t' << numSeqsPerProcessor << '\t' << startIndex << endl;
//align your part
- driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBad, outMPIBadAccnos, MPIPos, badSeqNames);
-
- if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBad); MPI_File_close(&outMPIBadAccnos); return 0; }
+ driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames);
+//cout << pid << " done" << endl;
+ if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; }
//send bad list
int badSize = badSeqNames.size();
//close files
MPI_File_close(&inMPI);
MPI_File_close(&outMPIGood);
- MPI_File_close(&outMPIBad);
MPI_File_close(&outMPIBadAccnos);
+ MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
#else
-
+
#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
if(processors == 1){
- ifstream inFASTA;
- openInputFile(fastafile, inFASTA);
- numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
- inFASTA.close();
+ numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
- lines.push_back(new linePair(0, numFastaSeqs));
-
- driver(lines[0], goodSeqFile, badSeqFile, badAccnosFile, fastafile, badSeqNames);
-
- if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
}else{
- vector<int> positions;
processIDS.resize(0);
- ifstream inFASTA;
- openInputFile(fastafile, inFASTA);
-
- string input;
- while(!inFASTA.eof()){
- input = getline(inFASTA);
- if (input.length() != 0) {
- if(input[0] == '>'){ long int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1); }
- }
- }
- inFASTA.close();
-
- numFastaSeqs = positions.size();
-
- int numSeqsPerProcessor = numFastaSeqs / processors;
-
- for (int i = 0; i < processors; i++) {
- long int startPos = positions[ i * numSeqsPerProcessor ];
- if(i == processors - 1){
- numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
- }
- lines.push_back(new linePair(startPos, numSeqsPerProcessor));
- }
-
- createProcesses(goodSeqFile, badSeqFile, badAccnosFile, fastafile, badSeqNames);
+ numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames);
rename((goodSeqFile + toString(processIDS[0]) + ".temp").c_str(), goodSeqFile.c_str());
- rename((badSeqFile + toString(processIDS[0]) + ".temp").c_str(), badSeqFile.c_str());
rename((badAccnosFile + toString(processIDS[0]) + ".temp").c_str(), badAccnosFile.c_str());
//append alignment and report files
for(int i=1;i<processors;i++){
- appendFiles((goodSeqFile + toString(processIDS[i]) + ".temp"), goodSeqFile);
+ m->appendFiles((goodSeqFile + toString(processIDS[i]) + ".temp"), goodSeqFile);
remove((goodSeqFile + toString(processIDS[i]) + ".temp").c_str());
-
- appendFiles((badSeqFile + toString(processIDS[i]) + ".temp"), badSeqFile);
- remove((badSeqFile + toString(processIDS[i]) + ".temp").c_str());
-
- appendFiles((badAccnosFile + toString(processIDS[i]) + ".temp"), badAccnosFile);
+
+ m->appendFiles((badAccnosFile + toString(processIDS[i]) + ".temp"), badAccnosFile);
remove((badAccnosFile + toString(processIDS[i]) + ".temp").c_str());
}
- if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
//read badSeqs in because root process doesnt know what other "bad" seqs the children found
ifstream inBad;
- int ableToOpen = openInputFile(badAccnosFile, inBad, "no error");
+ int ableToOpen = m->openInputFile(badAccnosFile, inBad, "no error");
if (ableToOpen == 0) {
badSeqNames.clear();
string tempName;
while (!inBad.eof()) {
- inBad >> tempName; gobble(inBad);
+ inBad >> tempName; m->gobble(inBad);
badSeqNames.insert(tempName);
}
inBad.close();
}
}
#else
- ifstream inFASTA;
- openInputFile(fastafile, inFASTA);
- numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
- inFASTA.close();
-
- lines.push_back(new linePair(0, numFastaSeqs));
-
- driver(lines[0], goodSeqFile, badSeqFile, badAccnosFile, fastafile, badSeqNames);
+ numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
- if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
#endif
badSeqNames.clear();
string tempName;
while (!iss.eof()) {
- iss >> tempName; gobble(iss);
+ iss >> tempName; m->gobble(iss);
badSeqNames.insert(tempName);
}
#endif
if(namefile != "" && groupfile != "") {
screenNameGroupFile(badSeqNames);
- if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
}else if(namefile != "") {
screenNameGroupFile(badSeqNames);
- if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
}else if(groupfile != "") { screenGroupFile(badSeqNames); } // this screens just the group
- if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
if(alignreport != "") { screenAlignReport(badSeqNames); }
- if (m->control_pressed) { remove(goodSeqFile.c_str()); remove(badSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
#ifdef USE_MPI
}
m->mothurOutEndLine();
m->mothurOut("Output File Names: "); m->mothurOutEndLine();
- m->mothurOut(goodSeqFile); m->mothurOutEndLine();
- m->mothurOut(badSeqFile); m->mothurOutEndLine();
- m->mothurOut(badAccnosFile); m->mothurOutEndLine();
+ m->mothurOut(goodSeqFile); m->mothurOutEndLine(); outputTypes["fasta"].push_back(goodSeqFile);
+ m->mothurOut(badAccnosFile); m->mothurOutEndLine(); outputTypes["accnos"].push_back(badAccnosFile);
for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
m->mothurOutEndLine();
m->mothurOutEndLine();
int ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
try {
ifstream inputNames;
- openInputFile(namefile, inputNames);
+ m->openInputFile(namefile, inputNames);
set<string> badSeqGroups;
string seqName, seqList, group;
set<string>::iterator it;
- string goodNameFile = outputDir + getRootName(getSimpleName(namefile)) + "good" + getExtension(namefile);
- string badNameFile = outputDir + getRootName(getSimpleName(namefile)) + "bad" + getExtension(namefile);
-
- outputNames.push_back(goodNameFile); outputNames.push_back(badNameFile);
-
- ofstream goodNameOut; openOutputFile(goodNameFile, goodNameOut);
- ofstream badNameOut; openOutputFile(badNameFile, badNameOut);
+ string goodNameFile = outputDir + m->getRootName(m->getSimpleName(namefile)) + "good" + m->getExtension(namefile);
+ outputNames.push_back(goodNameFile); outputTypes["name"].push_back(goodNameFile);
+ ofstream goodNameOut; m->openOutputFile(goodNameFile, goodNameOut);
+
while(!inputNames.eof()){
- if (m->control_pressed) { goodNameOut.close(); badNameOut.close(); inputNames.close(); remove(goodNameFile.c_str()); remove(badNameFile.c_str()); return 0; }
+ if (m->control_pressed) { goodNameOut.close(); inputNames.close(); remove(goodNameFile.c_str()); return 0; }
inputNames >> seqName >> seqList;
it = badSeqNames.find(seqName);
if(it != badSeqNames.end()){
badSeqNames.erase(it);
- badNameOut << seqName << '\t' << seqList << endl;
+
if(namefile != ""){
int start = 0;
for(int i=0;i<seqList.length();i++){
else{
goodNameOut << seqName << '\t' << seqList << endl;
}
- gobble(inputNames);
+ m->gobble(inputNames);
}
inputNames.close();
goodNameOut.close();
- badNameOut.close();
-
+
//we were unable to remove some of the bad sequences
if (badSeqNames.size() != 0) {
for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
if(groupfile != ""){
ifstream inputGroups;
- openInputFile(groupfile, inputGroups);
+ m->openInputFile(groupfile, inputGroups);
- string goodGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "good" + getExtension(groupfile);
- string badGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "bad" + getExtension(groupfile);
+ string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
+ outputNames.push_back(goodGroupFile); outputTypes["group"].push_back(goodGroupFile);
- outputNames.push_back(goodGroupFile); outputNames.push_back(badGroupFile);
-
- ofstream goodGroupOut; openOutputFile(goodGroupFile, goodGroupOut);
- ofstream badGroupOut; openOutputFile(badGroupFile, badGroupOut);
+ ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
while(!inputGroups.eof()){
- 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; }
+ if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodNameFile.c_str()); remove(goodGroupFile.c_str()); return 0; }
inputGroups >> seqName >> group;
if(it != badSeqGroups.end()){
badSeqGroups.erase(it);
- badGroupOut << seqName << '\t' << group << endl;
}
else{
goodGroupOut << seqName << '\t' << group << endl;
}
- gobble(inputGroups);
+ m->gobble(inputGroups);
}
inputGroups.close();
goodGroupOut.close();
- badGroupOut.close();
//we were unable to remove some of the bad sequences
if (badSeqGroups.size() != 0) {
for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {
- m->mothurOut("Your namefile does not include the sequence " + *it + " please correct.");
+ m->mothurOut("Your groupfile does not include the sequence " + *it + " please correct.");
m->mothurOutEndLine();
}
}
exit(1);
}
}
+//***************************************************************************************************************
+int ScreenSeqsCommand::getSummary(vector<unsigned long int>& positions){
+ try {
+
+ vector<int> startPosition;
+ vector<int> endPosition;
+ vector<int> seqLength;
+ vector<int> ambigBases;
+ vector<int> longHomoPolymer;
+
+ vector<unsigned long int> positions = m->divideFile(fastafile, processors);
+
+ for (int i = 0; i < (positions.size()-1); i++) {
+ lines.push_back(new linePair(positions[i], positions[(i+1)]));
+ }
+
+ int numSeqs = 0;
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ if(processors == 1){
+ numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
+ }else{
+ numSeqs = createProcessesCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile);
+ }
+
+ if (m->control_pressed) { return 0; }
+ #else
+ numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
+ if (m->control_pressed) { return 0; }
+ #endif
+
+ sort(startPosition.begin(), startPosition.end());
+ sort(endPosition.begin(), endPosition.end());
+ sort(seqLength.begin(), seqLength.end());
+ sort(ambigBases.begin(), ambigBases.end());
+ sort(longHomoPolymer.begin(), longHomoPolymer.end());
+
+ int criteriaPercentile = int(numSeqs * (criteria / (float) 100));
+
+ for (int i = 0; i < optimize.size(); i++) {
+ if (optimize[i] == "start") { startPos = startPosition[criteriaPercentile]; m->mothurOut("Optimizing start to " + toString(startPos) + "."); m->mothurOutEndLine(); }
+ else if (optimize[i] == "end") { int endcriteriaPercentile = int(numSeqs * ((100 - criteria) / (float) 100)); endPos = endPosition[endcriteriaPercentile]; m->mothurOut("Optimizing end to " + toString(endPos) + "."); m->mothurOutEndLine();}
+ else if (optimize[i] == "maxambig") { maxAmbig = ambigBases[criteriaPercentile]; m->mothurOut("Optimizing maxambig to " + toString(maxAmbig) + "."); m->mothurOutEndLine(); }
+ else if (optimize[i] == "maxhomop") { maxHomoP = longHomoPolymer[criteriaPercentile]; m->mothurOut("Optimizing maxhomop to " + toString(maxHomoP) + "."); m->mothurOutEndLine(); }
+ else if (optimize[i] == "minlength") { int mincriteriaPercentile = int(numSeqs * ((100 - criteria) / (float) 100)); minLength = seqLength[mincriteriaPercentile]; m->mothurOut("Optimizing minlength to " + toString(minLength) + "."); m->mothurOutEndLine(); }
+ else if (optimize[i] == "maxlength") { maxLength = seqLength[criteriaPercentile]; m->mothurOut("Optimizing maxlength to " + toString(maxLength) + "."); m->mothurOutEndLine(); }
+ }
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "getSummary");
+ exit(1);
+ }
+}
+/**************************************************************************************/
+int ScreenSeqsCommand::driverCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, linePair* filePos) {
+ try {
+
+ ifstream in;
+ m->openInputFile(filename, in);
+
+ in.seekg(filePos->start);
+
+ bool done = false;
+ int count = 0;
+
+ while (!done) {
+
+ if (m->control_pressed) { in.close(); return 1; }
+
+ Sequence current(in); m->gobble(in);
+
+ if (current.getName() != "") {
+ startPosition.push_back(current.getStartPos());
+ endPosition.push_back(current.getEndPos());
+ seqLength.push_back(current.getNumBases());
+ ambigBases.push_back(current.getAmbigBases());
+ longHomoPolymer.push_back(current.getLongHomoPolymer());
+ count++;
+ }
+
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ unsigned long int pos = in.tellg();
+ if ((pos == -1) || (pos >= filePos->end)) { break; }
+ #else
+ if (in.eof()) { break; }
+ #endif
+
+ }
+
+ in.close();
+
+ return count;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "driverCreateSummary");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+int ScreenSeqsCommand::createProcessesCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename) {
+ try {
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ int process = 1;
+ int num = 0;
+ processIDS.clear();
+
+ //loop through and create all the processes you want
+ while (process != processors) {
+ int pid = fork();
+
+ if (pid > 0) {
+ processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
+ process++;
+ }else if (pid == 0){
+ num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[process]);
+
+ //pass numSeqs to parent
+ ofstream out;
+ string tempFile = fastafile + toString(getpid()) + ".num.temp";
+ m->openOutputFile(tempFile, out);
+
+ out << num << endl;
+ for (int k = 0; k < startPosition.size(); k++) { out << startPosition[k] << '\t'; } out << endl;
+ for (int k = 0; k < endPosition.size(); k++) { out << endPosition[k] << '\t'; } out << endl;
+ for (int k = 0; k < seqLength.size(); k++) { out << seqLength[k] << '\t'; } out << endl;
+ for (int k = 0; k < ambigBases.size(); k++) { out << ambigBases[k] << '\t'; } out << endl;
+ for (int k = 0; k < longHomoPolymer.size(); k++) { out << longHomoPolymer[k] << '\t'; } out << endl;
+
+ out.close();
+
+ exit(0);
+ }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
+ }
+
+ num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
+
+ //force parent to wait until all the processes are done
+ for (int i=0;i<processIDS.size();i++) {
+ int temp = processIDS[i];
+ wait(&temp);
+ }
+
+ //parent reads in and combine Filter info
+ for (int i = 0; i < processIDS.size(); i++) {
+ string tempFilename = fastafile + toString(processIDS[i]) + ".num.temp";
+ ifstream in;
+ m->openInputFile(tempFilename, in);
+
+ int temp, tempNum;
+ in >> tempNum; m->gobble(in); num += tempNum;
+ for (int k = 0; k < tempNum; k++) { in >> temp; startPosition.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; endPosition.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; seqLength.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; ambigBases.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; longHomoPolymer.push_back(temp); } m->gobble(in);
+
+ in.close();
+ remove(tempFilename.c_str());
+ }
+
+ return num;
+#endif
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "createProcessesCreateSummary");
+ exit(1);
+ }
+}
//***************************************************************************************************************
int ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
try {
ifstream inputGroups;
- openInputFile(groupfile, inputGroups);
+ m->openInputFile(groupfile, inputGroups);
string seqName, group;
set<string>::iterator it;
- string goodGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "good" + getExtension(groupfile);
- string badGroupFile = outputDir + getRootName(getSimpleName(groupfile)) + "bad" + getExtension(groupfile);
-
- outputNames.push_back(goodGroupFile); outputNames.push_back(badGroupFile);
-
- ofstream goodGroupOut; openOutputFile(goodGroupFile, goodGroupOut);
- ofstream badGroupOut; openOutputFile(badGroupFile, badGroupOut);
+ string goodGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "good" + m->getExtension(groupfile);
+ outputNames.push_back(goodGroupFile); outputTypes["group"].push_back(goodGroupFile);
+ ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
while(!inputGroups.eof()){
- if (m->control_pressed) { goodGroupOut.close(); badGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); remove(badGroupFile.c_str()); return 0; }
+ if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); return 0; }
inputGroups >> seqName >> group;
it = badSeqNames.find(seqName);
if(it != badSeqNames.end()){
badSeqNames.erase(it);
- badGroupOut << seqName << '\t' << group << endl;
}
else{
goodGroupOut << seqName << '\t' << group << endl;
}
- gobble(inputGroups);
+ m->gobble(inputGroups);
}
- if (m->control_pressed) { goodGroupOut.close(); badGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); remove(badGroupFile.c_str()); return 0; }
+ if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); return 0; }
//we were unable to remove some of the bad sequences
if (badSeqNames.size() != 0) {
inputGroups.close();
goodGroupOut.close();
- badGroupOut.close();
- if (m->control_pressed) { remove(goodGroupFile.c_str()); remove(badGroupFile.c_str()); }
-
+ if (m->control_pressed) { remove(goodGroupFile.c_str()); }
return 0;
int ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
try {
ifstream inputAlignReport;
- openInputFile(alignreport, inputAlignReport);
+ m->openInputFile(alignreport, inputAlignReport);
string seqName, group;
set<string>::iterator it;
- string goodAlignReportFile = outputDir + getRootName(getSimpleName(alignreport)) + "good" + getExtension(alignreport);
- string badAlignReportFile = outputDir + getRootName(getSimpleName(alignreport)) + "bad" + getExtension(alignreport);
-
- outputNames.push_back(goodAlignReportFile); outputNames.push_back(badAlignReportFile);
-
- ofstream goodAlignReportOut; openOutputFile(goodAlignReportFile, goodAlignReportOut);
- ofstream badAlignReportOut; openOutputFile(badAlignReportFile, badAlignReportOut);
+ string goodAlignReportFile = outputDir + m->getRootName(m->getSimpleName(alignreport)) + "good" + m->getExtension(alignreport);
+ outputNames.push_back(goodAlignReportFile); outputTypes["alignreport"].push_back(goodAlignReportFile);
+ ofstream goodAlignReportOut; m->openOutputFile(goodAlignReportFile, goodAlignReportOut);
while (!inputAlignReport.eof()) { // need to copy header
char c = inputAlignReport.get();
goodAlignReportOut << c;
- badAlignReportOut << c;
if (c == 10 || c == 13){ break; }
}
while(!inputAlignReport.eof()){
- if (m->control_pressed) { goodAlignReportOut.close(); badAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); remove(badAlignReportFile.c_str()); return 0; }
+ if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); return 0; }
inputAlignReport >> seqName;
it = badSeqNames.find(seqName);
if(it != badSeqNames.end()){
badSeqNames.erase(it);
- badAlignReportOut << seqName << '\t' << line;
}
else{
goodAlignReportOut << seqName << '\t' << line;
}
- gobble(inputAlignReport);
+ m->gobble(inputAlignReport);
}
- if (m->control_pressed) { goodAlignReportOut.close(); badAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); remove(badAlignReportFile.c_str()); return 0; }
+ if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); return 0; }
//we were unable to remove some of the bad sequences
if (badSeqNames.size() != 0) {
for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
- m->mothurOut("Your file does not include the sequence " + *it + " please correct.");
+ m->mothurOut("Your alignreport file does not include the sequence " + *it + " please correct.");
m->mothurOutEndLine();
}
}
inputAlignReport.close();
goodAlignReportOut.close();
- badAlignReportOut.close();
- if (m->control_pressed) { remove(goodAlignReportFile.c_str()); remove(badAlignReportFile.c_str()); return 0; }
+ if (m->control_pressed) { remove(goodAlignReportFile.c_str()); return 0; }
return 0;
}
//**********************************************************************************************************************
-int ScreenSeqsCommand::driver(linePair* line, string goodFName, string badFName, string badAccnosFName, string filename, set<string>& badSeqNames){
+int ScreenSeqsCommand::driver(linePair* filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
try {
ofstream goodFile;
- openOutputFile(goodFName, goodFile);
-
- ofstream badFile;
- openOutputFile(badFName, badFile);
+ m->openOutputFile(goodFName, goodFile);
ofstream badAccnosFile;
- openOutputFile(badAccnosFName, badAccnosFile);
+ m->openOutputFile(badAccnosFName, badAccnosFile);
ifstream inFASTA;
- openInputFile(filename, inFASTA);
+ m->openInputFile(filename, inFASTA);
+
+ inFASTA.seekg(filePos->start);
- inFASTA.seekg(line->start);
+ bool done = false;
+ int count = 0;
- for(int i=0;i<line->numSeqs;i++){
-
+ while (!done) {
+
if (m->control_pressed) { return 0; }
- Sequence currSeq(inFASTA);
+ Sequence currSeq(inFASTA); m->gobble(inFASTA);
if (currSeq.getName() != "") {
bool goodSeq = 1; // innocent until proven guilty
if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos()) { goodSeq = 0; }
currSeq.printSequence(goodFile);
}
else{
- currSeq.printSequence(badFile);
badAccnosFile << currSeq.getName() << endl;
badSeqNames.insert(currSeq.getName());
}
+ count++;
}
- gobble(inFASTA);
+
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ unsigned long int pos = inFASTA.tellg();
+ if ((pos == -1) || (pos >= filePos->end)) { break; }
+ #else
+ if (inFASTA.eof()) { break; }
+ #endif
+
+ //report progress
+ if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
}
+ //report progress
+ if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); }
goodFile.close();
inFASTA.close();
- badFile.close();
badAccnosFile.close();
- return 1;
+ return count;
}
catch(exception& e) {
m->errorOut(e, "ScreenSeqsCommand", "driver");
}
//**********************************************************************************************************************
#ifdef USE_MPI
-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){
+int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badAccnosFile, vector<unsigned long int>& MPIPos, set<string>& badSeqNames){
try {
string outputString = "";
MPI_Status statusGood;
- MPI_Status statusBad;
MPI_Status statusBadAccnos;
MPI_Status status;
int pid;
delete buf2;
}
else{
- outputString = ">" + currSeq.getName() + "\n" + currSeq.getAligned() + "\n";
-
- //print bad seq to fasta
- length = outputString.length();
- char* buf2 = new char[length];
- memcpy(buf2, outputString.c_str(), length);
-
- MPI_File_write_shared(badFile, buf2, length, MPI_CHAR, &statusBad);
- delete buf2;
badSeqNames.insert(currSeq.getName());
#endif
/**************************************************************************************************/
-int ScreenSeqsCommand::createProcesses(string goodFileName, string badFileName, string badAccnos, string filename, set<string>& badSeqNames) {
+int ScreenSeqsCommand::createProcesses(string goodFileName, string badAccnos, string filename, set<string>& badSeqNames) {
try {
#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
int process = 0;
- int exitCommand = 1;
+ int num = 0;
//loop through and create all the processes you want
while (process != processors) {
processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
process++;
}else if (pid == 0){
- exitCommand = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
+ num = driver(lines[process], goodFileName + toString(getpid()) + ".temp", badAccnos + toString(getpid()) + ".temp", filename, badSeqNames);
+
+ //pass numSeqs to parent
+ ofstream out;
+ string tempFile = filename + toString(getpid()) + ".num.temp";
+ m->openOutputFile(tempFile, out);
+ out << num << endl;
+ out.close();
+
exit(0);
}else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
}
wait(&temp);
}
- return exitCommand;
+ for (int i = 0; i < processIDS.size(); i++) {
+ ifstream in;
+ string tempFile = filename + toString(processIDS[i]) + ".num.temp";
+ m->openInputFile(tempFile, in);
+ if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
+ in.close(); remove(tempFile.c_str());
+ }
+
+ return num;
#endif
}
catch(exception& e) {
projects / mothur.git / blobdiff