X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=pintail.cpp;h=b9f2434d19d4cc842bfac1dce4f6de5a0d1b742d;hp=a204294f515ee17f4e8e62f0045d8c67dfe6e406;hb=df7e3ff9f68ef157b0328a2d353c3258c5d45d89;hpb=dbd5da8043df1cb9f5ff7c6ddb5f550ea49b52c2 diff --git a/pintail.cpp b/pintail.cpp index a204294..b9f2434 100644 --- a/pintail.cpp +++ b/pintail.cpp @@ -9,95 +9,72 @@ #include "pintail.h" #include "ignoregaps.h" +#include "eachgapdist.h" +//******************************************************************************************************************** +//sorts lowest to highest +inline bool compareQuanMembers(quanMember left, quanMember right){ + return (left.score < right.score); +} //*************************************************************************************************************** -Pintail::Pintail(string filename, string temp) { fastafile = filename; templateFile = temp; } -//*************************************************************************************************************** - -Pintail::~Pintail() { +Pintail::Pintail(string filename, string temp, bool f, int p, string mask, string cons, string q, int win, int inc, string o) : Chimera() { try { - for (int i = 0; i < querySeqs.size(); i++) { delete querySeqs[i]; } - for (int i = 0; i < templateSeqs.size(); i++) { delete templateSeqs[i]; } + + fastafile = filename; + templateFileName = temp; templateSeqs = readSeqs(temp); + filter = f; + processors = p; + setMask(mask); + consfile = cons; + quanfile = q; + window = win; + increment = inc; + outputDir = o; + + distcalculator = new eachGapDist(); + decalc = new DeCalculator(); + + doPrep(); } catch(exception& e) { - errorOut(e, "Pintail", "~Pintail"); + m->errorOut(e, "Pintail", "Pintail"); exit(1); } -} + +} //*************************************************************************************************************** -void Pintail::print(ostream& out) { + +Pintail::~Pintail() { try { - for (int i = 0; i < querySeqs.size(); i++) { - - int index = ceil(deviation[i]); - - //is your DE value higher than the 95% - string chimera; - if (DE[i] > quantiles[index][4]) { chimera = "Yes"; } - else { chimera = "No"; } - - out << querySeqs[i]->getName() << '\t' << "div: " << deviation[i] << "\tstDev: " << DE[i] << "\tchimera flag: " << chimera << endl; - out << "Observed\t"; - - for (int j = 0; j < obsDistance[i].size(); j++) { out << obsDistance[i][j] << '\t'; } - out << endl; - - out << "Expected\t"; - - for (int m = 0; m < expectedDistance[i].size(); m++) { out << expectedDistance[i][m] << '\t'; } - out << endl; - - } + delete distcalculator; + delete decalc; } catch(exception& e) { - errorOut(e, "Pintail", "print"); + m->errorOut(e, "Pintail", "~Pintail"); exit(1); } } - //*************************************************************************************************************** -void Pintail::getChimeras() { +int Pintail::doPrep() { try { - //read in query sequences and subject sequences - mothurOut("Reading sequences and template file... "); cout.flush(); - querySeqs = readSeqs(fastafile); - templateSeqs = readSeqs(templateFile); - mothurOut("Done."); mothurOutEndLine(); - - int numSeqs = querySeqs.size(); - - obsDistance.resize(numSeqs); - expectedDistance.resize(numSeqs); - seqCoef.resize(numSeqs); - DE.resize(numSeqs); - Qav.resize(numSeqs); - bestfit.resize(numSeqs); - deviation.resize(numSeqs); - trimmed.resize(numSeqs); - windowSizes.resize(numSeqs, window); + mergedFilterString = ""; windowSizesTemplate.resize(templateSeqs.size(), window); - windowsForeachQuery.resize(numSeqs); quantiles.resize(100); //one for every percent mismatch + quantilesMembers.resize(100); //one for every percent mismatch - //break up file if needed - int linesPerProcess = processors / numSeqs; + //if the user does not enter a mask then you want to keep all the spots in the alignment + if (seqMask.length() == 0) { decalc->setAlignmentLength(templateSeqs[0]->getAligned().length()); } + else { decalc->setAlignmentLength(seqMask.length()); } - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - //find breakup of sequences for all times we will Parallelize - if (processors == 1) { lines.push_back(new linePair(0, numSeqs)); } - else { - //fill line pairs - for (int i = 0; i < (processors-1); i++) { - lines.push_back(new linePair((i*linesPerProcess), ((i*linesPerProcess) + linesPerProcess))); - } - //this is necessary to get remainder of processors / numSeqs so you don't miss any lines at the end - int i = processors - 1; - lines.push_back(new linePair((i*linesPerProcess), numSeqs)); - } - + decalc->setMask(seqMask); + + #ifdef USE_MPI + //do nothing + #else + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) //find breakup of templatefile for quantiles if (processors == 1) { templateLines.push_back(new linePair(0, templateSeqs.size())); } else { @@ -108,651 +85,454 @@ void Pintail::getChimeras() { } } #else - lines.push_back(new linePair(0, numSeqs)); templateLines.push_back(new linePair(0, templateSeqs.size())); #endif + #endif - distcalculator = new ignoreGaps(); - - - if (processors == 1) { - mothurOut("Finding closest sequence in template to each sequence... "); cout.flush(); - bestfit = findPairs(lines[0]->start, lines[0]->end); - - //ex.align matches from wigeon -/*for (int m = 0; m < templateSeqs.size(); m++) { - if (templateSeqs[m]->getName() == "159481") { bestfit[17] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "100137") { bestfit[16] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "112956") { bestfit[15] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "102326") { bestfit[14] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "66229") { bestfit[13] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "206276") { bestfit[12] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "63607") { bestfit[11] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "7056") { bestfit[10] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "7088") { bestfit[9] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "17553") { bestfit[8] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "131723") { bestfit[7] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "69013") { bestfit[6] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "24543") { bestfit[5] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "27824") { bestfit[4] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "1456") { bestfit[3] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "1456") { bestfit[2] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "141312") { bestfit[1] = *(templateSeqs[m]); } - if (templateSeqs[m]->getName() == "141312") { bestfit[0] = *(templateSeqs[m]); } - - -}*/ - - for (int j = 0; j < bestfit.size(); j++) { - //chops off beginning and end of sequences so they both start and end with a base - trimSeqs(querySeqs[j], bestfit[j], trimmed[j]); - } - mothurOut("Done."); mothurOutEndLine(); - - mothurOut("Finding window breaks... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - it = trimmed[i].begin(); -cout << "trimmed = " << it->first << '\t' << it->second << endl; - vector win = findWindows(querySeqs[i], it->first, it->second, windowSizes[i]); - windowsForeachQuery[i] = win; - } - mothurOut("Done."); mothurOutEndLine(); - - }else { createProcessesSpots(); } - - //find P - mothurOut("Getting conservation... "); cout.flush(); + m->mothurOut("Getting conservation... "); cout.flush(); if (consfile == "") { - mothurOut("Calculating probability of conservation for your template sequences. This can take a while... I will output the quantiles to a .prob file so that you can input them using the conservation parameter next time you run this command. Providing the .prob file will dramatically improve speed. "); cout.flush(); - probabilityProfile = calcFreq(templateSeqs); - mothurOut("Done."); mothurOutEndLine(); - }else { probabilityProfile = readFreq(); } + m->mothurOut("Calculating probability of conservation for your template sequences. This can take a while... I will output the frequency of the highest base in each position to a .freq file so that you can input them using the conservation parameter next time you run this command. Providing the .freq file will improve speed. "); cout.flush(); + probabilityProfile = decalc->calcFreq(templateSeqs, templateFileName); + if (m->control_pressed) { return 0; } + m->mothurOut("Done."); m->mothurOutEndLine(); + }else { probabilityProfile = readFreq(); m->mothurOut("Done."); } + m->mothurOutEndLine(); //make P into Q - for (int i = 0; i < probabilityProfile.size(); i++) { probabilityProfile[i] = 1 - probabilityProfile[i]; } - mothurOut("Done."); mothurOutEndLine(); + for (int i = 0; i < probabilityProfile.size(); i++) { probabilityProfile[i] = 1 - probabilityProfile[i]; } // - if (processors == 1) { + bool reRead = false; + //create filter if needed for later + if (filter) { - mothurOut("Calculating observed distance... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - cout << querySeqs[i]->getName() << '\t' << bestfit[i].getName() << " windows = " << windowsForeachQuery[i].size() << " size = " << windowSizes[i] << endl; - vector obsi = calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); - obsDistance[i] = obsi; - } - mothurOut("Done."); mothurOutEndLine(); - - - mothurOut("Finding variability... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector q = findQav(windowsForeachQuery[i], windowSizes[i]); - Qav[i] = q; - } - mothurOut("Done."); mothurOutEndLine(); - - - mothurOut("Calculating alpha... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - float alpha = getCoef(obsDistance[i], Qav[i]); - seqCoef.push_back(alpha); - } - mothurOut("Done."); mothurOutEndLine(); - - - mothurOut("Calculating expected distance... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector exp = calcExpected(Qav[i], seqCoef[i]); - expectedDistance[i] = exp; + //read in all query seqs + vector tempQuerySeqs = readSeqs(fastafile); + + vector temp; + //merge query seqs and template seqs + temp = templateSeqs; + for (int i = 0; i < tempQuerySeqs.size(); i++) { temp.push_back(tempQuerySeqs[i]); } + + if (seqMask != "") { + reRead = true; + //mask templates + for (int i = 0; i < temp.size(); i++) { + if (m->control_pressed) { + for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i]; } + return 0; + } + decalc->runMask(temp[i]); + } } - mothurOut("Done."); mothurOutEndLine(); - + + mergedFilterString = createFilter(temp, 0.5); - mothurOut("Finding deviation... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - float de = calcDE(obsDistance[i], expectedDistance[i]); - DE[i] = de; - - it = trimmed[i].begin(); - float dist = calcDist(querySeqs[i], bestfit[i], it->first, it->second); - deviation[i] = dist; + if (m->control_pressed) { + for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i]; } + return 0; } - mothurOut("Done."); mothurOutEndLine(); - } - else { createProcesses(); } + //reread template seqs + for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i]; } + } //quantiles are used to determine whether the de values found indicate a chimera //if you have to calculate them, its time intensive because you are finding the de and deviation values for each //combination of sequences in the template - if (quanfile != "") { quantiles = readQuantiles(); } - else { + if (quanfile != "") { + quantiles = readQuantiles(); + }else { + if ((!filter) && (seqMask != "")) { //if you didn't filter but you want to mask. if you filtered then you did mask first above. + reRead = true; + //mask templates + for (int i = 0; i < templateSeqs.size(); i++) { + if (m->control_pressed) { return 0; } + decalc->runMask(templateSeqs[i]); + } + } - mothurOut("Calculating quantiles for your template. This can take a while... I will output the quantiles to a .quan file that you can input them using the quantiles parameter next time you run this command. Providing the .quan file will dramatically improve speed. "); cout.flush(); + if (filter) { + reRead = true; + for (int i = 0; i < templateSeqs.size(); i++) { + if (m->control_pressed) { return 0; } + runFilter(templateSeqs[i]); + } + } + + m->mothurOut("Calculating quantiles for your template. This can take a while... I will output the quantiles to a .quan file that you can input them using the quantiles parameter next time you run this command. Providing the .quan file will dramatically improve speed. "); cout.flush(); if (processors == 1) { - quantiles = getQuantiles(0, templateSeqs.size()); + quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size()); }else { createProcessesQuan(); } + + if (m->control_pressed) { return 0; } - ofstream out4; - string o = getRootName(templateFile) + "quan"; + string noOutliers, outliers; - openOutputFile(o, out4); + if ((!filter) && (seqMask == "")) { + noOutliers = m->getRootName(m->getSimpleName(templateFileName)) + "pintail.quan"; + }else if ((!filter) && (seqMask != "")) { + noOutliers =m->getRootName(m->getSimpleName(templateFileName)) + "pintail.masked.quan"; + }else if ((filter) && (seqMask != "")) { + noOutliers = m->getRootName(m->getSimpleName(templateFileName)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastafile)) + "masked.quan"; + }else if ((filter) && (seqMask == "")) { + noOutliers = m->getRootName(m->getSimpleName(templateFileName)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastafile)) + "quan"; + } + + decalc->removeObviousOutliers(quantilesMembers, templateSeqs.size()); + + if (m->control_pressed) { return 0; } + + string outputString = "#" + m->getVersion() + "\n"; //adjust quantiles - for (int i = 0; i < quantiles.size(); i++) { - if (quantiles[i].size() == 0) { + for (int i = 0; i < quantilesMembers.size(); i++) { + vector temp; + + if (quantilesMembers[i].size() == 0) { //in case this is not a distance found in your template files for (int g = 0; g < 6; g++) { - quantiles[i].push_back(0.0); + temp.push_back(0.0); } }else{ - sort(quantiles[i].begin(), quantiles[i].end()); + sort(quantilesMembers[i].begin(), quantilesMembers[i].end()); - vector temp; //save 10% - temp.push_back(quantiles[i][int(quantiles[i].size() * 0.10)]); + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.10)]); //save 25% - temp.push_back(quantiles[i][int(quantiles[i].size() * 0.25)]); + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.25)]); //save 50% - temp.push_back(quantiles[i][int(quantiles[i].size() * 0.5)]); + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.5)]); //save 75% - temp.push_back(quantiles[i][int(quantiles[i].size() * 0.75)]); + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.75)]); //save 95% - temp.push_back(quantiles[i][int(quantiles[i].size() * 0.95)]); + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.95)]); //save 99% - temp.push_back(quantiles[i][int(quantiles[i].size() * 0.99)]); + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.99)]); - quantiles[i] = temp; } //output quan value - out4 << i+1 << '\t'; - for (int u = 0; u < quantiles[i].size(); u++) { out4 << quantiles[i][u] << '\t'; } - out4 << endl; - + outputString += toString(i+1) + "\t"; + for (int u = 0; u < temp.size(); u++) { outputString += toString(temp[u]) + "\t"; } + outputString += "\n"; + + quantiles[i] = temp; + } - mothurOut("Done."); mothurOutEndLine(); - } - - //free memory - for (int i = 0; i < lines.size(); i++) { delete lines[i]; } - for (int i = 0; i < templateLines.size(); i++) { delete templateLines[i]; } + printQuanFile(noOutliers, outputString); - delete distcalculator; - } - catch(exception& e) { - errorOut(e, "Pintail", "getChimeras"); - exit(1); - } -} -//*************************************************************************************************************** -//num is query's spot in querySeqs -void Pintail::trimSeqs(Sequence* query, Sequence subject, map& trim) { - try { - - string q = query->getAligned(); - string s = subject.getAligned(); - - int front = 0; - for (int i = 0; i < q.length(); i++) { - if (isalpha(q[i]) && isalpha(s[i])) { front = i; break; } + //free memory + quantilesMembers.clear(); + + m->mothurOut("Done."); m->mothurOutEndLine(); } - int back = 0; - for (int i = q.length(); i >= 0; i--) { - if (isalpha(q[i]) && isalpha(s[i])) { back = i; break; } + if (reRead) { + for (int i = 0; i < templateSeqs.size(); i++) { delete templateSeqs[i]; } + templateSeqs.clear(); + templateSeqs = readSeqs(templateFileName); } + - trim[front] = back; + //free memory + for (int i = 0; i < templateLines.size(); i++) { delete templateLines[i]; } + + return 0; } catch(exception& e) { - errorOut(e, "Pintail", "trimSeqs"); + m->errorOut(e, "Pintail", "doPrep"); exit(1); } } - //*************************************************************************************************************** - -vector Pintail::readFreq() { +Sequence Pintail::print(ostream& out, ostream& outAcc) { try { - - ifstream in; - openInputFile(consfile, in); - vector prob; + int index = ceil(deviation); - //read in probabilities and store in vector - int pos; float num; + //is your DE value higher than the 95% + string chimera; + if (index != 0) { //if index is 0 then its an exact match to a template seq + if (quantiles[index][4] == 0.0) { + chimera = "Your template does not include sequences that provide quantile values at distance " + toString(index); + }else { + if (DE > quantiles[index][4]) { chimera = "Yes"; } + else { chimera = "No"; } + } + }else{ chimera = "No"; } - while(!in.eof()){ - - in >> pos >> num; - - //do you want this spot - prob.push_back(num); - - gobble(in); + out << querySeq->getName() << '\t' << "div: " << deviation << "\tstDev: " << DE << "\tchimera flag: " << chimera << endl; + if (chimera == "Yes") { + m->mothurOut(querySeq->getName() + "\tdiv: " + toString(deviation) + "\tstDev: " + toString(DE) + "\tchimera flag: " + chimera); m->mothurOutEndLine(); + outAcc << querySeq->getName() << endl; } + out << "Observed\t"; - in.close(); - return prob; - - } - catch(exception& e) { - errorOut(e, "Pintail", "readFreq"); - exit(1); - } -} - -//*************************************************************************************************************** - -vector< vector > Pintail::readQuantiles() { - try { - - ifstream in; - openInputFile(quanfile, in); + for (int j = 0; j < obsDistance.size(); j++) { out << obsDistance[j] << '\t'; } + out << endl; - vector< vector > quan; - - int num; float ten, twentyfive, fifty, seventyfive, ninetyfive, ninetynine; + out << "Expected\t"; - while(!in.eof()){ - - in >> num >> ten >> twentyfive >> fifty >> seventyfive >> ninetyfive >> ninetynine; - - vector temp; - - temp.push_back(ten); - temp.push_back(twentyfive); - temp.push_back(fifty); - temp.push_back(seventyfive); - temp.push_back(ninetyfive); - temp.push_back(ninetynine); - - quan.push_back(temp); - - gobble(in); - } + for (int m = 0; m < expectedDistance.size(); m++) { out << expectedDistance[m] << '\t'; } + out << endl; - in.close(); - return quan; + return *querySeq; } catch(exception& e) { - errorOut(e, "Pintail", "readQuantiles"); + m->errorOut(e, "Pintail", "print"); exit(1); } } +#ifdef USE_MPI //*************************************************************************************************************** -//calculate the distances from each query sequence to all sequences in the template to find the closest sequence -vector Pintail::findPairs(int start, int end) { +Sequence Pintail::print(MPI_File& out, MPI_File& outAcc) { try { - vector seqsMatches; seqsMatches.resize(end-start); + string outputString = ""; + int index = ceil(deviation); + + //is your DE value higher than the 95% + string chimera; + if (index != 0) { //if index is 0 then its an exact match to a template seq + if (quantiles[index][4] == 0.0) { + chimera = "Your template does not include sequences that provide quantile values at distance " + toString(index); + }else { + if (DE > quantiles[index][4]) { chimera = "Yes"; } + else { chimera = "No"; } + } + }else{ chimera = "No"; } + + outputString += querySeq->getName() + "\tdiv: " + toString(deviation) + "\tstDev: " + toString(DE) + "\tchimera flag: " + chimera + "\n"; + if (chimera == "Yes") { + cout << querySeq->getName() << "\tdiv: " << toString(deviation) << "\tstDev: " << toString(DE) << "\tchimera flag: " << chimera << endl; + string outAccString = querySeq->getName() + "\n"; + + MPI_Status statusAcc; + int length = outAccString.length(); + char* buf = new char[length]; + memcpy(buf, outAccString.c_str(), length); + + MPI_File_write_shared(outAcc, buf, length, MPI_CHAR, &statusAcc); + delete buf; + + return *querySeq; + } + outputString += "Observed\t"; - for(int i = start; i < end; i++){ + for (int j = 0; j < obsDistance.size(); j++) { outputString += toString(obsDistance[j]) + "\t"; } + outputString += "\n"; - float smallest = 10000.0; - Sequence query = *(querySeqs[i]); + outputString += "Expected\t"; - for(int j = 0; j < templateSeqs.size(); j++){ - - Sequence temp = *(templateSeqs[j]); - - distcalculator->calcDist(query, temp); - float dist = distcalculator->getDist(); + for (int m = 0; m < expectedDistance.size(); m++) { outputString += toString(expectedDistance[m]) + "\t"; } + outputString += "\n"; + + MPI_Status status; + int length = outputString.length(); + char* buf2 = new char[length]; + memcpy(buf2, outputString.c_str(), length); - if (dist < smallest) { - seqsMatches[i] = *(templateSeqs[j]); - smallest = dist; - } - } - } + MPI_File_write_shared(out, buf2, length, MPI_CHAR, &status); + delete buf2; - return seqsMatches; - + return *querySeq; } catch(exception& e) { - errorOut(e, "Pintail", "findPairs"); + m->errorOut(e, "Pintail", "print"); exit(1); } } - +#endif //*************************************************************************************************************** -//find the window breaks for each sequence - this is so you can move ahead by bases. -vector Pintail::findWindows(Sequence* query, int front, int back, int& size) { +int Pintail::getChimeras(Sequence* query) { try { - - vector win; - - int cutoff = back - front; //back - front - - //if window is set to default - if (size == 0) { if (cutoff > 1200) { size = 300; } - else{ size = (cutoff / 4); } } - else if (size > (cutoff / 4)) { - mothurOut("You have selected to large a window size for sequence " + query->getName() + ". I will choose an appropriate window size."); mothurOutEndLine(); - size = (cutoff / 4); + querySeq = query; + trimmed.clear(); + windowSizes = window; + + //find pairs has to be done before a mask + bestfit = findPairs(query); + + if (m->control_pressed) { return 0; } + + //if they mask + if (seqMask != "") { + decalc->runMask(query); + decalc->runMask(bestfit); + } + + if (filter) { //must be done after a mask + runFilter(query); + runFilter(bestfit); } - - string seq = query->getAligned().substr(front, cutoff); - - //count bases - int numBases = 0; - for (int l = 0; l < seq.length(); l++) { if (isalpha(seq[l])) { numBases++; } } - - //save start of seq - win.push_back(front); - //move ahead increment bases at a time until all bases are in a window - int countBases = 0; - int totalBases = 0; //used to eliminate window of blanks at end of sequence - - seq = query->getAligned(); - for (int m = front; m < (back - size) ; m++) { - - //count number of bases you see - if (isalpha(seq[m])) { countBases++; totalBases++; } - - //if you have seen enough bases to make a new window - if (countBases >= increment) { - win.push_back(m); //save spot in alignment - countBases = 0; //reset bases you've seen in this window - } - //no need to continue if all your bases are in a window - if (totalBases == numBases) { break; } - } - - return win; - - } - catch(exception& e) { - errorOut(e, "Pintail", "findWindows"); - exit(1); - } -} + //trim seq + decalc->trimSeqs(query, bestfit, trimmed); + + //find windows + it = trimmed.begin(); + windowsForeachQuery = decalc->findWindows(query, it->first, it->second, windowSizes, increment); -//*************************************************************************************************************** -vector Pintail::calcObserved(Sequence* query, Sequence subject, vector window, int size) { - try { + //find observed distance + obsDistance = decalc->calcObserved(query, bestfit, windowsForeachQuery, windowSizes); - vector temp; -//cout << "query length = " << query->getAligned().length() << '\t' << " subject length = " << subject.getAligned().length() << endl; - for (int m = 0; m < window.size(); m++) { - - string seqFrag = query->getAligned().substr(window[m], size); - string seqFragsub = subject.getAligned().substr(window[m], size); - //cout << "start point = " << window[m] << " end point = " << window[m]+size << endl; - int diff = 0; - for (int b = 0; b < seqFrag.length(); b++) { - if (seqFrag[b] != seqFragsub[b]) { diff++; } - } - - //percentage of mismatched bases - float dist; - dist = diff / (float) seqFrag.length() * 100; + if (m->control_pressed) { return 0; } - temp.push_back(dist); - } - - return temp; - } - catch(exception& e) { - errorOut(e, "Pintail", "calcObserved"); - exit(1); - } -} -//*************************************************************************************************************** -float Pintail::calcDist(Sequence* query, Sequence subject, int front, int back) { - try { + Qav = decalc->findQav(windowsForeachQuery, windowSizes, probabilityProfile); - //so you only look at the trimmed part of the sequence - int cutoff = back - front; - - //from first startpoint with length back-front - string seqFrag = query->getAligned().substr(front, cutoff); - string seqFragsub = subject.getAligned().substr(front, cutoff); - - int diff = 0; - for (int b = 0; b < seqFrag.length(); b++) { - if (seqFrag[b] != seqFragsub[b]) { diff++; } - } - - //percentage of mismatched bases - float dist = diff / (float) seqFrag.length() * 100; - - return dist; - } - catch(exception& e) { - errorOut(e, "Pintail", "calcDist"); - exit(1); - } -} + if (m->control_pressed) { return 0; } -//*************************************************************************************************************** -vector Pintail::calcExpected(vector qav, float coef) { - try { + //find alpha + seqCoef = decalc->getCoef(obsDistance, Qav); - //for each window - vector queryExpected; - - for (int m = 0; m < qav.size(); m++) { - - float expected = qav[m] * coef; - - queryExpected.push_back(expected); - } - - return queryExpected; - - } - catch(exception& e) { - errorOut(e, "Pintail", "calcExpected"); - exit(1); - } -} -//*************************************************************************************************************** -float Pintail::calcDE(vector obs, vector exp) { - try { + //calculating expected distance + expectedDistance = decalc->calcExpected(Qav, seqCoef); - //for each window - float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2 - for (int m = 0; m < obs.size(); m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } - - float de = sqrt((sum / (obs.size() - 1))); - - return de; + if (m->control_pressed) { return 0; } + + //finding de + DE = decalc->calcDE(obsDistance, expectedDistance); + + if (m->control_pressed) { return 0; } + + //find distance between query and closest match + it = trimmed.begin(); + deviation = decalc->calcDist(query, bestfit, it->first, it->second); + + delete bestfit; + + return 0; } catch(exception& e) { - errorOut(e, "Pintail", "calcDE"); + m->errorOut(e, "Pintail", "getChimeras"); exit(1); } } //*************************************************************************************************************** -vector Pintail::calcFreq(vector seqs) { +vector Pintail::readFreq() { try { + //read in probabilities and store in vector + int pos; float num; vector prob; - string freqfile = getRootName(templateFile) + "prob"; - ofstream outFreq; + set h = decalc->getPos(); //positions of bases in masking sequence - openOutputFile(freqfile, outFreq); + #ifdef USE_MPI - //at each position in the sequence - for (int i = 0; i < seqs[0]->getAligned().length(); i++) { - - vector freq; freq.resize(4,0); - int gaps = 0; + MPI_File inMPI; + MPI_Offset size; + MPI_Status status; + + //char* inFileName = new char[consfile.length()]; + //memcpy(inFileName, consfile.c_str(), consfile.length()); + + char inFileName[1024]; + strcpy(inFileName, consfile.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; + delete buffer; + + if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size); } + istringstream iss (tempBuf,istringstream::in); + + //read version + string line = m->getline(iss); m->gobble(iss); + + while(!iss.eof()) { + iss >> pos >> num; + + if (h.count(pos) > 0) { + float Pi; + Pi = (num - 0.25) / 0.75; - //find the frequency of each nucleotide - for (int j = 0; j < seqs.size(); j++) { - - char value = seqs[j]->getAligned()[i]; - - if(toupper(value) == 'A') { freq[0]++; } - else if(toupper(value) == 'T' || toupper(value) == 'U') { freq[1]++; } - else if(toupper(value) == 'G') { freq[2]++; } - else if(toupper(value) == 'C') { freq[3]++; } - else { gaps++; } + //cannot have probability less than 0. + if (Pi < 0) { Pi = 0.0; } + + //do you want this spot + prob.push_back(Pi); } - //find base with highest frequency - int highest = 0; - for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } + m->gobble(iss); + } + + MPI_File_close(&inMPI); + + #else + + ifstream in; + m->openInputFile(consfile, in); + + //read version + string line = m->getline(in); m->gobble(in); + + while(!in.eof()){ - float highFreq; - //subtract gaps to "ignore them" - if ( (seqs.size() - gaps) == 0 ) { highFreq = 1.0; } - else { highFreq = highest / (float) (seqs.size() - gaps); } - - float Pi; - Pi = (highFreq - 0.25) / 0.75; + in >> pos >> num; - //cannot have probability less than 0. - if (Pi < 0) { Pi = 0.0; } + if (h.count(pos) > 0) { + float Pi; + Pi = (num - 0.25) / 0.75; - //saves this for later - outFreq << i+1 << '\t' << Pi << endl; + //cannot have probability less than 0. + if (Pi < 0) { Pi = 0.0; } + + //do you want this spot + prob.push_back(Pi); + } - prob.push_back(Pi); + m->gobble(in); } + in.close(); - outFreq.close(); - - return prob; - - } - catch(exception& e) { - errorOut(e, "Pintail", "calcFreq"); - exit(1); - } -} -//*************************************************************************************************************** -vector Pintail::findQav(vector window, int size) { - try { - vector averages; - - //for each window find average - for (int m = 0; m < window.size(); m++) { - - float average = 0.0; - - //while you are in the window for this sequence - for (int j = window[m]; j < (window[m]+size); j++) { average += probabilityProfile[j]; } - - average = average / size; + #endif - //save this windows average - averages.push_back(average); - } - - return averages; - } - catch(exception& e) { - errorOut(e, "Pintail", "findQav"); - exit(1); - } -} - -//*************************************************************************************************************** -vector< vector > Pintail::getQuantiles(int start, int end) { - try { - vector< vector > quan; - - //percentage of mismatched pairs 1 to 100 - quan.resize(100); - - - //for each sequence - for(int i = start; i < end; i++){ + return prob; - mothurOut("Processing template sequence " + toString(i)); mothurOutEndLine(); - Sequence* query = templateSeqs[i]; - - //compare to every other sequence in template - for(int j = 0; j < i; j++){ - - Sequence subject = *(templateSeqs[j]); - - map trim; - trimSeqs(query, subject, trim); - - it = trim.begin(); - int front = it->first; int back = it->second; - - //reset window for each new comparison - windowSizesTemplate[i] = window; - - vector win = findWindows(query, front, back, windowSizesTemplate[i]); - - vector obsi = calcObserved(query, subject, win, windowSizesTemplate[i]); - - vector q = findQav(win, windowSizesTemplate[i]); - - float alpha = getCoef(obsi, q); - - vector exp = calcExpected(q, alpha); - - float de = calcDE(obsi, exp); - - float dist = calcDist(query, subject, front, back); - - dist = ceil(dist); - - //dist-1 because vector indexes start at 0. - quan[dist-1].push_back(de); - - } - } - - return quan; - } catch(exception& e) { - errorOut(e, "Pintail", "findQav"); + m->errorOut(e, "Pintail", "readFreq"); exit(1); } } //*************************************************************************************************************** -float Pintail::getCoef(vector obs, vector qav) { +//calculate the distances from each query sequence to all sequences in the template to find the closest sequence +Sequence* Pintail::findPairs(Sequence* q) { try { - - //find average prob for this seqs windows - float probAverage = 0.0; - for (int j = 0; j < qav.size(); j++) { probAverage += qav[j]; } - probAverage = probAverage / (float) qav.size(); - //find observed average - float obsAverage = 0.0; - for (int j = 0; j < obs.size(); j++) { obsAverage += obs[j]; } - obsAverage = obsAverage / (float) obs.size(); + Sequence* seqsMatches; - - float coef = obsAverage / probAverage; - - return coef; + seqsMatches = decalc->findClosest(q, templateSeqs); + return seqsMatches; + } catch(exception& e) { - errorOut(e, "Pintail", "getCoef"); + m->errorOut(e, "Pintail", "findPairs"); exit(1); } } -/**************************************************************************************************/ - -void Pintail::createProcessesSpots() { +//************************************************************************************************** +void Pintail::createProcessesQuan() { try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + int process = 1; vector processIDS; - vector< vector > win; win.resize(querySeqs.size()); - + //loop through and create all the processes you want while (process != processors) { int pid = fork(); @@ -762,233 +542,234 @@ void Pintail::createProcessesSpots() { process++; }else if (pid == 0){ - mothurOut("Finding pairs for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); - vector tempbest; - tempbest = findPairs(lines[process]->start, lines[process]->end); - int count = 0; - for (int i = lines[process]->start; i < lines[process]->end; i++) { - bestfit[i] = tempbest[count]; - - //chops off beginning and end of sequences so they both start and end with a base - trimSeqs(querySeqs[i], bestfit[i], trimmed[i]); - count++; - } - mothurOut("Done finding pairs for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); + quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, templateLines[process]->start, templateLines[process]->end); - mothurOut("Finding window breaks for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); - for (int i = lines[process]->start; i < lines[process]->end; i++) { - vector temp = findWindows(querySeqs[i], it->first, it->second, windowSizes[i]); - win[i] = temp; + //write out data to file so parent can read it + ofstream out; + string s = toString(getpid()) + ".temp"; + m->openOutputFile(s, out); + + //output observed distances + for (int i = 0; i < quantilesMembers.size(); i++) { + out << quantilesMembers[i].size() << '\t'; + for (int j = 0; j < quantilesMembers[i].size(); j++) { + out << quantilesMembers[i][j] << '\t'; + } + out << endl; } - mothurOut("Done finding window breaks for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); + out.close(); + + exit(0); + }else { + m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); + for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } exit(0); - }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } + } } + //parent does its part + quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, templateLines[0]->start, templateLines[0]->end); + //force parent to wait until all the processes are done - for (int i=0;istart, lines[0]->end); - for (int j = 0; j < bestfit.size(); j++) { - //chops off beginning and end of sequences so they both start and end with a base - trimSeqs(querySeqs[j], bestfit[j], j); - } - for (int i = lines[0]->start; i < lines[0]->end; i++) { - it = trimmed[i].begin(); - map win = findWindows(querySeqs[i], it->first, it->second, windowSizes[i]); - windows[i] = win; + //get data created by processes + for (int i=0;i<(processors-1);i++) { + ifstream in; + string s = toString(processIDS[i]) + ".temp"; + m->openInputFile(s, in); + + vector< vector > quan; + quan.resize(100); + + //get quantiles + for (int h = 0; h < quan.size(); h++) { + int num; + in >> num; + + m->gobble(in); + + vector q; float w; + for (int j = 0; j < num; j++) { + in >> w; + q.push_back(w); + } + + quan[h] = q; + m->gobble(in); + } + + + //save quan in quantiles + for (int j = 0; j < quan.size(); j++) { + //put all values of q[i] into quan[i] + for (int l = 0; l < quan[j].size(); l++) { quantilesMembers[j].push_back(quan[j][l]); } + //quantilesMembers[j].insert(quantilesMembers[j].begin(), quan[j].begin(), quan[j].end()); + } + + in.close(); + m->mothurRemove(s); } +#else + quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size()); #endif } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesSpots"); + m->errorOut(e, "Pintail", "createProcessesQuan"); exit(1); } } - - -/**************************************************************************************************/ - -void Pintail::createProcesses() { +//*************************************************************************************************************** +vector< vector > Pintail::readQuantiles() { try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector processIDS; + int num; + float ten, twentyfive, fifty, seventyfive, ninetyfive, ninetynine; + + vector< vector > quan; + vector temp; temp.resize(6, 0); - vector< vector > exp; exp.resize(querySeqs.size()); - vector de; de.resize(querySeqs.size()); - vector< vector > obs; obs.resize(querySeqs.size()); - vector dev; dev.resize(querySeqs.size()); + //to fill 0 + quan.push_back(temp); + + #ifdef USE_MPI + MPI_File inMPI; + MPI_Offset size; + MPI_Status status; - //loop through and create all the processes you want - while (process != processors) { - int pid = fork(); - - if (pid > 0) { - processIDS.push_back(pid); - process++; - }else if (pid == 0){ - - mothurOut("Calculating observed, expected and de values for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); - for (int i = lines[process]->start; i < lines[process]->end; i++) { - - vector obsi = calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); - obs[i] = obsi; - - //calc Qav - vector q = findQav(windowsForeachQuery[i], windowSizes[i]); - - //get alpha - float alpha = getCoef(obsDistance[i], q); - - //find expected - vector exp = calcExpected(q, alpha); - expectedDistance[i] = exp; - - //get de and deviation - float dei = calcDE(obsi, exp); - de[i] = dei; - - it = trimmed[i].begin(); - float dist = calcDist(querySeqs[i], bestfit[i], it->first, it->second); - dev[i] = dist; - } - mothurOut("Done calculating observed, expected and de values for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); - - exit(0); - }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } - } + //char* inFileName = new char[quanfile.length()]; + //memcpy(inFileName, quanfile.c_str(), quanfile.length()); - //force parent to wait until all the processes are done - for (int i=0;i size) { tempBuf = tempBuf.substr(0, size); } + istringstream iss (tempBuf,istringstream::in); + delete buffer; - obsDistance = obs; - expectedDistance = exp; - DE = de; - deviation = dev; + //read version + string line = m->getline(iss); m->gobble(iss); -#else - mothurOut("Calculating observed distance... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector obsi = calcObserved(querySeqs[i], bestfit[i], windows[i], windowSizes[i]); - obsDistance[i] = obsi; - } - mothurOut("Done."); mothurOutEndLine(); + while(!iss.eof()) { + iss >> num >> ten >> twentyfive >> fifty >> seventyfive >> ninetyfive >> ninetynine; + temp.clear(); + temp.push_back(ten); + temp.push_back(twentyfive); + temp.push_back(fifty); + temp.push_back(seventyfive); + temp.push_back(ninetyfive); + temp.push_back(ninetynine); - mothurOut("Finding variability... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector q = findQav(windows[i], windowSizes[i]); - Qav[i] = q; - } - mothurOut("Done."); mothurOutEndLine(); - - + quan.push_back(temp); - mothurOut("Calculating alpha... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - float alpha = getCoef(obsDistance[i], Qav[i]); - seqCoef.push_back(alpha); - } - mothurOut("Done."); mothurOutEndLine(); - + m->gobble(iss); + } + + MPI_File_close(&inMPI); + #else + + ifstream in; + m->openInputFile(quanfile, in); - mothurOut("Calculating expected distance... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector exp = calcExpected(Qav[i], seqCoef[i]); - expectedDistance[i] = exp; - } - mothurOut("Done."); mothurOutEndLine(); + //read version + string line = m->getline(in); m->gobble(in); + while(!in.eof()){ + in >> num >> ten >> twentyfive >> fifty >> seventyfive >> ninetyfive >> ninetynine; - mothurOut("Finding deviation... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - float de = calcDE(obsDistance[i], expectedDistance[i]); - DE[i] = de; - - it = trimmed[i].begin(); - float dist = calcDist(querySeqs[i], bestfit[i], it->first, it->second); - deviation[i] = dist; - } - mothurOut("Done."); mothurOutEndLine(); - -#endif + temp.clear(); + + temp.push_back(ten); + temp.push_back(twentyfive); + temp.push_back(fifty); + temp.push_back(seventyfive); + temp.push_back(ninetyfive); + temp.push_back(ninetynine); + + quan.push_back(temp); + + m->gobble(in); + } + in.close(); + #endif + + return quan; + } catch(exception& e) { - errorOut(e, "Pintail", "createProcesses"); + m->errorOut(e, "Pintail", "readQuantiles"); exit(1); } } +//***************************************************************************************************************/ +void Pintail::printQuanFile(string file, string outputString) { + try { + + #ifdef USE_MPI + + MPI_File outQuan; + MPI_Status status; + + int pid; + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are -/**************************************************************************************************/ + int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; -void Pintail::createProcessesQuan() { - try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector processIDS; - vector< vector > quan; quan.resize(100); - - //loop through and create all the processes you want - while (process != processors) { - int pid = fork(); + //char* FileName = new char[file.length()]; + //memcpy(FileName, file.c_str(), file.length()); - if (pid > 0) { - processIDS.push_back(pid); - process++; - }else if (pid == 0){ - - vector< vector > q = getQuantiles(templateLines[process]->start, templateLines[process]->end); - - for (int i = 0; i < q.size(); i++) { - //put all values of q[i] into quan[i] - quan[i].insert(quan[i].begin(), q[i].begin(), q[i].end()); - } + char FileName[1024]; + strcpy(FileName, file.c_str()); + + if (pid == 0) { + MPI_File_open(MPI_COMM_SELF, FileName, outMode, MPI_INFO_NULL, &outQuan); //comm, filename, mode, info, filepointer - for (int i = 0; i < quan.size(); i++) { - cout << i+1 << '\t'; - for (int j = 0; j < quan[i].size(); j++) { cout << quan[i][j] << '\t'; } - cout << endl; - } + int length = outputString.length(); + char* buf = new char[length]; + memcpy(buf, outputString.c_str(), length); + + MPI_File_write(outQuan, buf, length, MPI_CHAR, &status); + delete buf; - exit(0); - }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } - } - - //force parent to wait until all the processes are done - for (int i=0;iopenOutputFile(file, outQuan); + + outQuan << outputString; + + outQuan.close(); + #endif } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesQuan"); + m->errorOut(e, "Pintail", "printQuanFile"); exit(1); } } +//***************************************************************************************************************/ -//***************************************************************************************************************