X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=pintail.cpp;h=a204294f515ee17f4e8e62f0045d8c67dfe6e406;hb=dbd5da8043df1cb9f5ff7c6ddb5f550ea49b52c2;hp=4b812a7c628d2816a0d8215cf0b80862b5471626;hpb=0f297d3c211acbdfd4d941ef193a6702c5046dd1;p=mothur.git diff --git a/pintail.cpp b/pintail.cpp index 4b812a7..a204294 100644 --- a/pintail.cpp +++ b/pintail.cpp @@ -12,15 +12,7 @@ //*************************************************************************************************************** -Pintail::Pintail(string name) { - try { - fastafile = name; - } - catch(exception& e) { - errorOut(e, "Pintail", "Pintail"); - exit(1); - } -} +Pintail::Pintail(string filename, string temp) { fastafile = filename; templateFile = temp; } //*************************************************************************************************************** Pintail::~Pintail() { @@ -37,24 +29,27 @@ Pintail::~Pintail() { void Pintail::print(ostream& out) { try { - for (itCoef = DE.begin(); itCoef != DE.end(); itCoef++) { + 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 << itCoef->first->getName() << '\t' << itCoef->second << endl; + out << querySeqs[i]->getName() << '\t' << "div: " << deviation[i] << "\tstDev: " << DE[i] << "\tchimera flag: " << chimera << endl; out << "Observed\t"; - itObsDist = obsDistance.find(itCoef->first); - for (int i = 0; i < itObsDist->second.size(); i++) { out << itObsDist->second[i] << '\t'; } + for (int j = 0; j < obsDistance[i].size(); j++) { out << obsDistance[i][j] << '\t'; } out << endl; out << "Expected\t"; - itExpDist = expectedDistance.find(itCoef->first); - for (int i = 0; i < itExpDist->second.size(); i++) { out << itExpDist->second[i] << '\t'; } + for (int m = 0; m < expectedDistance[i].size(); m++) { out << expectedDistance[i][m] << '\t'; } out << endl; } - - } catch(exception& e) { errorOut(e, "Pintail", "print"); @@ -66,8 +61,6 @@ void Pintail::print(ostream& out) { void Pintail::getChimeras() { try { - distCalculator = new ignoreGaps(); - //read in query sequences and subject sequences mothurOut("Reading sequences and template file... "); cout.flush(); querySeqs = readSeqs(fastafile); @@ -76,116 +69,330 @@ void Pintail::getChimeras() { int numSeqs = querySeqs.size(); - //if window is set to default - if (window == 0) { if (querySeqs[0]->getAligned().length() > 800) { setWindow(200); } - else{ setWindow((querySeqs[0]->getAligned().length() / 4)); } } - else if (window > (querySeqs[0]->getAligned().length() / 4)) { - mothurOut("You have selected to large a window size for you sequences. I will choose a smaller window."); mothurOutEndLine(); - setWindow((querySeqs[0]->getAligned().length() / 4)); - } - - //calculate number of iters - iters = (querySeqs[0]->getAligned().length() - window + 1) / increment; -cout << "length = " << querySeqs[0]->getAligned().length() << " window = " << window << " increment = " << increment << " iters = " << iters << endl; + 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); + windowSizesTemplate.resize(templateSeqs.size(), window); + windowsForeachQuery.resize(numSeqs); + quantiles.resize(100); //one for every percent mismatch + + //break up file if needed int linesPerProcess = processors / numSeqs; - //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))); + #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)); } - //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)); - } - - //map query sequences to their most similiar sequences in the template - Parallelized - mothurOut("Finding closest sequence in template to each sequence... "); cout.flush(); - if (processors == 1) { findPairs(lines[0]->start, lines[0]->end); } - else { createProcessesPairs(); } - mothurOut("Done."); mothurOutEndLine(); + + //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 - //find Oqs for each sequence - the observed distance in each window - Parallelized - mothurOut("Calculating observed percentage differences for each sequence... "); cout.flush(); - if (processors == 1) { calcObserved(lines[0]->start, lines[0]->end); } - else { createProcessesObserved(); } - mothurOut("Done."); mothurOutEndLine(); + 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("Calculating expected percentage differences for each sequence... "); cout.flush(); - vector probabilityProfile = calcFreq(templateSeqs); + 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(); } //make P into Q for (int i = 0; i < probabilityProfile.size(); i++) { probabilityProfile[i] = 1 - probabilityProfile[i]; } - - //find Qav - averageProbability = findQav(probabilityProfile); - - //find Coefficient - maps a sequence to its coefficient - seqCoef = getCoef(averageProbability); - - //find Eqs for each sequence - the expected distance in each window - Parallelized - if (processors == 1) { calcExpected(lines[0]->start, lines[0]->end); } - else { createProcessesExpected(); } mothurOut("Done."); mothurOutEndLine(); - //find deviation - Parallelized - mothurOut("Finding deviation from expected... "); cout.flush(); - if (processors == 1) { calcDE(lines[0]->start, lines[0]->end); } - else { createProcessesDE(); } - mothurOut("Done."); mothurOutEndLine(); + if (processors == 1) { + + 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; + } + 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]); + DE[i] = de; - //free memory - for (int i = 0; i < lines.size(); i++) { delete lines[i]; } - delete distCalculator; + it = trimmed[i].begin(); + float dist = 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 { + + 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()); + }else { createProcessesQuan(); } + + ofstream out4; + string o = getRootName(templateFile) + "quan"; + + openOutputFile(o, out4); + + //adjust quantiles + for (int i = 0; i < quantiles.size(); i++) { + if (quantiles[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); + } + }else{ + + sort(quantiles[i].begin(), quantiles[i].end()); + + vector temp; + //save 10% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.10)]); + //save 25% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.25)]); + //save 50% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.5)]); + //save 75% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.75)]); + //save 95% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.95)]); + //save 99% + temp.push_back(quantiles[i][int(quantiles[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; + } + + 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]; } + + 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; } + } + + int back = 0; + for (int i = q.length(); i >= 0; i--) { + if (isalpha(q[i]) && isalpha(s[i])) { back = i; break; } + } + + trim[front] = back; + + } + catch(exception& e) { + errorOut(e, "Pintail", "trimSeqs"); + exit(1); + } +} //*************************************************************************************************************** -vector Pintail::readSeqs(string file) { +vector Pintail::readFreq() { try { ifstream in; - openInputFile(file, in); - vector container; + openInputFile(consfile, in); + + vector prob; + + //read in probabilities and store in vector + int pos; float num; - //read in seqs and store in vector while(!in.eof()){ - Sequence* current = new Sequence(in); - if (current->getAligned() == "") { current->setAligned(current->getUnaligned()); } - //takes out stuff is needed - current->setUnaligned(current->getUnaligned()); + in >> pos >> num; - container.push_back(current); + //do you want this spot + prob.push_back(num); gobble(in); } in.close(); - return container; + return prob; } catch(exception& e) { - errorOut(e, "Pintail", "readSeqs"); + errorOut(e, "Pintail", "readFreq"); exit(1); } } +//*************************************************************************************************************** + +vector< vector > Pintail::readQuantiles() { + try { + + ifstream in; + openInputFile(quanfile, in); + + vector< vector > quan; + + 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); + + quan.push_back(temp); + + gobble(in); + } + + in.close(); + return quan; + + } + catch(exception& e) { + errorOut(e, "Pintail", "readQuantiles"); + exit(1); + } +} //*************************************************************************************************************** //calculate the distances from each query sequence to all sequences in the template to find the closest sequence -void Pintail::findPairs(int start, int end) { +vector Pintail::findPairs(int start, int end) { try { + vector seqsMatches; seqsMatches.resize(end-start); + for(int i = start; i < end; i++){ float smallest = 10000.0; @@ -195,124 +402,170 @@ void Pintail::findPairs(int start, int end) { Sequence temp = *(templateSeqs[j]); - distCalculator->calcDist(query, temp); - float dist = distCalculator->getDist(); + distcalculator->calcDist(query, temp); + float dist = distcalculator->getDist(); if (dist < smallest) { - - bestfit[querySeqs[i]] = templateSeqs[j]; + seqsMatches[i] = *(templateSeqs[j]); smallest = dist; } } } + return seqsMatches; + } catch(exception& e) { errorOut(e, "Pintail", "findPairs"); exit(1); } } + //*************************************************************************************************************** -void Pintail::calcObserved(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) { try { - - - for(int i = start; i < end; i++){ - itBest = bestfit.find(querySeqs[i]); - Sequence* query; - Sequence* subject; + vector win; - if (itBest != bestfit.end()) { - query = itBest->first; - subject = itBest->second; - }else{ mothurOut("Error in calcObserved"); mothurOutEndLine(); } -//cout << query->getName() << '\t' << subject->getName() << endl; + int cutoff = back - front; //back - front - int startpoint = 0; - for (int m = 0; m < iters; m++) { - - string seqFrag = query->getAligned().substr(startpoint, window); - string seqFragsub = subject->getAligned().substr(startpoint, window); - - int diff = 0; - for (int b = 0; b < seqFrag.length(); b++) { - - //if this is not a gap - if ((isalpha(seqFrag[b])) && (isalpha(seqFragsub[b]))) { - //and they are different - penalize - if (seqFrag[b] != seqFragsub[b]) { diff++; } - } - } - - //percentage of mismatched bases - float dist = diff / (float)seqFrag.length(); + //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 + + seq = query->getAligned(); + for (int m = front; m < (back - size) ; m++) { - obsDistance[query].push_back(dist); + //count number of bases you see + if (isalpha(seq[m])) { countBases++; totalBases++; } - startpoint += increment; + //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", "calcObserved"); + errorOut(e, "Pintail", "findWindows"); exit(1); } } //*************************************************************************************************************** -void Pintail::calcExpected(int start, int end) { +vector Pintail::calcObserved(Sequence* query, Sequence subject, vector window, int size) { try { - - //for each sequence - for(int i = start; i < end; i++){ - - itCoef = seqCoef.find(querySeqs[i]); - float coef = itCoef->second; - - //for each window - vector queryExpected; - for (int m = 0; m < iters; m++) { - float expected = averageProbability[m] * coef; - queryExpected.push_back(expected); -//cout << "average variabilty over window = " << averageProbability[m] << " coef = " << coef << " ei = " << expected << '\t' << "window = " << m << endl; + 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; + + temp.push_back(dist); + } - expectedDistance[querySeqs[i]] = queryExpected; + 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; + return dist; } catch(exception& e) { - errorOut(e, "Pintail", "calcExpected"); + errorOut(e, "Pintail", "calcDist"); exit(1); } } + //*************************************************************************************************************** -void Pintail::calcDE(int start, int end) { +vector Pintail::calcExpected(vector qav, float coef) { try { - - //for each sequence - for(int i = start; i < end; i++){ + //for each window + vector queryExpected; - itObsDist = obsDistance.find(querySeqs[i]); - vector obs = itObsDist->second; + for (int m = 0; m < qav.size(); m++) { + + float expected = qav[m] * coef; + + queryExpected.push_back(expected); + } - itExpDist = expectedDistance.find(querySeqs[i]); - vector exp = itExpDist->second; -// cout << "difference between obs and exp = " << abs(obs[m] - exp[m]) << endl; - //for each window - float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2 - for (int m = 0; m < iters; m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } + 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 < obs.size(); m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } - float de = sqrt((sum / (iters - 1))); + float de = sqrt((sum / (obs.size() - 1))); - DE[querySeqs[i]] = de; - } - + return de; } catch(exception& e) { errorOut(e, "Pintail", "calcDE"); @@ -326,10 +579,14 @@ vector Pintail::calcFreq(vector seqs) { try { vector prob; + string freqfile = getRootName(templateFile) + "prob"; + ofstream outFreq; + + openOutputFile(freqfile, outFreq); //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; @@ -337,7 +594,7 @@ vector Pintail::calcFreq(vector seqs) { 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]++; } @@ -347,19 +604,27 @@ vector Pintail::calcFreq(vector seqs) { //find base with highest frequency int highest = 0; - for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } + for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } + + 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; - //add in gaps - so you can effectively "ignore them" - highest += gaps; + //cannot have probability less than 0. + if (Pi < 0) { Pi = 0.0; } - float highFreq = highest / (float) seqs.size(); + //saves this for later + outFreq << i+1 << '\t' << Pi << endl; - float Pi; - Pi = (highFreq - 0.25) / 0.75; - prob.push_back(Pi); } + outFreq.close(); + return prob; } @@ -369,25 +634,24 @@ vector Pintail::calcFreq(vector seqs) { } } //*************************************************************************************************************** -vector Pintail::findQav(vector prob) { +vector Pintail::findQav(vector window, int size) { try { - vector averages; - + vector averages; + //for each window find average - int startpoint = 0; - for (int m = 0; m < iters; m++) { - + for (int m = 0; m < window.size(); m++) { + float average = 0.0; - for (int i = startpoint; i < (startpoint+window); i++) { average += prob[i]; } - - average = average / window; -//cout << average << endl; + + //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); - - startpoint += increment; } - + return averages; } catch(exception& e) { @@ -395,88 +659,99 @@ vector Pintail::findQav(vector prob) { exit(1); } } + //*************************************************************************************************************** -map Pintail::getCoef(vector prob) { +vector< vector > Pintail::getQuantiles(int start, int end) { try { - map coefs; + vector< vector > quan; - //find average prob - float probAverage = 0.0; - for (int i = 0; i < prob.size(); i++) { probAverage += prob[i]; } - probAverage = probAverage / (float) prob.size(); -cout << "(sum of ai) / m = " << probAverage << endl; - //find a coef for each sequence - map >::iterator it; - for (it = obsDistance.begin(); it != obsDistance.end(); it++) { - - vector temp = it->second; - Sequence* tempSeq = it->first; - - //find observed average - float obsAverage = 0.0; - for (int i = 0; i < temp.size(); i++) { obsAverage += temp[i]; } - obsAverage = obsAverage / (float) temp.size(); -cout << tempSeq->getName() << '\t' << obsAverage << endl; - float coef = obsAverage / probAverage; -cout << tempSeq->getName() << '\t' << "coef = " << coef << endl; - //save this sequences coefficient - coefs[tempSeq] = coef; - } + //percentage of mismatched pairs 1 to 100 + quan.resize(100); + + //for each sequence + for(int i = start; i < end; i++){ + + 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; - return coefs; } catch(exception& e) { - errorOut(e, "Pintail", "getCoef"); + errorOut(e, "Pintail", "findQav"); exit(1); } } -/**************************************************************************************************/ - -void Pintail::createProcessesPairs() { +//*************************************************************************************************************** +float Pintail::getCoef(vector obs, vector qav) { try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector processIDS; + + //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(); - //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){ - findPairs(lines[process]->start, lines[process]->end); - exit(0); - }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } - } + //find observed average + float obsAverage = 0.0; + for (int j = 0; j < obs.size(); j++) { obsAverage += obs[j]; } + obsAverage = obsAverage / (float) obs.size(); - //force parent to wait until all the processes are done - for (int i=0;istart, lines[0]->end); -#endif + float coef = obsAverage / probAverage; + + return coef; } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesPairs"); + errorOut(e, "Pintail", "getCoef"); exit(1); } } - /**************************************************************************************************/ -void Pintail::createProcessesObserved() { +void Pintail::createProcessesSpots() { 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) { @@ -486,7 +761,27 @@ void Pintail::createProcessesObserved() { processIDS.push_back(pid); process++; }else if (pid == 0){ - calcObserved(lines[process]->start, lines[process]->end); + + 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(); + + 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; + } + mothurOut("Done finding window breaks 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); } } @@ -497,25 +792,43 @@ void Pintail::createProcessesObserved() { wait(&temp); } + windowsForeachQuery = win; #else - calcObserved(lines[0]->start, lines[0]->end); + 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", "createProcessesObserved"); + errorOut(e, "Pintail", "createProcessesSpots"); exit(1); } } -//*************************************************************************************************************** -void Pintail::createProcessesExpected() { +/**************************************************************************************************/ + +void Pintail::createProcesses() { try { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) 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(); @@ -524,7 +837,33 @@ void Pintail::createProcessesExpected() { processIDS.push_back(pid); process++; }else if (pid == 0){ - calcExpected(lines[process]->start, lines[process]->end); + + 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); } } @@ -535,25 +874,77 @@ void Pintail::createProcessesExpected() { wait(&temp); } + obsDistance = obs; + expectedDistance = exp; + DE = de; + deviation = dev; + #else - calcExpected(lines[0]->start, lines[0]->end); + 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(); + + + + 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(); + + + + 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; + } + 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]); + 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 } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesExpected"); + errorOut(e, "Pintail", "createProcesses"); exit(1); } } + /**************************************************************************************************/ -void Pintail::createProcessesDE() { +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(); @@ -562,7 +953,20 @@ void Pintail::createProcessesDE() { processIDS.push_back(pid); process++; }else if (pid == 0){ - calcDE(lines[process]->start, lines[process]->end); + + 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()); + } + + 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; + } + exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } } @@ -573,17 +977,18 @@ void Pintail::createProcessesDE() { wait(&temp); } + quantiles = quan; #else - calcDE(lines[0]->start, lines[0]->end); - + quantiles = getQuantiles(0, templateSeqs.size()); #endif } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesDE"); + errorOut(e, "Pintail", "createProcessesQuan"); exit(1); } } + //***************************************************************************************************************