X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=decalc.cpp;h=98545f8e24a72dca332bf7571fc0e5e9e6936988;hp=b2e8be61f8df0aff52753058532a3ef42c6b6f0f;hb=cf9987b67aa49777a4c91c2d21f96e58bf17aa82;hpb=4f8291bc0a482c68dc739ef5aa7a311fc74f5f43 diff --git a/decalc.cpp b/decalc.cpp index b2e8be6..98545f8 100644 --- a/decalc.cpp +++ b/decalc.cpp @@ -8,11 +8,16 @@ */ #include "decalc.h" +#include "chimera.h" +#include "dist.h" +#include "eachgapdist.h" +#include "ignoregaps.h" +#include "eachgapdist.h" //*************************************************************************************************************** -void DeCalculator::setMask(string m) { +void DeCalculator::setMask(string ms) { try { - seqMask = m; + seqMask = ms; int count = 0; maskMap.clear(); @@ -35,7 +40,7 @@ void DeCalculator::setMask(string m) { } } catch(exception& e) { - errorOut(e, "DeCalculator", "setMask"); + m->errorOut(e, "DeCalculator", "setMask"); exit(1); } } @@ -57,7 +62,7 @@ void DeCalculator::runMask(Sequence* seq) { seq->setUnaligned(tempQuery); } catch(exception& e) { - errorOut(e, "DeCalculator", "runMask"); + m->errorOut(e, "DeCalculator", "runMask"); exit(1); } } @@ -85,7 +90,7 @@ void DeCalculator::trimSeqs(Sequence* query, Sequence* subject, map& t } catch(exception& e) { - errorOut(e, "DeCalculator", "trimSeqs"); + m->errorOut(e, "DeCalculator", "trimSeqs"); exit(1); } } @@ -101,7 +106,7 @@ vector DeCalculator::findWindows(Sequence* query, int front, int back, int 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(); + m->mothurOut("You have selected too large a window size for sequence " + query->getName() + ". I will choose an appropriate window size."); m->mothurOutEndLine(); size = (cutoff / 4); } @@ -145,7 +150,7 @@ vector DeCalculator::findWindows(Sequence* query, int front, int back, int */ //this follows wigeon, but we may want to consider that it chops off the end values if the sequence cannot be evenly divided into steps - for (int m = front; m < (back - size) ; m+=increment) { win.push_back(m); } + for (int i = front; i < (back - size) ; i+=increment) { win.push_back(i); } @@ -153,7 +158,7 @@ vector DeCalculator::findWindows(Sequence* query, int front, int back, int } catch(exception& e) { - errorOut(e, "DeCalculator", "findWindows"); + m->errorOut(e, "DeCalculator", "findWindows"); exit(1); } } @@ -164,10 +169,10 @@ vector DeCalculator::calcObserved(Sequence* query, Sequence* subject, vec vector temp; //int gaps = 0; - for (int m = 0; m < window.size(); m++) { + for (int i = 0; i < window.size(); i++) { - string seqFrag = query->getAligned().substr(window[m], size); - string seqFragsub = subject->getAligned().substr(window[m], size); + string seqFrag = query->getAligned().substr(window[i], size); + string seqFragsub = subject->getAligned().substr(window[i], size); int diff = 0; for (int b = 0; b < seqFrag.length(); b++) { @@ -195,7 +200,7 @@ vector DeCalculator::calcObserved(Sequence* query, Sequence* subject, vec return temp; } catch(exception& e) { - errorOut(e, "DeCalculator", "calcObserved"); + m->errorOut(e, "DeCalculator", "calcObserved"); exit(1); } } @@ -227,7 +232,7 @@ float DeCalculator::calcDist(Sequence* query, Sequence* subject, int front, int return dist; } catch(exception& e) { - errorOut(e, "DeCalculator", "calcDist"); + m->errorOut(e, "DeCalculator", "calcDist"); exit(1); } } @@ -239,9 +244,9 @@ vector DeCalculator::calcExpected(vector qav, float coef) { //for each window vector queryExpected; - for (int m = 0; m < qav.size(); m++) { + for (int j = 0; j < qav.size(); j++) { - float expected = qav[m] * coef; + float expected = qav[j] * coef; queryExpected.push_back(expected); } @@ -250,7 +255,7 @@ vector DeCalculator::calcExpected(vector qav, float coef) { } catch(exception& e) { - errorOut(e, "DeCalculator", "calcExpected"); + m->errorOut(e, "DeCalculator", "calcExpected"); exit(1); } } @@ -259,12 +264,12 @@ float DeCalculator::calcDE(vector obs, vector exp) { try { //for each window - float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2 + float sum = 0.0; //sum = sum from 1 to i of (oi-ei)^2 int numZeros = 0; - for (int m = 0; m < obs.size(); m++) { + for (int j = 0; j < obs.size(); j++) { - //if (obs[m] != 0.0) { - sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); + //if (obs[j] != 0.0) { + sum += ((obs[j] - exp[j]) * (obs[j] - exp[j])); //}else { numZeros++; } } @@ -274,7 +279,7 @@ float DeCalculator::calcDE(vector obs, vector exp) { return de; } catch(exception& e) { - errorOut(e, "DeCalculator", "calcDE"); + m->errorOut(e, "DeCalculator", "calcDE"); exit(1); } } @@ -285,10 +290,12 @@ vector DeCalculator::calcFreq(vector seqs, string filename) { try { vector prob; - string freqfile = getRootName(filename) + "freq"; + string freqfile = m->getRootName(filename) + "freq"; ofstream outFreq; - openOutputFile(freqfile, outFreq); + m->openOutputFile(freqfile, outFreq); + + outFreq << "#" << m->getVersion() << endl; string length = toString(seqs.size()); //if there are 5000 seqs in the template then set precision to 3 int precision = length.length() - 1; @@ -316,7 +323,7 @@ vector DeCalculator::calcFreq(vector seqs, string filename) { //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 j = 0; j < freq.size(); j++) { if (freq[j] > highest) { highest = freq[j]; } } float highFreq = highest / (float) (seqs.size()); @@ -340,7 +347,7 @@ vector DeCalculator::calcFreq(vector seqs, string filename) { } catch(exception& e) { - errorOut(e, "DeCalculator", "calcFreq"); + m->errorOut(e, "DeCalculator", "calcFreq"); exit(1); } } @@ -350,13 +357,13 @@ vector DeCalculator::findQav(vector window, int size, vector vector averages; //for each window find average - for (int m = 0; m < window.size(); m++) { + for (int i = 0; i < window.size(); i++) { float average = 0.0; //while you are in the window for this sequence int count = 0; - for (int j = window[m]; j < (window[m]+size); j++) { + for (int j = window[i]; j < (window[i]+size); j++) { average += probabilityProfile[j]; count++; } @@ -370,26 +377,23 @@ vector DeCalculator::findQav(vector window, int size, vector return averages; } catch(exception& e) { - errorOut(e, "DeCalculator", "findQav"); + m->errorOut(e, "DeCalculator", "findQav"); exit(1); } } //*************************************************************************************************************** //seqs have already been masked -vector< vector > DeCalculator::getQuantiles(vector seqs, vector windowSizesTemplate, int window, vector probProfile, int increment, int start, int end) { +vector< vector > DeCalculator::getQuantiles(vector seqs, vector windowSizesTemplate, int window, vector probProfile, int increment, int start, int end) { try { - vector< vector > quan; + vector< vector > quan; //percentage of mismatched pairs 1 to 100 quan.resize(100); -//ofstream o; -//string out = "getQuantiles.out"; -//openOutputFile(out, o); - + //for each sequence for(int i = start; i < end; i++){ - mothurOut("Processing sequence " + toString(i)); mothurOutEndLine(); + m->mothurOut("Processing sequence " + toString(i)); m->mothurOutEndLine(); Sequence* query = new Sequence(seqs[i]->getName(), seqs[i]->getAligned()); //compare to every other sequence in template @@ -397,6 +401,8 @@ vector< vector > DeCalculator::getQuantiles(vector seqs, Sequence* subject = new Sequence(seqs[j]->getName(), seqs[j]->getAligned()); + if (m->control_pressed) { delete query; delete subject; return quan; } + map trim; map::iterator it; @@ -421,25 +427,25 @@ vector< vector > DeCalculator::getQuantiles(vector seqs, float de = calcDE(obsi, exp); float dist = calcDist(query, subject, front, back); - //o << i << '\t' << j << '\t' << dist << '\t' << de << endl; + //cout << i << '\t' << j << '\t' << dist << '\t' << de << endl; dist = ceil(dist); - quanMember newScore(de, i, j); - - //dist-1 because vector indexes start at 0. - quan[dist-1].push_back(newScore); + //quanMember newScore(de, i, j); + quan[dist].push_back(de); + delete subject; } delete query; } + return quan; } catch(exception& e) { - errorOut(e, "DeCalculator", "getQuantiles"); + m->errorOut(e, "DeCalculator", "getQuantiles"); exit(1); } } @@ -450,27 +456,27 @@ inline bool compareQuanMembers(quanMember left, quanMember right){ } //*************************************************************************************************************** //this was going to be used by pintail to increase the sensitivity of the chimera detection, but it wasn't quite right. may want to revisit in the future... -void DeCalculator::removeObviousOutliers(vector< vector >& quantiles, int num) { +void DeCalculator::removeObviousOutliers(vector< vector >& quantiles, int num) { try { for (int i = 0; i < quantiles.size(); i++) { - - //find mean of this quantile score - sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers); - float high = quantiles[i][int(quantiles[i].size() * 0.99)].score; - float low = quantiles[i][int(quantiles[i].size() * 0.01)].score; - - vector temp; + //find mean of this quantile score + sort(quantiles[i].begin(), quantiles[i].end()); - //look at each value in quantiles to see if it is an outlier - for (int j = 0; j < quantiles[i].size(); j++) { - //is this score between 1 and 99% - if ((quantiles[i][j].score > low) && (quantiles[i][j].score < high)) { - temp.push_back(quantiles[i][j]); + vector temp; + if (quantiles[i].size() != 0) { + float high = quantiles[i][int(quantiles[i].size() * 0.99)]; + float low = quantiles[i][int(quantiles[i].size() * 0.01)]; + + //look at each value in quantiles to see if it is an outlier + for (int j = 0; j < quantiles[i].size(); j++) { + //is this score between 1 and 99% + if ((quantiles[i][j] > low) && (quantiles[i][j] < high)) { + temp.push_back(quantiles[i][j]); + } } } - quantiles[i] = temp; } @@ -481,7 +487,7 @@ void DeCalculator::removeObviousOutliers(vector< vector >& quantiles //while you still have guys to eliminate while (contributions.size() > 0) { - mothurOut("Removing scores contributed by sequence " + toString(largestContrib) + " in your template file."); mothurOutEndLine(); + m->mothurOut("Removing scores contributed by sequence " + toString(largestContrib) + " in your template file."); m->mothurOutEndLine(); //remove from quantiles all scores that were made using this largestContrib for (int i = 0; i < quantiles.size(); i++) { @@ -512,7 +518,7 @@ cout << "high = " << high << endl; for(int i = 0; i < marked.size(); i++) { if (marked[i] > high) { - mothurOut("Removing scores contributed by sequence " + toString(marked[i]) + " in your template file."); mothurOutEndLine(); + m->mothurOut("Removing scores contributed by sequence " + toString(marked[i]) + " in your template file."); m->mothurOutEndLine(); for (int i = 0; i < quantiles.size(); i++) { removeContrib(marked[i], quantiles[i]); } @@ -556,7 +562,7 @@ cout << "high = " << high << endl; } catch(exception& e) { - errorOut(e, "DeCalculator", "removeObviousOutliers"); + m->errorOut(e, "DeCalculator", "removeObviousOutliers"); exit(1); } } @@ -588,12 +594,12 @@ cout << largest->second << '\t' << largest->first->score << '\t' << largest->fir } catch(exception& e) { - errorOut(e, "DeCalculator", "sortContrib"); + m->errorOut(e, "DeCalculator", "sortContrib"); exit(1); } } -//*************************************************************************************************************** +*************************************************************************************************************** //used by removeObviousOutliers which was attempt to increase sensitivity of chimera detection...not currently used... int DeCalculator::findLargestContrib(vector seen) { try{ @@ -614,11 +620,11 @@ int DeCalculator::findLargestContrib(vector seen) { } catch(exception& e) { - errorOut(e, "DeCalculator", "findLargestContribs"); + m->errorOut(e, "DeCalculator", "findLargestContribs"); exit(1); } } -//*************************************************************************************************************** +*************************************************************************************************************** void DeCalculator::removeContrib(int bad, vector& quan) { try{ @@ -633,7 +639,7 @@ void DeCalculator::removeContrib(int bad, vector& quan) { } catch(exception& e) { - errorOut(e, "DeCalculator", "removeContrib"); + m->errorOut(e, "DeCalculator", "removeContrib"); exit(1); } } @@ -651,7 +657,7 @@ float DeCalculator::findAverage(vector myVector) { } catch(exception& e) { - errorOut(e, "DeCalculator", "findAverage"); + m->errorOut(e, "DeCalculator", "findAverage"); exit(1); } } @@ -671,10 +677,328 @@ float DeCalculator::getCoef(vector obs, vector qav) { return coef; } catch(exception& e) { - errorOut(e, "DeCalculator", "getCoef"); + m->errorOut(e, "DeCalculator", "getCoef"); + exit(1); + } +} +//*************************************************************************************************************** +//gets closest matches to each end, since chimeras will most likely have different parents on each end +vector DeCalculator::findClosest(Sequence querySeq, vector& thisTemplate, vector& thisFilteredTemplate, int numWanted, int minSim) { + try { + //indexes.clear(); + + vector seqsMatches; + + vector distsLeft; + vector distsRight; + + Dist* distcalculator = new eachGapDist(); + + string queryUnAligned = querySeq.getUnaligned(); + int numBases = int(queryUnAligned.length() * 0.33); + + string leftQuery = ""; //first 1/3 of the sequence + string rightQuery = ""; //last 1/3 of the sequence + string queryAligned = querySeq.getAligned(); + + //left side + bool foundFirstBase = false; + int baseCount = 0; + int leftSpot = 0; + int firstBaseSpot = 0; + for (int i = 0; i < queryAligned.length(); i++) { + //if you are a base + if (isalpha(queryAligned[i])) { + baseCount++; + if (!foundFirstBase) { foundFirstBase = true; firstBaseSpot = i; } + } + + //eliminate opening .'s + if (foundFirstBase) { leftQuery += queryAligned[i]; } + //if you have 1/3 + if (baseCount >= numBases) { leftSpot = i; break; } //first 1/3 + } + + //right side - count through another 1/3, so you are at last third + baseCount = 0; + int rightSpot = 0; + for (int i = leftSpot; i < queryAligned.length(); i++) { + //if you are a base + if (isalpha(queryAligned[i])) { baseCount++; } + //if you have 1/3 + if (baseCount > numBases + 1) { rightSpot = i; break; } //last 1/3 + } + + //trim end + //find last position in query that is a non gap character + int lastBaseSpot = queryAligned.length()-1; + for (int j = queryAligned.length()-1; j >= 0; j--) { + if (isalpha(queryAligned[j])) { + lastBaseSpot = j; + break; + } + } + rightQuery = queryAligned.substr(rightSpot, (lastBaseSpot-rightSpot+1)); //sequence from pos spot to end + + Sequence queryLeft(querySeq.getName(), leftQuery); + Sequence queryRight(querySeq.getName(), rightQuery); + +//cout << querySeq->getName() << '\t' << leftSpot << '\t' << rightSpot << '\t' << firstBaseSpot << '\t' << lastBaseSpot << endl; +//cout << queryUnAligned.length() << '\t' << queryLeft.getUnaligned().length() << '\t' << queryRight.getUnaligned().length() << endl; + for(int j = 0; j < thisFilteredTemplate.size(); j++){ + + string dbAligned = thisFilteredTemplate[j]->getAligned(); + string leftDB = dbAligned.substr(firstBaseSpot, (leftSpot-firstBaseSpot+1)); //first 1/3 of the sequence + string rightDB = dbAligned.substr(rightSpot, (lastBaseSpot-rightSpot+1)); //last 1/3 of the sequence + + Sequence dbLeft(thisFilteredTemplate[j]->getName(), leftDB); + Sequence dbRight(thisFilteredTemplate[j]->getName(), rightDB); + + distcalculator->calcDist(queryLeft, dbLeft); + float distLeft = distcalculator->getDist(); + + distcalculator->calcDist(queryRight, dbRight); + float distRight = distcalculator->getDist(); + + SeqDist subjectLeft; + subjectLeft.seq = NULL; + subjectLeft.dist = distLeft; + subjectLeft.index = j; + + distsLeft.push_back(subjectLeft); + + SeqDist subjectRight; + subjectRight.seq = NULL; + subjectRight.dist = distRight; + subjectRight.index = j; + + distsRight.push_back(subjectRight); + + } + + delete distcalculator; + + //sort by smallest distance + sort(distsRight.begin(), distsRight.end(), compareSeqDist); + sort(distsLeft.begin(), distsLeft.end(), compareSeqDist); + + + //merge results + map seen; + map::iterator it; + + vector dists; + float lastRight = distsRight[0].dist; + float lastLeft = distsLeft[0].dist; + + float maxDist = 1.0 - (minSim / 100.0); + + for (int i = 0; i < numWanted+1; i++) { + if (m->control_pressed) { return seqsMatches; } + + //add left if you havent already + it = seen.find(thisTemplate[distsLeft[i].index]->getName()); + if (it == seen.end() && distsLeft[i].dist <= maxDist) { + dists.push_back(distsLeft[i]); + seen[thisTemplate[distsLeft[i].index]->getName()] = thisTemplate[distsLeft[i].index]->getName(); + lastLeft = distsLeft[i].dist; +// cout << "loop-left\t" << db[distsLeft[i].index]->getName() << '\t' << distsLeft[i].dist << endl; + } + + //add right if you havent already + it = seen.find(thisTemplate[distsRight[i].index]->getName()); + if (it == seen.end() && distsRight[i].dist <= maxDist) { + dists.push_back(distsRight[i]); + seen[thisTemplate[distsRight[i].index]->getName()] = thisTemplate[distsRight[i].index]->getName(); + lastRight = distsRight[i].dist; +// cout << "loop-right\t" << db[distsRight[i].index]->getName() << '\t' << distsRight[i].dist << endl; + } + + if (i == numWanted) { break; } + + } + + //are we still above the minimum similarity cutoff + if ((lastLeft >= minSim) || (lastRight >= minSim)) { + //add in ties from left + int i = numWanted; + while (i < distsLeft.size()) { + if (distsLeft[i].dist == lastLeft) { dists.push_back(distsLeft[i]); } + else { break; } + i++; + } + + //add in ties from right + i = numWanted; + while (i < distsRight.size()) { + if (distsRight[i].dist == lastRight) { dists.push_back(distsRight[i]); } + else { break; } + i++; + } + } + + //cout << numWanted << endl; + for (int i = 0; i < dists.size(); i++) { +// cout << db[dists[i].index]->getName() << '\t' << dists[i].dist << endl; + + if ((thisTemplate[dists[i].index]->getName() != querySeq.getName()) && (((1.0-dists[i].dist)*100) >= minSim)) { + Sequence temp(thisTemplate[dists[i].index]->getName(), thisTemplate[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. + //cout << querySeq->getName() << '\t' << thisTemplate[dists[i].index]->getName() << '\t' << dists[i].dist << endl; + seqsMatches.push_back(temp); + } + + } + + return seqsMatches; + } + catch(exception& e) { + m->errorOut(e, "DeCalculator", "findClosest"); exit(1); } } //*************************************************************************************************************** +Sequence* DeCalculator::findClosest(Sequence* querySeq, vector db) { + try { + + Sequence* seqsMatch; + + Dist* distcalculator = new eachGapDist(); + int index = 0; + int smallest = 1000000; + + for(int j = 0; j < db.size(); j++){ + + distcalculator->calcDist(*querySeq, *db[j]); + float dist = distcalculator->getDist(); + + if (dist < smallest) { + smallest = dist; + index = j; + } + } + + delete distcalculator; + + seqsMatch = new Sequence(db[index]->getName(), db[index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. + + return seqsMatch; + } + catch(exception& e) { + m->errorOut(e, "DeCalculator", "findClosest"); + exit(1); + } +} +/***************************************************************************************************************/ +map DeCalculator::trimSeqs(Sequence& query, vector& topMatches) { + try { + + int frontPos = 0; //should contain first position in all seqs that is not a gap character + int rearPos = query.getAligned().length(); + + //********find first position in topMatches that is a non gap character***********// + //find first position all query seqs that is a non gap character + for (int i = 0; i < topMatches.size(); i++) { + + string aligned = topMatches[i].getAligned(); + int pos = 0; + + //find first spot in this seq + for (int j = 0; j < aligned.length(); j++) { + if (isalpha(aligned[j])) { + pos = j; + break; + } + } + + //save this spot if it is the farthest + if (pos > frontPos) { frontPos = pos; } + } + + + string aligned = query.getAligned(); + int pos = 0; + + //find first position in query that is a non gap character + for (int j = 0; j < aligned.length(); j++) { + if (isalpha(aligned[j])) { + pos = j; + break; + } + } + + //save this spot if it is the farthest + if (pos > frontPos) { frontPos = pos; } + + + //********find last position in topMatches that is a non gap character***********// + for (int i = 0; i < topMatches.size(); i++) { + + string aligned = topMatches[i].getAligned(); + int pos = aligned.length(); + + //find first spot in this seq + for (int j = aligned.length()-1; j >= 0; j--) { + if (isalpha(aligned[j])) { + pos = j; + break; + } + } + + //save this spot if it is the farthest + if (pos < rearPos) { rearPos = pos; } + } + + + aligned = query.getAligned(); + pos = aligned.length(); + + //find last position in query that is a non gap character + for (int j = aligned.length()-1; j >= 0; j--) { + if (isalpha(aligned[j])) { + pos = j; + break; + } + } + + //save this spot if it is the farthest + if (pos < rearPos) { rearPos = pos; } + + map trimmedPos; + //check to make sure that is not whole seq + if ((rearPos - frontPos - 1) <= 0) { + query.setAligned(""); + //trim topMatches + for (int i = 0; i < topMatches.size(); i++) { + topMatches[i].setAligned(""); + } + + }else { + + //trim query + string newAligned = query.getAligned(); + newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1)); + query.setAligned(newAligned); + + //trim topMatches + for (int i = 0; i < topMatches.size(); i++) { + newAligned = topMatches[i].getAligned(); + newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1)); + topMatches[i].setAligned(newAligned); + } + + for (int i = 0; i < newAligned.length(); i++) { + trimmedPos[i] = i+frontPos; + } + } + return trimmedPos; + } + catch(exception& e) { + m->errorOut(e, "DeCalculator", "trimSequences"); + exit(1); + } + +} +//*************************************************************************************************************** +