X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=chimeracheckrdp.cpp;h=42b531282c5f020664277a9e907699a3bf227658;hp=3e6dacdee85dbb75fb9b6535f259eead8a4f2509;hb=1a20e24ee786195ab0e1cccd4f5aede7a88f3f4e;hpb=74844a60d80c6dd06e3fb02ee9b928424f9019b0 diff --git a/chimeracheckrdp.cpp b/chimeracheckrdp.cpp index 3e6dacd..42b5312 100644 --- a/chimeracheckrdp.cpp +++ b/chimeracheckrdp.cpp @@ -10,7 +10,30 @@ #include "chimeracheckrdp.h" //*************************************************************************************************************** -ChimeraCheckRDP::ChimeraCheckRDP(string filename, string o) { fastafile = filename; outputDir = o; } +ChimeraCheckRDP::ChimeraCheckRDP(string filename, string temp, string n, bool s, int inc, int k, string o) : Chimera() { + try { + fastafile = filename; + templateFileName = temp; + name = n; + svg = s; + increment = inc; + kmerSize = k; + outputDir = o; + + templateDB = new AlignmentDB(templateFileName, "kmer", kmerSize, 0.0,0.0,0.0,0.0, rand()); + m->mothurOutEndLine(); + + kmer = new Kmer(kmerSize); + + if (name != "") { + readName(name); //fills name map with names of seqs the user wants to have .svg for. + } + } + catch(exception& e) { + m->errorOut(e, "ChimeraCheckRDP", "ChimeraCheckRDP"); + exit(1); + } +} //*************************************************************************************************************** ChimeraCheckRDP::~ChimeraCheckRDP() { @@ -19,12 +42,12 @@ ChimeraCheckRDP::~ChimeraCheckRDP() { delete kmer; } catch(exception& e) { - m->errorOut(e, "ChimeraCheckRDP", "~AlignSim"); + m->errorOut(e, "ChimeraCheckRDP", "~ChimeraCheckRDP"); exit(1); } } //*************************************************************************************************************** -void ChimeraCheckRDP::print(ostream& out, ostream& outAcc) { +Sequence ChimeraCheckRDP::print(ostream& out, ostream& outAcc) { try { m->mothurOut("Processing: " + querySeq->getName()); m->mothurOutEndLine(); @@ -48,29 +71,58 @@ void ChimeraCheckRDP::print(ostream& out, ostream& outAcc) { makeSVGpic(IS); //zeros out negative results } } + + return *querySeq; } catch(exception& e) { m->errorOut(e, "ChimeraCheckRDP", "print"); exit(1); } } +#ifdef USE_MPI //*************************************************************************************************************** -void ChimeraCheckRDP::doPrep() { +Sequence ChimeraCheckRDP::print(MPI_File& out, MPI_File& outAcc) { try { - templateDB = new AlignmentDB(templateFileName, "kmer", kmerSize, 0.0,0.0,0.0,0.0); - m->mothurOutEndLine(); - kmer = new Kmer(kmerSize); + cout << "Processing: " << querySeq->getName() << endl; - if (name != "") { - readName(name); //fills name map with names of seqs the user wants to have .svg for. + string outString = ""; + + outString += querySeq->getName() + "\nIS scores: \t"; + + for (int k = 0; k < IS.size(); k++) { + outString += toString(IS[k].score) + "\t"; + } + outString += "\n"; + + MPI_Status status; + int length = outString.length(); + char* buf = new char[length]; + memcpy(buf, outString.c_str(), length); + + MPI_File_write_shared(out, buf, length, MPI_CHAR, &status); + delete buf; + + if (svg) { + if (name != "") { //if user has specific names + map::iterator it = names.find(querySeq->getName()); + + if (it != names.end()) { //user wants pic of this + makeSVGpic(IS); //zeros out negative results + } + }else{//output them all + makeSVGpic(IS); //zeros out negative results + } } + + return *querySeq; } catch(exception& e) { - m->errorOut(e, "ChimeraCheckRDP", "doPrep"); + m->errorOut(e, "ChimeraCheckRDP", "print"); exit(1); } } +#endif //*************************************************************************************************************** int ChimeraCheckRDP::getChimeras(Sequence* query) { try { @@ -119,6 +171,8 @@ vector ChimeraCheckRDP::findIS() { //for each window for (int f = start; f < (seq.length() - start); f+=increment) { + + if (m->control_pressed) { return isValues; } if ((f - kmerSize) < 0) { m->mothurOut("Your sequence is too short for your kmerSize."); m->mothurOutEndLine(); exit(1); } @@ -200,18 +254,53 @@ vector ChimeraCheckRDP::findIS() { //*************************************************************************************************************** void ChimeraCheckRDP::readName(string namefile) { try{ - ifstream in; - openInputFile(namefile, in); + string name; + + #ifdef USE_MPI + MPI_File inMPI; + MPI_Offset size; + MPI_Status status; + + //char* inFileName = new char[namefile.length()]; + //memcpy(inFileName, namefile.c_str(), namefile.length()); + + char inFileName[1024]; + strcpy(inFileName, namefile.c_str()); + + MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI); + MPI_File_get_size(inMPI, &size); + + //delete inFileName; + + char* buffer = new char[size]; + MPI_File_read(inMPI, buffer, size, MPI_CHAR, &status); + + string tempBuf = buffer; + if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size); } + istringstream iss (tempBuf,istringstream::in); + delete buffer; + + while(!iss.eof()) { + iss >> name; m->gobble(iss); + names[name] = name; + } + + MPI_File_close(&inMPI); + + #else + + ifstream in; + m->openInputFile(namefile, in); + while (!in.eof()) { - - in >> name; - + in >> name; m->gobble(in); names[name] = name; - - gobble(in); } + in.close(); + + #endif } catch(exception& e) { @@ -226,25 +315,26 @@ int ChimeraCheckRDP::calcKmers(map query, map subject) { try{ int common = 0; - map::iterator small; - map::iterator large; + map::iterator smallone; + map::iterator largeone; + if (query.size() < subject.size()) { - for (small = query.begin(); small != query.end(); small++) { - large = subject.find(small->first); + for (smallone = query.begin(); smallone != query.end(); smallone++) { + largeone = subject.find(smallone->first); //if you found it they have that kmer in common - if (large != subject.end()) { common++; } + if (largeone != subject.end()) { common++; } } }else { - for (small = subject.begin(); small != subject.end(); small++) { - large = query.find(small->first); + for (smallone = subject.begin(); smallone != subject.end(); smallone++) { + largeone = query.find(smallone->first); //if you found it they have that kmer in common - if (large != query.end()) { common++; } + if (largeone != query.end()) { common++; } } } @@ -256,14 +346,93 @@ int ChimeraCheckRDP::calcKmers(map query, map subject) { exit(1); } } +#ifdef USE_MPI +//*************************************************************************************************************** +void ChimeraCheckRDP::makeSVGpic(vector info) { + try{ + + string file = outputDir + querySeq->getName() + ".chimeracheck.svg"; + + MPI_File outSVG; + int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; + + //char* FileName = new char[file.length()]; + //memcpy(FileName, file.c_str(), file.length()); + + char FileName[1024]; + strcpy(FileName, file.c_str()); + MPI_File_open(MPI_COMM_SELF, FileName, outMode, MPI_INFO_NULL, &outSVG); //comm, filename, mode, info, filepointer + + //delete FileName; + + int width = (info.size()*5) + 150; + + string outString = ""; + + outString += "\n"; + outString += "\n"; + outString += "Plotted IS values for " + querySeq->getName() + "\n"; + + outString += "\n"; + outString += "\n"; + + outString += "" + toString(info[0].midpoint) + "\n"; + outString += "" + toString(info[info.size()-1].midpoint) + "\n"; + outString += "Base Positions\n"; + + outString += "0\n"; + + outString += "IS\n"; + + + //find max is score + float biggest = 0.0; + for (int i = 0; i < info.size(); i++) { + if (info[i].score > biggest) { + biggest = info[i].score; + } + } + + outString += "" + toString(biggest) + "\n"; + + int scaler2 = 500 / biggest; + + + outString += " "; + for (int i = 0; i < info.size(); i++) { + if(info[i].score < 0) { info[i].score = 0; } + outString += toString(((i*5) + 75)) + "," + toString((600 - (info[i].score * scaler2))) + " "; + } + + outString += "\"/> "; + outString += "\n\n"; + + MPI_Status status; + int length = outString.length(); + char* buf2 = new char[length]; + memcpy(buf2, outString.c_str(), length); + + MPI_File_write(outSVG, buf2, length, MPI_CHAR, &status); + delete buf2; + + MPI_File_close(&outSVG); + + } + catch(exception& e) { + m->errorOut(e, "ChimeraCheckRDP", "makeSVGpic"); + exit(1); + } +} +#else //*************************************************************************************************************** void ChimeraCheckRDP::makeSVGpic(vector info) { try{ string file = outputDir + querySeq->getName() + ".chimeracheck.svg"; ofstream outsvg; - openOutputFile(file, outsvg); + m->openOutputFile(file, outsvg); int width = (info.size()*5) + 150; @@ -314,6 +483,7 @@ void ChimeraCheckRDP::makeSVGpic(vector info) { exit(1); } } -//*************************************************************************************************************** +#endif +//***************************************************************************************************************/