*/
#include "screenseqscommand.h"
-#include "sequence.hpp"
+
//**********************************************************************************************************************
vector<string> ScreenSeqsCommand::setParameters(){
CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
+ CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pqfile);
CommandParameter palignreport("alignreport", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(palignreport);
+ CommandParameter ptax("taxonomy", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(ptax);
CommandParameter pstart("start", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pstart);
CommandParameter pend("end", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pend);
CommandParameter pmaxambig("maxambig", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pmaxambig);
try {
string helpString = "";
helpString += "The screen.seqs command reads a fastafile and creates .....\n";
- helpString += "The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group, optimize, criteria and processors.\n";
+ helpString += "The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, group, qfile, alignreport, taxonomy, optimize, criteria and processors.\n";
helpString += "The fasta parameter is required.\n";
+ helpString += "The alignreport and taxonomy parameters allow you to remove bad seqs from taxonomy and alignreport files.\n";
helpString += "The start parameter .... The default is -1.\n";
helpString += "The end parameter .... The default is -1.\n";
helpString += "The maxambig parameter allows you to set the maximum number of ambigious bases allowed. The default is -1.\n";
helpString += "screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig, \n";
helpString += "maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n";
helpString += "Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n";
- helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n";
+ helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
return helpString;
}
catch(exception& e) {
outputTypes["group"] = tempOutNames;
outputTypes["alignreport"] = tempOutNames;
outputTypes["accnos"] = tempOutNames;
+ outputTypes["qfile"] = tempOutNames;
+ outputTypes["taxonomy"] = tempOutNames;
}
catch(exception& e) {
m->errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
//allow user to run help
if(option == "help") { help(); abort = true; calledHelp = true; }
+ else if(option == "citation") { citation(); abort = true; calledHelp = true;}
else {
vector<string> myArray = setParameters();
outputTypes["group"] = tempOutNames;
outputTypes["alignreport"] = tempOutNames;
outputTypes["accnos"] = tempOutNames;
+ outputTypes["qfile"] = tempOutNames;
+ outputTypes["taxonomy"] = 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 the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["alignreport"] = inputDir + it->second; }
}
+
+ it = parameters.find("qfile");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["qfile"] = inputDir + it->second; }
+ }
+
+ it = parameters.find("taxonomy");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["taxonomy"] = inputDir + it->second; }
+ }
}
//check for required parameters
if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
else { m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
}
- else if (fastafile == "not open") { abort = true; }
+ else if (fastafile == "not open") { abort = true; }
+ else { m->setFastaFile(fastafile); }
groupfile = validParameter.validFile(parameters, "group", true);
if (groupfile == "not open") { abort = true; }
else if (groupfile == "not found") { groupfile = ""; }
+ else { m->setGroupFile(groupfile); }
+
+ qualfile = validParameter.validFile(parameters, "qfile", true);
+ if (qualfile == "not open") { abort = true; }
+ else if (qualfile == "not found") { qualfile = ""; }
+ else { m->setQualFile(qualfile); }
namefile = validParameter.validFile(parameters, "name", true);
if (namefile == "not open") { namefile = ""; abort = true; }
else if (namefile == "not found") { namefile = ""; }
-
+ else { m->setNameFile(namefile); }
+
alignreport = validParameter.validFile(parameters, "alignreport", true);
if (alignreport == "not open") { abort = true; }
- else if (alignreport == "not found") { alignreport = ""; }
+ else if (alignreport == "not found") { alignreport = ""; }
+
+ taxonomy = validParameter.validFile(parameters, "taxonomy", true);
+ if (taxonomy == "not open") { abort = true; }
+ else if (taxonomy == "not found") { taxonomy = ""; }
//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"){
// ...at some point should added some additional type checking...
string temp;
temp = validParameter.validFile(parameters, "start", false); if (temp == "not found") { temp = "-1"; }
- convert(temp, startPos);
+ m->mothurConvert(temp, startPos);
temp = validParameter.validFile(parameters, "end", false); if (temp == "not found") { temp = "-1"; }
- convert(temp, endPos);
+ m->mothurConvert(temp, endPos);
temp = validParameter.validFile(parameters, "maxambig", false); if (temp == "not found") { temp = "-1"; }
- convert(temp, maxAmbig);
+ m->mothurConvert(temp, maxAmbig);
temp = validParameter.validFile(parameters, "maxhomop", false); if (temp == "not found") { temp = "-1"; }
- convert(temp, maxHomoP);
+ m->mothurConvert(temp, maxHomoP);
temp = validParameter.validFile(parameters, "minlength", false); if (temp == "not found") { temp = "-1"; }
- convert(temp, minLength);
+ m->mothurConvert(temp, minLength);
temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found") { temp = "-1"; }
- convert(temp, maxLength);
+ m->mothurConvert(temp, maxLength);
temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
m->setProcessors(temp);
- convert(temp, processors);
+ m->mothurConvert(temp, processors);
temp = validParameter.validFile(parameters, "optimize", false); //optimizing trumps the optimized values original value
if (temp == "not found"){ temp = "none"; }
if (optimize.size() == 1) { if (optimize[0] == "none") { optimize.clear(); } }
temp = validParameter.validFile(parameters, "criteria", false); if (temp == "not found"){ temp = "90"; }
- convert(temp, criteria);
+ m->mothurConvert(temp, criteria);
+
+ if (namefile == "") {
+ vector<string> files; files.push_back(fastafile);
+ parser.getNameFile(files);
+ }
}
}
if (abort == true) { if (calledHelp) { return 0; } return 2; }
//if the user want to optimize we need to know the 90% mark
- vector<unsigned long int> positions;
+ vector<unsigned long long> positions;
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
//use the namefile to optimize correctly
if (namefile != "") { nameMap = m->readNames(namefile); }
getSummary(positions);
}
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)]));
- }
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ positions = m->divideFile(fastafile, processors);
+ for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }
+ #else
+ if(processors == 1){ lines.push_back(linePair(0, 1000)); }
+ else {
+ int numFastaSeqs = 0;
+ positions = m->setFilePosFasta(fastafile, numFastaSeqs);
+ if (positions.size() < processors) { processors = positions.size(); }
+
+ //figure out how many sequences you have to process
+ int numSeqsPerProcessor = numFastaSeqs / processors;
+ for (int i = 0; i < processors; i++) {
+ int startIndex = i * numSeqsPerProcessor;
+ if(i == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor; }
+ lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
+ }
+ }
+ #endif
}
-
+
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;
int start = time(NULL);
-
+
#ifdef USE_MPI
int pid, numSeqsPerProcessor;
int tag = 2001;
- vector<unsigned long int> MPIPos;
+ vector<unsigned long long> MPIPos;
MPI_Status status;
MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
MPI_Comm_size(MPI_COMM_WORLD, &processors);
-
+
MPI_File inMPI;
MPI_File outMPIGood;
MPI_File outMPIBadAccnos;
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, 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++) {
-
//get bad lists
int badSize;
MPI_Recv(&badSize, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
- /*for (int j = 0; j < badSize; j++) {
- int length;
- MPI_Recv(&length, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status); //recv the length of the name
- char* buf2 = new char[length]; //make space to recieve it
- MPI_Recv(buf2, length, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); //get name
-
- string tempBuf = buf2;
- if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
- delete buf2;
-
- badSeqNames.insert(tempBuf);
- }*/
}
}else{ //you are a child process
MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
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, 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();
MPI_Send(&badSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
-
- /*
- set<string>::iterator it;
- for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
- string name = *it;
- int length = name.length();
- char* buf2 = new char[length];
- memcpy(buf2, name.c_str(), length);
-
- MPI_Send(&length, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
- MPI_Send(buf2, length, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
- }*/
}
//close files
MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
#else
-
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- if(processors == 1){
- numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
-
- if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
-
- }else{
- processIDS.resize(0);
-
- numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames);
-
- rename((goodSeqFile + toString(processIDS[0]) + ".temp").c_str(), goodSeqFile.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++){
- m->appendFiles((goodSeqFile + toString(processIDS[i]) + ".temp"), goodSeqFile);
- remove((goodSeqFile + toString(processIDS[i]) + ".temp").c_str());
-
- m->appendFiles((badAccnosFile + toString(processIDS[i]) + ".temp"), badAccnosFile);
- remove((badAccnosFile + toString(processIDS[i]) + ".temp").c_str());
- }
-
- 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 = m->openInputFile(badAccnosFile, inBad, "no error");
-
- if (ableToOpen == 0) {
- badSeqNames.clear();
- string tempName;
- while (!inBad.eof()) {
- inBad >> tempName; m->gobble(inBad);
- badSeqNames.insert(tempName);
- }
- inBad.close();
- }
- }
- #else
- numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames);
-
- if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
-
- #endif
-
+ if(processors == 1){ numFastaSeqs = driver(lines[0], goodSeqFile, badAccnosFile, fastafile, badSeqNames); }
+ else{ numFastaSeqs = createProcesses(goodSeqFile, badAccnosFile, fastafile, badSeqNames); }
+
+ if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
#endif
#ifdef USE_MPI
if(namefile != "" && groupfile != "") {
screenNameGroupFile(badSeqNames);
- if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
}else if(namefile != "") {
screenNameGroupFile(badSeqNames);
- if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
}else if(groupfile != "") { screenGroupFile(badSeqNames); } // this screens just the group
- if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
if(alignreport != "") { screenAlignReport(badSeqNames); }
+ if(qualfile != "") { screenQual(badSeqNames); }
+ if(taxonomy != "") { screenTaxonomy(badSeqNames); }
- if (m->control_pressed) { remove(goodSeqFile.c_str()); return 0; }
+ if (m->control_pressed) { m->mothurRemove(goodSeqFile); return 0; }
#ifdef USE_MPI
}
if (itTypes != outputTypes.end()) {
if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
}
+
+ itTypes = outputTypes.find("qfile");
+ if (itTypes != outputTypes.end()) {
+ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setQualFile(current); }
+ }
+
+ itTypes = outputTypes.find("taxonomy");
+ if (itTypes != outputTypes.end()) {
+ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
+ }
m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
m->mothurOutEndLine();
ofstream goodNameOut; m->openOutputFile(goodNameFile, goodNameOut);
while(!inputNames.eof()){
- if (m->control_pressed) { goodNameOut.close(); inputNames.close(); remove(goodNameFile.c_str()); return 0; }
+ if (m->control_pressed) { goodNameOut.close(); inputNames.close(); m->mothurRemove(goodNameFile); return 0; }
inputNames >> seqName >> seqList;
it = badSeqNames.find(seqName);
ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
while(!inputGroups.eof()){
- if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodNameFile.c_str()); remove(goodGroupFile.c_str()); return 0; }
+ if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodNameFile); m->mothurRemove(goodGroupFile); return 0; }
inputGroups >> seqName >> group;
}
}
//***************************************************************************************************************
-int ScreenSeqsCommand::getSummary(vector<unsigned long int>& positions){
+int ScreenSeqsCommand::getSummary(vector<unsigned long long>& positions){
try {
vector<int> startPosition;
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)]));
- }
+ vector<unsigned long long> positions;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ positions = m->divideFile(fastafile, processors);
+ for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }
+#else
+ if(processors == 1){ lines.push_back(linePair(0, 1000)); }
+ else {
+ int numFastaSeqs = 0;
+ positions = m->setFilePosFasta(fastafile, numFastaSeqs);
+ if (positions.size() < processors) { processors = positions.size(); }
+
+ //figure out how many sequences you have to process
+ int numSeqsPerProcessor = numFastaSeqs / processors;
+ for (int i = 0; i < processors; i++) {
+ int startIndex = i * numSeqsPerProcessor;
+ if(i == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor; }
+ lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
+ }
+ }
+#endif
+#ifdef USE_MPI
+ int pid;
+ MPI_Comm_rank(MPI_COMM_WORLD, &pid);
+
+ if (pid == 0) {
+ driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
+#else
int numSeqs = 0;
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
if(processors == 1){
numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
}else{
}
if (m->control_pressed) { return 0; }
- #else
- numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
- if (m->control_pressed) { return 0; }
- #endif
-
+ //#else
+ // numSeqs = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[0]);
+ // if (m->control_pressed) { return 0; }
+ //#endif
+#endif
sort(startPosition.begin(), startPosition.end());
sort(endPosition.begin(), endPosition.end());
sort(seqLength.begin(), seqLength.end());
else if (optimize[i] == "maxlength") { maxLength = seqLength[criteriaPercentile]; m->mothurOut("Optimizing maxlength to " + toString(maxLength) + "."); m->mothurOutEndLine(); }
}
+#ifdef USE_MPI
+ }
+
+ MPI_Status status;
+ MPI_Comm_rank(MPI_COMM_WORLD, &pid);
+ MPI_Comm_size(MPI_COMM_WORLD, &processors);
+
+ if (pid == 0) {
+ //send file positions to all processes
+ for(int i = 1; i < processors; i++) {
+ MPI_Send(&startPos, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ MPI_Send(&endPos, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ MPI_Send(&maxAmbig, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ MPI_Send(&maxHomoP, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ MPI_Send(&minLength, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ MPI_Send(&maxLength, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ }
+ }else {
+ MPI_Recv(&startPos, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ MPI_Recv(&endPos, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ MPI_Recv(&maxAmbig, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ MPI_Recv(&maxHomoP, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ MPI_Recv(&minLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ MPI_Recv(&maxLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ }
+ MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
+#endif
return 0;
}
catch(exception& e) {
}
}
/**************************************************************************************/
-int ScreenSeqsCommand::driverCreateSummary(vector<int>& startPosition, vector<int>& endPosition, vector<int>& seqLength, vector<int>& ambigBases, vector<int>& longHomoPolymer, string filename, linePair* filePos) {
+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);
+ in.seekg(filePos.start);
bool done = false;
int count = 0;
count++;
}
-
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- unsigned long int pos = in.tellg();
- if ((pos == -1) || (pos >= filePos->end)) { break; }
+ //if((count) % 100 == 0){ m->mothurOut("Optimizing sequence: " + toString(count)); m->mothurOutEndLine(); }
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ unsigned long long pos = in.tellg();
+ if ((pos == -1) || (pos >= filePos.end)) { break; }
#else
if (in.eof()) { break; }
#endif
/**************************************************************************************************/
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 process = 1;
int num = 0;
- processIDS.clear();
-
+ vector<int> processIDS;
+
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
for (int k = 0; k < tempNum; k++) { in >> temp; longHomoPolymer.push_back(temp); } m->gobble(in);
in.close();
- remove(tempFilename.c_str());
+ m->mothurRemove(tempFilename);
}
- return num;
+
+#else
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+ //Windows version shared memory, so be careful when passing variables through the seqSumData struct.
+ //Above fork() will clone, so memory is separate, but that's not the case with windows,
+ //Taking advantage of shared memory to allow both threads to add info to vectors.
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ vector<sumData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+
+ //Create processor worker threads.
+ for( int i=0; i<processors-1; i++ ){
+
+ // Allocate memory for thread data.
+ sumData* tempSum = new sumData(filename, m, lines[i].start, lines[i].end, namefile, nameMap);
+ pDataArray.push_back(tempSum);
+
+ //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
+ //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
+ hThreadArray[i] = CreateThread(NULL, 0, MySumThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
+ }
+
+ //do your part
+ num = driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, fastafile, lines[processors-1]);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ num += pDataArray[i]->count;
+ for (int k = 0; k < pDataArray[i]->startPosition.size(); k++) { startPosition.push_back(pDataArray[i]->startPosition[k]); }
+ for (int k = 0; k < pDataArray[i]->endPosition.size(); k++) { endPosition.push_back(pDataArray[i]->endPosition[k]); }
+ for (int k = 0; k < pDataArray[i]->seqLength.size(); k++) { seqLength.push_back(pDataArray[i]->seqLength[k]); }
+ for (int k = 0; k < pDataArray[i]->ambigBases.size(); k++) { ambigBases.push_back(pDataArray[i]->ambigBases[k]); }
+ for (int k = 0; k < pDataArray[i]->longHomoPolymer.size(); k++) { longHomoPolymer.push_back(pDataArray[i]->longHomoPolymer[k]); }
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+
#endif
+ return num;
}
catch(exception& e) {
m->errorOut(e, "ScreenSeqsCommand", "createProcessesCreateSummary");
ofstream goodGroupOut; m->openOutputFile(goodGroupFile, goodGroupOut);
while(!inputGroups.eof()){
- if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); return 0; }
+ if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodGroupFile); return 0; }
inputGroups >> seqName >> group;
it = badSeqNames.find(seqName);
m->gobble(inputGroups);
}
- if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); remove(goodGroupFile.c_str()); return 0; }
+ if (m->control_pressed) { goodGroupOut.close(); inputGroups.close(); m->mothurRemove(goodGroupFile); return 0; }
//we were unable to remove some of the bad sequences
if (badSeqNames.size() != 0) {
inputGroups.close();
goodGroupOut.close();
- if (m->control_pressed) { remove(goodGroupFile.c_str()); }
+ if (m->control_pressed) { m->mothurRemove(goodGroupFile); }
return 0;
}
while(!inputAlignReport.eof()){
- if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); return 0; }
+ if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); m->mothurRemove(goodAlignReportFile); return 0; }
inputAlignReport >> seqName;
it = badSeqNames.find(seqName);
m->gobble(inputAlignReport);
}
- if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); remove(goodAlignReportFile.c_str()); return 0; }
+ if (m->control_pressed) { goodAlignReportOut.close(); inputAlignReport.close(); m->mothurRemove(goodAlignReportFile); return 0; }
//we were unable to remove some of the bad sequences
if (badSeqNames.size() != 0) {
inputAlignReport.close();
goodAlignReportOut.close();
- if (m->control_pressed) { remove(goodAlignReportFile.c_str()); return 0; }
+ if (m->control_pressed) { m->mothurRemove(goodAlignReportFile); return 0; }
return 0;
exit(1);
}
+}
+//***************************************************************************************************************
+
+int ScreenSeqsCommand::screenTaxonomy(set<string> badSeqNames){
+ try {
+ ifstream input;
+ m->openInputFile(taxonomy, input);
+ string seqName, tax;
+ set<string>::iterator it;
+
+ string goodTaxFile = outputDir + m->getRootName(m->getSimpleName(taxonomy)) + "good" + m->getExtension(taxonomy);
+ outputNames.push_back(goodTaxFile); outputTypes["taxonomy"].push_back(goodTaxFile);
+ ofstream goodTaxOut; m->openOutputFile(goodTaxFile, goodTaxOut);
+
+ while(!input.eof()){
+ if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); return 0; }
+
+ input >> seqName >> tax;
+ it = badSeqNames.find(seqName);
+
+ if(it != badSeqNames.end()){ badSeqNames.erase(it); }
+ else{
+ goodTaxOut << seqName << '\t' << tax << endl;
+ }
+ m->gobble(input);
+ }
+
+ if (m->control_pressed) { goodTaxOut.close(); input.close(); m->mothurRemove(goodTaxFile); 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 taxonomy file does not include the sequence " + *it + " please correct.");
+ m->mothurOutEndLine();
+ }
+ }
+
+ input.close();
+ goodTaxOut.close();
+
+ if (m->control_pressed) { m->mothurRemove(goodTaxFile); return 0; }
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "screenTaxonomy");
+ exit(1);
+ }
+
+}
+//***************************************************************************************************************
+
+int ScreenSeqsCommand::screenQual(set<string> badSeqNames){
+ try {
+ ifstream in;
+ m->openInputFile(qualfile, in);
+ set<string>::iterator it;
+
+ string goodQualFile = outputDir + m->getRootName(m->getSimpleName(qualfile)) + "good" + m->getExtension(qualfile);
+ outputNames.push_back(goodQualFile); outputTypes["qfile"].push_back(goodQualFile);
+ ofstream goodQual; m->openOutputFile(goodQualFile, goodQual);
+
+ while(!in.eof()){
+
+ if (m->control_pressed) { goodQual.close(); in.close(); m->mothurRemove(goodQualFile); return 0; }
+
+ string saveName = "";
+ string name = "";
+ string scores = "";
+
+ in >> name;
+
+ if (name.length() != 0) {
+ saveName = name.substr(1);
+ while (!in.eof()) {
+ char c = in.get();
+ if (c == 10 || c == 13){ break; }
+ else { name += c; }
+ }
+ m->gobble(in);
+ }
+
+ while(in){
+ char letter= in.get();
+ if(letter == '>'){ in.putback(letter); break; }
+ else{ scores += letter; }
+ }
+
+ m->gobble(in);
+
+ it = badSeqNames.find(saveName);
+
+ if(it != badSeqNames.end()){
+ badSeqNames.erase(it);
+ }else{
+ goodQual << name << endl << scores;
+ }
+
+ m->gobble(in);
+ }
+
+ in.close();
+ goodQual.close();
+
+ //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 qual file does not include the sequence " + *it + " please correct.");
+ m->mothurOutEndLine();
+ }
+ }
+
+ if (m->control_pressed) { m->mothurRemove(goodQualFile); return 0; }
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "screenQual");
+ exit(1);
+ }
+
}
//**********************************************************************************************************************
-int ScreenSeqsCommand::driver(linePair* filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
+int ScreenSeqsCommand::driver(linePair filePos, string goodFName, string badAccnosFName, string filename, set<string>& badSeqNames){
try {
ofstream goodFile;
m->openOutputFile(goodFName, goodFile);
ifstream inFASTA;
m->openInputFile(filename, inFASTA);
- inFASTA.seekg(filePos->start);
+ inFASTA.seekg(filePos.start);
bool done = false;
int count = 0;
count++;
}
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- unsigned long int pos = inFASTA.tellg();
- if ((pos == -1) || (pos >= filePos->end)) { break; }
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ unsigned long long pos = inFASTA.tellg();
+ if ((pos == -1) || (pos >= filePos.end)) { break; }
#else
if (inFASTA.eof()) { break; }
#endif
}
//**********************************************************************************************************************
#ifdef USE_MPI
-int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badAccnosFile, vector<unsigned long int>& MPIPos, set<string>& badSeqNames){
+int ScreenSeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& goodFile, MPI_File& badAccnosFile, vector<unsigned long long>& MPIPos, set<string>& badSeqNames){
try {
string outputString = "";
MPI_Status statusGood;
int ScreenSeqsCommand::createProcesses(string goodFileName, string badAccnos, string filename, set<string>& badSeqNames) {
try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- int process = 0;
+
+ vector<int> processIDS;
+ int process = 1;
int num = 0;
-
+
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
}
}
+ num = driver(lines[0], goodFileName, badAccnos, filename, badSeqNames);
+
//force parent to wait until all the processes are done
- for (int i=0;i<processors;i++) {
+ for (int i=0;i<processIDS.size();i++) {
int temp = processIDS[i];
wait(&temp);
}
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());
+ in.close(); m->mothurRemove(tempFile);
+
+ m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
+ m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
+
+ m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
+ m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
}
- return num;
-#endif
+ //read badSeqs in because root process doesnt know what other "bad" seqs the children found
+ ifstream inBad;
+ int ableToOpen = m->openInputFile(badAccnos, inBad, "no error");
+
+ if (ableToOpen == 0) {
+ badSeqNames.clear();
+ string tempName;
+ while (!inBad.eof()) {
+ inBad >> tempName; m->gobble(inBad);
+ badSeqNames.insert(tempName);
+ }
+ inBad.close();
+ }
+#else
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+ //Windows version shared memory, so be careful when passing variables through the sumScreenData struct.
+ //Above fork() will clone, so memory is separate, but that's not the case with windows,
+ //Taking advantage of shared memory to allow both threads to add info to badSeqNames.
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ vector<sumScreenData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+
+ //Create processor worker threads.
+ for( int i=0; i<processors-1; i++ ){
+
+ string extension = "";
+ if (i!=0) {extension += toString(i) + ".temp"; processIDS.push_back(i); }
+
+ // Allocate memory for thread data.
+ sumScreenData* tempSum = new sumScreenData(startPos, endPos, maxAmbig, maxHomoP, minLength, maxLength, filename, m, lines[i].start, lines[i].end,goodFileName+extension, badAccnos+extension);
+ pDataArray.push_back(tempSum);
+
+ //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
+ hThreadArray[i] = CreateThread(NULL, 0, MySumScreenThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
+ }
+
+ //do your part
+ num = driver(lines[processors-1], (goodFileName+toString(processors-1)+".temp"), (badAccnos+toString(processors-1)+".temp"), filename, badSeqNames);
+ processIDS.push_back(processors-1);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ num += pDataArray[i]->count;
+ for (set<string>::iterator it = pDataArray[i]->badSeqNames.begin(); it != pDataArray[i]->badSeqNames.end(); it++) { badSeqNames.insert(*it); }
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+
+ for (int i = 0; i < processIDS.size(); i++) {
+ m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
+ m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
+
+ m->appendFiles((badAccnos + toString(processIDS[i]) + ".temp"), badAccnos);
+ m->mothurRemove((badAccnos + toString(processIDS[i]) + ".temp"));
+ }
+
+#endif
+
+ return num;
+
}
catch(exception& e) {
m->errorOut(e, "ScreenSeqsCommand", "createProcesses");
projects / mothur.git / blobdiff