X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=pintail.cpp;h=90c97c5e77d50f478aeb1b4bb7b1fa28713ee468;hb=818dcf9513122fa52616a7722f4893f77685a621;hp=7ae3ff6dbb96224197dff97db4b435880738d87b;hpb=98ea55ba2d46a429f031086b4b3272780d0ec894;p=mothur.git diff --git a/pintail.cpp b/pintail.cpp index 7ae3ff6..90c97c5 100644 --- a/pintail.cpp +++ b/pintail.cpp @@ -18,104 +18,40 @@ inline bool compareQuanMembers(quanMember left, quanMember right){ } //*************************************************************************************************************** -Pintail::Pintail(string filename, string temp) { fastafile = filename; templateFile = temp; } +Pintail::Pintail(string filename, string o) { + fastafile = filename; outputDir = o; + distcalculator = new eachGapDist(); + decalc = new DeCalculator(); +} //*************************************************************************************************************** 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"); - exit(1); - } -} -//*************************************************************************************************************** -void Pintail::print(ostream& out) { try { - mothurOutEndLine(); - - for (int i = 0; i < querySeqs.size(); i++) { - - 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'; } - 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"); + errorOut(e, "Pintail", "~Pintail"); exit(1); } } - //*************************************************************************************************************** -void Pintail::getChimeras() { +void 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); - 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 = numSeqs / processors ; + //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()); } + + decalc->setMask(seqMask); #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)); - } - //find breakup of templatefile for quantiles if (processors == 1) { templateLines.push_back(new linePair(0, templateSeqs.size())); } else { @@ -126,148 +62,78 @@ void Pintail::getChimeras() { } } #else - lines.push_back(new linePair(0, numSeqs)); templateLines.push_back(new linePair(0, templateSeqs.size())); #endif - - 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); - 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); + probabilityProfile = decalc->calcFreq(templateSeqs, outputDir + getSimpleName(templateFileName)); mothurOut("Done."); mothurOutEndLine(); - }else { probabilityProfile = readFreq(); } - + }else { probabilityProfile = readFreq(); mothurOut("Done."); } + mothurOutEndLine(); + //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]); - } - - } + bool reRead = false; + //create filter if needed for later if (filter) { - vector temp = templateSeqs; - for (int i = 0; i < querySeqs.size(); i++) { temp.push_back(querySeqs[i]); } - - createFilter(temp); + + //read in all query seqs + ifstream in; + openInputFile(fastafile, in); - runFilter(querySeqs); - runFilter(templateSeqs); - runFilter(bestfit); - } - - - if (processors == 1) { - - for (int j = 0; j < bestfit.size(); j++) { - decalc->trimSeqs(querySeqs[j], bestfit[j], trimmed[j]); + vector tempQuerySeqs; + while(!in.eof()){ + Sequence* s = new Sequence(in); + gobble(in); + + if (s->getName() != "") { tempQuerySeqs.push_back(s); } } + in.close(); - mothurOut("Finding window breaks... "); cout.flush(); - 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; - } - mothurOut("Done."); mothurOutEndLine(); - - }else { createProcessesSpots(); } - - if (processors == 1) { - - mothurOut("Calculating observed distance... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector obsi = decalc->calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); + vector temp; + //merge query seqs and template seqs + temp = templateSeqs; + for (int i = 0; i < tempQuerySeqs.size(); i++) { temp.push_back(tempQuerySeqs[i]); } - obsDistance[i] = obsi; + if (seqMask != "") { + reRead = true; + //mask templates + for (int i = 0; i < temp.size(); i++) { + decalc->runMask(temp[i]); + } } - mothurOut("Done."); mothurOutEndLine(); - - - mothurOut("Finding variability... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector q = decalc->findQav(windowsForeachQuery[i], windowSizes[i], probabilityProfile); - Qav[i] = q; - } - mothurOut("Done."); mothurOutEndLine(); - + mergedFilterString = createFilter(temp, 0.5); - mothurOut("Calculating alpha... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; i++) { - float alpha = decalc->getCoef(obsDistance[i], Qav[i]); - seqCoef[i] = alpha; - } - mothurOut("Done."); mothurOutEndLine(); + //reread template seqs + for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i]; } + } - mothurOut("Calculating expected distance... "); cout.flush(); - for (int i = lines[0]->start; i < lines[0]->end; 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 = decalc->calcDE(obsDistance[i], expectedDistance[i]); - DE[i] = de; - - it = trimmed[i].begin(); - float dist = decalc->calcDist(querySeqs[i], bestfit[i], it->first, it->second); - deviation[i] = dist; - } - mothurOut("Done."); mothurOutEndLine(); - - } - else { createProcesses(); } - //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++) { + decalc->runMask(templateSeqs[i]); + } + } + + if (filter) { + reRead = true; + for (int i = 0; i < templateSeqs.size(); i++) { + runFilter(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 (processors == 1) { @@ -279,62 +145,21 @@ void Pintail::getChimeras() { string noOutliers, outliers; if ((!filter) && (seqMask == "")) { - noOutliers = getRootName(templateFile) + "pintail.quan"; - }else if ((filter) && (seqMask == "")) { - noOutliers = getRootName(templateFile) + "pintail.filtered.quan"; + noOutliers = outputDir + getRootName(getSimpleName(templateFileName)) + "pintail.quan"; }else if ((!filter) && (seqMask != "")) { - noOutliers = getRootName(templateFile) + "pintail.masked.quan"; + noOutliers = outputDir + getRootName(getSimpleName(templateFileName)) + "pintail.masked.quan"; }else if ((filter) && (seqMask != "")) { - noOutliers = getRootName(templateFile) + "pintail.filtered.masked.quan"; + noOutliers = outputDir + getRootName(getSimpleName(templateFileName)) + "pintail.filtered." + getSimpleName(getRootName(fastafile)) + "masked.quan"; + }else if ((filter) && (seqMask == "")) { + noOutliers = outputDir + getRootName(getSimpleName(templateFileName)) + "pintail.filtered." + getSimpleName(getRootName(fastafile)) + "quan"; } - //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; - - } - - out4.close();*/ - + + + decalc->removeObviousOutliers(quantilesMembers, templateSeqs.size()); - openOutputFile(noOutliers, out5); - + openOutputFile(noOutliers, out5); //adjust quantiles for (int i = 0; i < quantilesMembers.size(); i++) { vector temp; @@ -373,15 +198,113 @@ void Pintail::getChimeras() { } mothurOut("Done."); mothurOutEndLine(); - } - + + if (reRead) { + for (int i = 0; i < templateSeqs.size(); i++) { delete templateSeqs[i]; } + templateSeqs.clear(); + templateSeqs = readSeqs(templateFileName); + } + + //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; + for (int i = 0; i < templateLines.size(); i++) { delete templateLines[i]; } + + } + catch(exception& e) { + errorOut(e, "Pintail", "doPrep"); + exit(1); + } +} +//*************************************************************************************************************** +void Pintail::print(ostream& out, ostream& outAcc) { + try { + 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"; } + + out << querySeq->getName() << '\t' << "div: " << deviation << "\tstDev: " << DE << "\tchimera flag: " << chimera << endl; + if (chimera == "Yes") { + mothurOut(querySeq->getName() + "\tdiv: " + toString(deviation) + "\tstDev: " + toString(DE) + "\tchimera flag: " + chimera); mothurOutEndLine(); + outAcc << querySeq->getName() << endl; + } + out << "Observed\t"; + + for (int j = 0; j < obsDistance.size(); j++) { out << obsDistance[j] << '\t'; } + out << endl; + + out << "Expected\t"; + + for (int m = 0; m < expectedDistance.size(); m++) { out << expectedDistance[m] << '\t'; } + out << endl; + + + } + catch(exception& e) { + errorOut(e, "Pintail", "print"); + exit(1); + } +} + +//*************************************************************************************************************** +int Pintail::getChimeras(Sequence* query) { + try { + querySeq = query; + trimmed.clear(); + windowSizes = window; + + //find pairs has to be done before a mask + bestfit = findPairs(query); + + //if they mask + if (seqMask != "") { + decalc->runMask(query); + decalc->runMask(bestfit); + } + + if (filter) { //must be done after a mask + runFilter(query); + runFilter(bestfit); + } + + + //trim seq + decalc->trimSeqs(query, bestfit, trimmed); + + //find windows + it = trimmed.begin(); + windowsForeachQuery = decalc->findWindows(query, it->first, it->second, windowSizes, increment); + + //find observed distance + obsDistance = decalc->calcObserved(query, bestfit, windowsForeachQuery, windowSizes); + + Qav = decalc->findQav(windowsForeachQuery, windowSizes, probabilityProfile); + + //find alpha + seqCoef = decalc->getCoef(obsDistance, Qav); + + //calculating expected distance + expectedDistance = decalc->calcExpected(Qav, seqCoef); + + //finding de + DE = decalc->calcDE(obsDistance, expectedDistance); + + //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", "getChimeras"); @@ -433,36 +356,12 @@ vector Pintail::readFreq() { //*************************************************************************************************************** //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::findPairs(Sequence* q) { try { - 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]); - - distcalculator->calcDist(query, temp); - float dist = distcalculator->getDist(); - - if (dist < smallest) { - match = templateSeqs[j]; - smallest = dist; - } - } - - //make copy so trimSeqs doesn't overtrim - Sequence* copy = new Sequence(match->getName(), match->getAligned()); - - seqsMatches.push_back(copy); - } + Sequence* seqsMatches; + seqsMatches = decalc->findClosest(q, templateSeqs); return seqsMatches; } @@ -471,458 +370,7 @@ vector Pintail::findPairs(int start, int end) { exit(1); } } - /**************************************************************************************************/ - -void Pintail::createProcessesSpots() { - try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector processIDS; - - //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){ - - for (int j = lines[process]->start; j < lines[process]->end; j++) { - - //chops off beginning and end of sequences so they both start and end with a base - map trim; - - 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(); - - //write out data to file so parent can read it - ofstream out; - string s = toString(getpid()) + ".temp"; - openOutputFile(s, out); - - //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; - } - - //output windowSizes - for (int i = lines[process]->start; i < lines[process]->end; i++) { - out << windowSizes[i] << '\t'; - } - out << endl; - - //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(); - - 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;iend - 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); - } - - windowsForeachQuery[count] = win; - count++; - gobble(in); - } - - gobble(in); - count = lines[i]->start; - for (int m = 0; m < size; m++) { - int num; - in >> num; - - windowSizes[count] = num; - count++; - } - - gobble(in); - - count = lines[i]->start; - for (int m = 0; m < size; m++) { - int front, back; - in >> front >> back; - - map t; - - t[front] = back; - - 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]); - } - - 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; - } - -#endif - } - catch(exception& e) { - errorOut(e, "Pintail", "createProcessesSpots"); - exit(1); - } -} -/**************************************************************************************************/ - -void Pintail::createProcessesPairs() { - try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector processIDS; - - //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("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); - - //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); } - } - - //force parent to wait until all the processes are done - 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); -#endif - } - catch(exception& e) { - errorOut(e, "Pintail", "createProcessesPairs"); - exit(1); - } -} -/**************************************************************************************************/ - -void Pintail::createProcesses() { - try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector processIDS; - - //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 = decalc->calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); - obsDistance[i] = obsi; - - //calc Qav - vector q = decalc->findQav(windowsForeachQuery[i], windowSizes[i], probabilityProfile); - - //get alpha - float alpha = decalc->getCoef(obsDistance[i], q); - - //find expected - vector exp = decalc->calcExpected(q, alpha); - expectedDistance[i] = exp; - - //get de and deviation - float dei = decalc->calcDE(obsi, exp); - DE[i] = dei; - - it = trimmed[i].begin(); - 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); } - } - - //force parent to wait until all the processes are done - 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 = 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 = 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 = decalc->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 = 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 = decalc->calcDE(obsDistance[i], expectedDistance[i]); - DE[i] = de; - - it = trimmed[i].begin(); - float dist = decalc->calcDist(querySeqs[i], bestfit[i], it->first, it->second); - deviation[i] = dist; - } - mothurOut("Done."); mothurOutEndLine(); - -#endif - } - catch(exception& e) { - errorOut(e, "Pintail", "createProcesses"); - exit(1); - } -} - - -/**************************************************************************************************/ - void Pintail::createProcessesQuan() { try { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)