X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=pintail.cpp;h=7ae3ff6dbb96224197dff97db4b435880738d87b;hb=98ea55ba2d46a429f031086b4b3272780d0ec894;hp=72be67948e757f8b8d09722baf51b5a1da720c12;hpb=2ad3477dcd6a01628712b32e767353f917f2a95a;p=mothur.git diff --git a/pintail.cpp b/pintail.cpp index 72be679..7ae3ff6 100644 --- a/pintail.cpp +++ b/pintail.cpp @@ -9,7 +9,13 @@ #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; } @@ -19,6 +25,7 @@ Pintail::~Pintail() { try { for (int i = 0; i < querySeqs.size(); i++) { delete querySeqs[i]; } for (int i = 0; i < templateSeqs.size(); i++) { delete templateSeqs[i]; } + for (int i = 0; i < bestfit.size(); i++) { delete bestfit[i]; } } catch(exception& e) { errorOut(e, "Pintail", "~Pintail"); @@ -29,9 +36,25 @@ Pintail::~Pintail() { void Pintail::print(ostream& out) { try { + mothurOutEndLine(); + for (int i = 0; i < querySeqs.size(); i++) { - out << querySeqs[i]->getName() << '\t' << "div: " << deviation[i] << "\tstDev: " << DE[i] << endl; + int index = ceil(deviation[i]); + + //is your DE value higher than the 95% + string chimera; + if (quantiles[index][4] == 0.0) { + chimera = "Your template does not include sequences that provide quantile values at distance " + toString(index); + }else { + 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; + if (chimera == "Yes") { + mothurOut(querySeqs[i]->getName() + "\tdiv: " + toString(deviation[i]) + "\tstDev: " + toString(DE[i]) + "\tchimera flag: " + chimera); mothurOutEndLine(); + } out << "Observed\t"; for (int j = 0; j < obsDistance[i].size(); j++) { out << obsDistance[i][j] << '\t'; } @@ -71,10 +94,14 @@ void Pintail::getChimeras() { deviation.resize(numSeqs); trimmed.resize(numSeqs); windowSizes.resize(numSeqs, window); - windows.resize(numSeqs); + windowSizesTemplate.resize(templateSeqs.size(), window); + windowsForeachQuery.resize(numSeqs); + h.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; + int linesPerProcess = numSeqs / processors ; #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) //find breakup of sequences for all times we will Parallelize @@ -88,106 +115,147 @@ void Pintail::getChimeras() { int i = processors - 1; lines.push_back(new linePair((i*linesPerProcess), numSeqs)); } + + //find breakup of templatefile for quantiles + if (processors == 1) { templateLines.push_back(new linePair(0, templateSeqs.size())); } + else { + for (int i = 0; i < processors; i++) { + templateLines.push_back(new linePair()); + templateLines[i]->start = int (sqrt(float(i)/float(processors)) * templateSeqs.size()); + templateLines[i]->end = int (sqrt(float(i+1)/float(processors)) * templateSeqs.size()); + } + } #else lines.push_back(new linePair(0, numSeqs)); + templateLines.push_back(new linePair(0, templateSeqs.size())); #endif - distcalculator = new ignoreGaps(); - - + distcalculator = new eachGapDist(); + decalc = new DeCalculator(); + + //if the user does enter a mask then you want to keep all the spots in the alignment + if (seqMask.length() == 0) { decalc->setAlignmentLength(querySeqs[0]->getAligned().length()); } + else { decalc->setAlignmentLength(seqMask.length()); } + + decalc->setMask(seqMask); + + //find pairs if (processors == 1) { mothurOut("Finding closest sequence in template to each sequence... "); cout.flush(); bestfit = findPairs(lines[0]->start, lines[0]->end); -/*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]); } - - -}*/ + mothurOut("Done."); mothurOutEndLine(); + }else { createProcessesPairs(); } + +//string o = "closestmatch.eachgap.fasta"; +//ofstream out7; +//openOutputFile(o, out7); + +//for (int i = 0; i < bestfit.size(); i++) { + //out7 << ">" << querySeqs[i]->getName() << "-"<< bestfit[i]->getName() << endl; + //out7 << bestfit[i]->getAligned() << endl; +//} +//out7.close(); + //find P + 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 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, templateFile); + mothurOut("Done."); mothurOutEndLine(); + }else { probabilityProfile = readFreq(); } + + //make P into Q + for (int i = 0; i < probabilityProfile.size(); i++) { probabilityProfile[i] = 1 - probabilityProfile[i]; } //cout << i << '\t' << probabilityProfile[i] << endl; + mothurOut("Done."); mothurOutEndLine(); + + //mask sequences if the user wants to + if (seqMask != "") { + //mask querys + for (int i = 0; i < querySeqs.size(); i++) { + decalc->runMask(querySeqs[i]); + } + + //mask templates + for (int i = 0; i < templateSeqs.size(); i++) { + decalc->runMask(templateSeqs[i]); + } + for (int i = 0; i < bestfit.size(); i++) { + decalc->runMask(bestfit[i]); + } + + } + + if (filter) { + vector temp = templateSeqs; + for (int i = 0; i < querySeqs.size(); i++) { temp.push_back(querySeqs[i]); } + + createFilter(temp); + + runFilter(querySeqs); + runFilter(templateSeqs); + runFilter(bestfit); + } + + + if (processors == 1) { + 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); + decalc->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(); - vector win = findWindows(querySeqs[i], it->first, it->second, windowSizes[i]); - windows[i] = win; + vector win = decalc->findWindows(querySeqs[i], it->first, it->second, windowSizes[i], increment); + windowsForeachQuery[i] = win; } - - - } else { createProcessesSpots(); } - - //find P - if (consfile == "") { probabilityProfile = calcFreq(templateSeqs); } - else { probabilityProfile = readFreq(); } + mothurOut("Done."); mothurOutEndLine(); - //make P into Q - for (int i = 0; i < probabilityProfile.size(); i++) { probabilityProfile[i] = 1 - probabilityProfile[i]; } + }else { createProcessesSpots(); } if (processors == 1) { 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]); + vector obsi = decalc->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(windows[i], windowSizes[i]); + vector q = decalc->findQav(windowsForeachQuery[i], windowSizes[i], probabilityProfile); + 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); + float alpha = decalc->getCoef(obsDistance[i], Qav[i]); + seqCoef[i] = 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]); + vector exp = decalc->calcExpected(Qav[i], seqCoef[i]); expectedDistance[i] = exp; } mothurOut("Done."); mothurOutEndLine(); - mothurOut("Finding deviation... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { - float de = calcDE(obsDistance[i], expectedDistance[i]); + float de = decalc->calcDE(obsDistance[i], expectedDistance[i]); DE[i] = de; - + it = trimmed[i].begin(); - float dist = calcDist(querySeqs[i], bestfit[i], it->first, it->second); + float dist = decalc->calcDist(querySeqs[i], bestfit[i], it->first, it->second); deviation[i] = dist; } mothurOut("Done."); mothurOutEndLine(); @@ -195,93 +263,128 @@ void Pintail::getChimeras() { } else { createProcesses(); } - delete distcalculator; - //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 != "") { readQuantiles(); } + if (quanfile != "") { quantiles = readQuantiles(); } else { + + 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(lines[0]->start, lines[0]->end); + quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size()); }else { createProcessesQuan(); } - } - - - //free memory - for (int i = 0; i < lines.size(); i++) { delete lines[i]; } + + ofstream out4, out5; + string noOutliers, outliers; + + if ((!filter) && (seqMask == "")) { + noOutliers = getRootName(templateFile) + "pintail.quan"; + }else if ((filter) && (seqMask == "")) { + noOutliers = getRootName(templateFile) + "pintail.filtered.quan"; + }else if ((!filter) && (seqMask != "")) { + noOutliers = getRootName(templateFile) + "pintail.masked.quan"; + }else if ((filter) && (seqMask != "")) { + noOutliers = getRootName(templateFile) + "pintail.filtered.masked.quan"; + } - } - catch(exception& e) { - errorOut(e, "Pintail", "getChimeras"); - exit(1); - } -} - -//*************************************************************************************************************** - -vector Pintail::readSeqs(string file) { - try { - - ifstream in; - openInputFile(file, in); - vector container; - - //read in seqs and store in vector - while(!in.eof()){ - Sequence* current = new Sequence(in); + //outliers = getRootName(templateFile) + "pintail.quanYESOUTLIERS"; + + /*openOutputFile(outliers, out4); + + //adjust quantiles + 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++) { + temp.push_back(0.0); + } + }else{ + + sort(quantilesMembers[i].begin(), quantilesMembers[i].end(), compareQuanMembers); + + //save 10% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.10)].score); + //save 25% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.25)].score); + //save 50% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.5)].score); + //save 75% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.75)].score); + //save 95% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.95)].score); + //save 99% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.99)].score); + + } + + //output quan value + out4 << i+1 << '\t'; + for (int u = 0; u < temp.size(); u++) { out4 << temp[u] << '\t'; } + out4 << endl; + + quantiles[i] = temp; + + } - //if (current->getAligned() == "") { current->setAligned(current->getUnaligned()); } - //takes out stuff is needed - //current->setUnaligned(current->getUnaligned()); + out4.close();*/ - container.push_back(current); + decalc->removeObviousOutliers(quantilesMembers, templateSeqs.size()); - gobble(in); - } - - in.close(); - return container; - - } - catch(exception& e) { - errorOut(e, "Pintail", "readSeqs"); - exit(1); - } -} + openOutputFile(noOutliers, out5); + + //adjust quantiles + 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++) { + temp.push_back(0.0); + } + }else{ + + sort(quantilesMembers[i].begin(), quantilesMembers[i].end(), compareQuanMembers); + + //save 10% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.10)].score); + //save 25% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.25)].score); + //save 50% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.5)].score); + //save 75% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.75)].score); + //save 95% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.95)].score); + //save 99% + temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.99)].score); + + } + + //output quan value + out5 << i+1 << '\t'; + for (int u = 0; u < temp.size(); u++) { out5 << temp[u] << '\t'; } + out5 << endl; + + quantiles[i] = temp; + + } -//*************************************************************************************************************** -//num is query's spot in querySeqs -map Pintail::trimSeqs(Sequence* query, Sequence subject, int num) { - 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; } - } - - int back = 0; - for (int i = q.length(); i >= 0; i--) { - if (isalpha(q[i]) && isalpha(s[i])) { back = i; break; } - } - - //if num = -1 then you are calling this from quantiles - if (num != -1) { - trimmed[num][front] = back; - return trimmed[num]; + mothurOut("Done."); mothurOutEndLine(); + } - - map temp; - temp[front] = back; - return temp; - + + //free memory + for (int i = 0; i < lines.size(); i++) { delete lines[i]; } + for (int i = 0; i < templateLines.size(); i++) { delete templateLines[i]; } + + delete distcalculator; + delete decalc; } catch(exception& e) { - errorOut(e, "Pintail", "trimSeqs"); + errorOut(e, "Pintail", "getChimeras"); exit(1); } } @@ -295,6 +398,7 @@ vector Pintail::readFreq() { openInputFile(consfile, in); vector prob; + set h = decalc->getPos(); //positions of bases in masking sequence //read in probabilities and store in vector int pos; float num; @@ -303,77 +407,43 @@ vector Pintail::readFreq() { in >> pos >> num; - //do you want this spot - prob.push_back(num); + if (h.count(pos) > 0) { + float Pi; + Pi = (num - 0.25) / 0.75; - gobble(in); - } - - in.close(); - return prob; - - } - catch(exception& e) { - errorOut(e, "Pintail", "readFreq"); - exit(1); - } -} + //cannot have probability less than 0. + if (Pi < 0) { Pi = 0.0; } -//*************************************************************************************************************** - -vector< vector > Pintail::readQuantiles() { - try { - - ifstream in; - openInputFile(quanfile, in); - - vector< vector > quan; - - //read in probabilities and store in vector - int num; float ten, twentyfive, fifty, seventyfive, ninetyfive, ninetynine; - - 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); - - //do you want this spot - quan.push_back(temp); + //do you want this spot + prob.push_back(Pi); + } gobble(in); } in.close(); - return quan; + return prob; } catch(exception& e) { - errorOut(e, "Pintail", "readQuantiles"); + errorOut(e, "Pintail", "readFreq"); exit(1); } } - //*************************************************************************************************************** //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) { +vector Pintail::findPairs(int start, int end) { try { - vector seqsMatches; seqsMatches.resize(end-start); + vector seqsMatches; for(int i = start; i < end; i++){ float smallest = 10000.0; Sequence query = *(querySeqs[i]); - + Sequence* match; + for(int j = 0; j < templateSeqs.size(); j++){ Sequence temp = *(templateSeqs[j]); @@ -382,10 +452,15 @@ vector Pintail::findPairs(int start, int end) { float dist = distcalculator->getDist(); if (dist < smallest) { - seqsMatches[i] = *(templateSeqs[j]); + match = templateSeqs[j]; smallest = dist; } } + + //make copy so trimSeqs doesn't overtrim + Sequence* copy = new Sequence(match->getName(), match->getAligned()); + + seqsMatches.push_back(copy); } return seqsMatches; @@ -397,316 +472,160 @@ vector Pintail::findPairs(int start, int end) { } } -//*************************************************************************************************************** -//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) { +/**************************************************************************************************/ + +void Pintail::createProcessesSpots() { try { +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + int process = 0; + vector processIDS; - 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); - } - - 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 + //loop through and create all the processes you want + while (process != processors) { + int pid = fork(); - seq = query->getAligned(); - for (int m = front; m < (back - size) ; m++) { - - //count number of bases you see - if (isalpha(seq[m])) { countBases++; totalBases++; } + if (pid > 0) { + processIDS.push_back(pid); + process++; + }else if (pid == 0){ - //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 - } + for (int j = lines[process]->start; j < lines[process]->end; j++) { - //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); - } -} + //chops off beginning and end of sequences so they both start and end with a base + map trim; -//*************************************************************************************************************** -vector Pintail::calcObserved(Sequence* query, Sequence subject, vector window, int size) { - try { - - vector temp; - - int startpoint = 0; - for (int m = 0; m < windows.size(); m++) { - - string seqFrag = query->getAligned().substr(window[startpoint], size); - string seqFragsub = subject.getAligned().substr(window[startpoint], size); - - 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; - - temp.push_back(dist); + decalc->trimSeqs(querySeqs[j], bestfit[j], trim); + trimmed[j] = trim; + + } + + 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++) { + it = trimmed[i].begin(); + windowsForeachQuery[i] = decalc->findWindows(querySeqs[i], it->first, it->second, windowSizes[i], increment); + } + mothurOut("Done finding window breaks for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); - startpoint++; - } - - return temp; - } - catch(exception& e) { - errorOut(e, "Pintail", "calcObserved"); - exit(1); - } -} -//*************************************************************************************************************** -float Pintail::calcDist(Sequence* query, Sequence subject, int front, int back) { - try { - - //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; + //write out data to file so parent can read it + ofstream out; + string s = toString(getpid()) + ".temp"; + openOutputFile(s, out); - return dist; - } - catch(exception& e) { - errorOut(e, "Pintail", "calcDist"); - exit(1); - } -} - -//*************************************************************************************************************** -vector Pintail::calcExpected(vector qav, float coef) { - try { - - //for each window - vector queryExpected; + //output windowsForeachQuery + for (int i = lines[process]->start; i < lines[process]->end; i++) { + out << windowsForeachQuery[i].size() << '\t'; + for (int j = 0; j < windowsForeachQuery[i].size(); j++) { + out << windowsForeachQuery[i][j] << '\t'; + } + out << endl; + } - for (int m = 0; m < qav.size(); m++) { + //output windowSizes + for (int i = lines[process]->start; i < lines[process]->end; i++) { + out << windowSizes[i] << '\t'; + } + out << endl; - float expected = qav[m] * coef; + //output trimmed values + for (int i = lines[process]->start; i < lines[process]->end; i++) { + it = trimmed[i].begin(); + out << it->first << '\t' << it->second << endl; + } + out.close(); - queryExpected.push_back(expected); + exit(0); + }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } } - - return queryExpected; - - } - catch(exception& e) { - errorOut(e, "Pintail", "calcExpected"); - exit(1); - } -} -//*************************************************************************************************************** -float Pintail::calcDE(vector obs, vector exp) { - try { - - //for each window - float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2 - for (int m = 0; m < obsDistance.size(); m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } - - float de = sqrt((sum / (obsDistance.size() - 1))); - - return de; - } - catch(exception& e) { - errorOut(e, "Pintail", "calcDE"); - exit(1); - } -} - -//*************************************************************************************************************** - -vector Pintail::calcFreq(vector seqs) { - try { - - vector prob; - string freqfile = getRootName(templateFile) + "probability"; - ofstream outFreq; - openOutputFile(freqfile, outFreq); + //force parent to wait until all the processes are done + for (int i=0;igetAligned().length(); i++) { + //get data created by processes + for (int i=0;i freq; freq.resize(4,0); - int gaps = 0; + int size = lines[i]->end - lines[i]->start; + + int count = lines[i]->start; + for (int m = 0; m < size; m++) { + int num; + in >> num; + + vector win; int w; + for (int j = 0; j < num; j++) { + in >> w; + win.push_back(w); + } - //find the frequency of each nucleotide - for (int j = 0; j < seqs.size(); j++) { - - char value = seqs[j]->getAligned()[i]; + windowsForeachQuery[count] = win; + count++; + gobble(in); + } + + gobble(in); + count = lines[i]->start; + for (int m = 0; m < size; m++) { + int num; + in >> num; - 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++; } + windowSizes[count] = num; + count++; } - //find base with highest frequency - int highest = 0; - for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } - - float highFreq; - if ( (seqs.size() - gaps) == 0 ) { highFreq = 1.0; } - else { highFreq = highest / (float) (seqs.size() - gaps); } - //highFreq = highest / (float) seqs.size(); - //cout << i << '\t' << highFreq << endl; - - float Pi; - Pi = (highFreq - 0.25) / 0.75; - - //cannot have probability less than 0. - if (Pi < 0) { Pi = 0.0; } - - //saves this for later - outFreq << i+1 << '\t' << Pi << endl; + gobble(in); - prob.push_back(Pi); - } - - 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 < windows.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; - - //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; - - //for each sequence - for(int i = start; i < end; i++){ - - Sequence* query = templateSeqs[i]; + count = lines[i]->start; + for (int m = 0; m < size; m++) { + int front, back; + in >> front >> back; - //compare to every other sequence in template - for(int j = 0; j < i; j++){ + map t; - Sequence subject = *(templateSeqs[j]); + t[front] = back; - map trim = trimSeqs(query, subject, -1); + trimmed[count] = t; + count++; - - - - + gobble(in); } + + in.close(); + remove(s.c_str()); + } - + +#else + for (int j = 0; j < bestfit.size(); j++) { + //chops off beginning and end of sequences so they both start and end with a base + decalc->trimSeqs(querySeqs[j], bestfit[j], trimmed[j]); } - return quan; - - } - catch(exception& e) { - errorOut(e, "Pintail", "findQav"); - exit(1); - } -} -//*************************************************************************************************************** -float Pintail::getCoef(vector obs, vector qav) { - 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(); - + for (int i = lines[0]->start; i < lines[0]->end; i++) { + it = trimmed[i].begin(); + vector win = decalc->findWindows(querySeqs[i], it->first, it->second, windowSizes[i], increment); + windowsForeachQuery[i] = win; + } - float coef = obsAverage / probAverage; - - return coef; +#endif } catch(exception& e) { - errorOut(e, "Pintail", "getCoef"); + errorOut(e, "Pintail", "createProcessesSpots"); exit(1); } } /**************************************************************************************************/ -void Pintail::createProcessesSpots() { +void Pintail::createProcessesPairs() { try { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) int process = 0; vector processIDS; - vector< vector > win; win.resize(querySeqs.size()); //loop through and create all the processes you want while (process != processors) { @@ -717,22 +636,23 @@ void Pintail::createProcessesSpots() { process++; }else if (pid == 0){ - 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], i); - count++; - } + mothurOut("Finding pairs for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); + bestfit = findPairs(lines[process]->start, lines[process]->end); + mothurOut("Done finding pairs for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); + //write out data to file so parent can read it + ofstream out; + string s = toString(getpid()) + ".temp"; + openOutputFile(s, out); - 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; + //output range and size + out << bestfit.size() << endl; + + //output pairs + for (int i = 0; i < bestfit.size(); i++) { + out << ">" << bestfit[i]->getName() << endl << bestfit[i]->getAligned() << endl; } + out.close(); exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } @@ -744,29 +664,39 @@ void Pintail::createProcessesSpots() { wait(&temp); } - windows = win; + //get data created by processes + for (int i=0;i> size; gobble(in); + + //get pairs + int count = lines[i]->start; + for (int m = 0; m < size; m++) { + Sequence* temp = new Sequence(in); + bestfit[count] = temp; + + count++; + gobble(in); + } + + in.close(); + remove(s.c_str()); + } + + #else bestfit = findPairs(lines[0]->start, 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; - } - #endif } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesSpots"); + errorOut(e, "Pintail", "createProcessesPairs"); exit(1); } } - - /**************************************************************************************************/ void Pintail::createProcesses() { @@ -775,12 +705,6 @@ void Pintail::createProcesses() { int process = 0; vector processIDS; - 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()); - - //loop through and create all the processes you want while (process != processors) { int pid = fork(); @@ -790,30 +714,74 @@ void Pintail::createProcesses() { process++; }else if (pid == 0){ - //calc obs + 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], windows[i], windowSizes[i]); - obs[i] = obsi; + + vector obsi = decalc->calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); + obsDistance[i] = obsi; //calc Qav - vector q = findQav(windows[i], windowSizes[i]); + vector q = decalc->findQav(windowsForeachQuery[i], windowSizes[i], probabilityProfile); //get alpha - float alpha = getCoef(obsDistance[i], q); + float alpha = decalc->getCoef(obsDistance[i], q); //find expected - vector exp = calcExpected(q, alpha); + vector exp = decalc->calcExpected(q, alpha); expectedDistance[i] = exp; //get de and deviation - float dei = calcDE(obsi, exp); - de[i] = dei; + float dei = decalc->calcDE(obsi, exp); + DE[i] = dei; it = trimmed[i].begin(); - float dist = calcDist(querySeqs[i], bestfit[i], it->first, it->second); - dev[i] = dist; + float dist = decalc->calcDist(querySeqs[i], bestfit[i], it->first, it->second); + deviation[i] = dist; + } + mothurOut("Done calculating observed, expected and de values for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); + + //write out data to file so parent can read it + ofstream out; + string s = toString(getpid()) + ".temp"; + openOutputFile(s, out); + + int size = lines[process]->end - lines[process]->start; + out << size << endl; + + //output observed distances + for (int i = lines[process]->start; i < lines[process]->end; i++) { + out << obsDistance[i].size() << '\t'; + for (int j = 0; j < obsDistance[i].size(); j++) { + out << obsDistance[i][j] << '\t'; + } + out << endl; + } + + + //output expected distances + for (int i = lines[process]->start; i < lines[process]->end; i++) { + out << expectedDistance[i].size() << '\t'; + for (int j = 0; j < expectedDistance[i].size(); j++) { + out << expectedDistance[i][j] << '\t'; + } + out << endl; + } + + + //output de values + for (int i = lines[process]->start; i < lines[process]->end; i++) { + out << DE[i] << '\t'; } + out << endl; + //output de values + for (int i = lines[process]->start; i < lines[process]->end; i++) { + out << deviation[i] << '\t'; + } + out << endl; + + out.close(); + exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } } @@ -824,15 +792,82 @@ void Pintail::createProcesses() { wait(&temp); } - obsDistance = obs; - expectedDistance = exp; - DE = de; - deviation = dev; - + //get data created by processes + for (int i=0;i> size; gobble(in); + + //get observed distances + int count = lines[i]->start; + for (int m = 0; m < size; m++) { + int num; + in >> num; + + vector obs; float w; + for (int j = 0; j < num; j++) { + in >> w; + obs.push_back(w); + } + + obsDistance[count] = obs; + count++; + gobble(in); + } + + gobble(in); + + //get expected distances + count = lines[i]->start; + for (int m = 0; m < size; m++) { + int num; + in >> num; + + vector exp; float w; + for (int j = 0; j < num; j++) { + in >> w; + exp.push_back(w); + } + + expectedDistance[count] = exp; + count++; + gobble(in); + } + + gobble(in); + + count = lines[i]->start; + for (int m = 0; m < size; m++) { + float num; + in >> num; + + DE[count] = num; + count++; + } + + gobble(in); + + count = lines[i]->start; + for (int m = 0; m < size; m++) { + float num; + in >> num; + + deviation[count] = num; + count++; + } + + in.close(); + remove(s.c_str()); + } + + #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]); + vector obsi = decalc->calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); obsDistance[i] = obsi; } mothurOut("Done."); mothurOutEndLine(); @@ -841,7 +876,7 @@ void Pintail::createProcesses() { mothurOut("Finding variability... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector q = findQav(windows[i], windowSizes[i]); + vector q = decalc->findQav(windowsForeachQuery[i], windowSizes[i], probabilityProfile); Qav[i] = q; } mothurOut("Done."); mothurOutEndLine(); @@ -850,7 +885,7 @@ void Pintail::createProcesses() { mothurOut("Calculating alpha... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { - float alpha = getCoef(obsDistance[i], Qav[i]); + float alpha = decalc->getCoef(obsDistance[i], Qav[i]); seqCoef.push_back(alpha); } mothurOut("Done."); mothurOutEndLine(); @@ -859,7 +894,7 @@ void Pintail::createProcesses() { mothurOut("Calculating expected distance... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector exp = calcExpected(Qav[i], seqCoef[i]); + vector exp = decalc->calcExpected(Qav[i], seqCoef[i]); expectedDistance[i] = exp; } mothurOut("Done."); mothurOutEndLine(); @@ -868,11 +903,11 @@ void Pintail::createProcesses() { mothurOut("Finding deviation... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { - float de = calcDE(obsDistance[i], expectedDistance[i]); + float de = decalc->calcDE(obsDistance[i], expectedDistance[i]); DE[i] = de; it = trimmed[i].begin(); - float dist = calcDist(querySeqs[i], bestfit[i], it->first, it->second); + float dist = decalc->calcDist(querySeqs[i], bestfit[i], it->first, it->second); deviation[i] = dist; } mothurOut("Done."); mothurOutEndLine(); @@ -903,6 +938,24 @@ void Pintail::createProcessesQuan() { process++; }else if (pid == 0){ + quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, templateLines[process]->start, templateLines[process]->end); + + //write out data to file so parent can read it + ofstream out; + string s = toString(getpid()) + ".temp"; + 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].score << '\t' << quantilesMembers[i][j].member1 << '\t' << quantilesMembers[i][j].member2 << '\t'; + } + out << endl; + } + + out.close(); exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } @@ -913,9 +966,50 @@ void Pintail::createProcessesQuan() { int temp = processIDS[i]; wait(&temp); } + + //get data created by processes + for (int i=0;i > quan; + quan.resize(100); + + //get quantiles + for (int m = 0; m < quan.size(); m++) { + int num; + in >> num; + + gobble(in); + + vector q; float w; int b, n; + for (int j = 0; j < num; j++) { + in >> w >> b >> n; + //cout << w << '\t' << b << '\t' n << endl; + quanMember newMember(w, b, n); + q.push_back(newMember); + } +//cout << "here" << endl; + quan[m] = q; +//cout << "now here" << endl; + 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(); + remove(s.c_str()); + } #else - + quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size()); #endif } catch(exception& e) {