X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=decalc.cpp;h=4e80c65d29aa7c8d6a2354db95027c3365c7962d;hb=9946a1b4b50969d08ce059b248bdeecafbf989ac;hp=05897517d061540e6548b9bf2c51c9c97443248d;hpb=f7fbd74ffedcf62217109c22e828453eaefa1458;p=mothur.git diff --git a/decalc.cpp b/decalc.cpp index 0589751..4e80c65 100644 --- a/decalc.cpp +++ b/decalc.cpp @@ -8,25 +8,39 @@ */ #include "decalc.h" +#include "chimera.h" +#include "dist.h" +#include "eachgapdist.h" +#include "ignoregaps.h" + //*************************************************************************************************************** -void DeCalculator::setMask(string m) { +void DeCalculator::setMask(string ms) { try { - seqMask = m; + seqMask = ms; + int count = 0; + maskMap.clear(); if (seqMask.length() != 0) { //whereever there is a base in the mask, save that value is query and subject for (int i = 0; i < seqMask.length(); i++) { if (isalpha(seqMask[i])) { h.insert(i); + maskMap[count] = i; + count++; + } } }else { - for (int i = 0; i < alignLength; i++) { h.insert(i); } + for (int i = 0; i < alignLength; i++) { + h.insert(i); + maskMap[count] = i; + count++; + } } } catch(exception& e) { - errorOut(e, "DeCalculator", "setMask"); + m->errorOut(e, "DeCalculator", "setMask"); exit(1); } } @@ -48,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); } } @@ -76,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); } } @@ -92,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); } @@ -136,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); } @@ -144,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); } } @@ -155,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++) { @@ -186,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); } } @@ -218,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); } } @@ -230,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); } @@ -241,7 +255,7 @@ vector DeCalculator::calcExpected(vector qav, float coef) { } catch(exception& e) { - errorOut(e, "DeCalculator", "calcExpected"); + m->errorOut(e, "DeCalculator", "calcExpected"); exit(1); } } @@ -250,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++; } } @@ -265,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); } } @@ -307,7 +321,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()); @@ -331,7 +345,7 @@ vector DeCalculator::calcFreq(vector seqs, string filename) { } catch(exception& e) { - errorOut(e, "DeCalculator", "calcFreq"); + m->errorOut(e, "DeCalculator", "calcFreq"); exit(1); } } @@ -341,13 +355,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++; } @@ -361,7 +375,7 @@ vector DeCalculator::findQav(vector window, int size, vector return averages; } catch(exception& e) { - errorOut(e, "DeCalculator", "findQav"); + m->errorOut(e, "DeCalculator", "findQav"); exit(1); } } @@ -373,14 +387,11 @@ vector< vector > DeCalculator::getQuantiles(vector seqs, //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 @@ -388,6 +399,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; @@ -412,25 +425,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); - + quan[dist].push_back(newScore); + delete subject; } delete query; } + return quan; } catch(exception& e) { - errorOut(e, "DeCalculator", "getQuantiles"); + m->errorOut(e, "DeCalculator", "getQuantiles"); exit(1); } } @@ -445,23 +458,23 @@ void DeCalculator::removeObviousOutliers(vector< vector >& quantiles 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; - - //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]); + if (quantiles[i].size() != 0) { + float high = quantiles[i][int(quantiles[i].size() * 0.99)].score; + float low = quantiles[i][int(quantiles[i].size() * 0.01)].score; + + //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]); + } } } - quantiles[i] = temp; } @@ -472,7 +485,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++) { @@ -503,7 +516,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]); } @@ -547,7 +560,7 @@ cout << "high = " << high << endl; } catch(exception& e) { - errorOut(e, "DeCalculator", "removeObviousOutliers"); + m->errorOut(e, "DeCalculator", "removeObviousOutliers"); exit(1); } } @@ -579,7 +592,7 @@ cout << largest->second << '\t' << largest->first->score << '\t' << largest->fir } catch(exception& e) { - errorOut(e, "DeCalculator", "sortContrib"); + m->errorOut(e, "DeCalculator", "sortContrib"); exit(1); } } @@ -605,7 +618,7 @@ int DeCalculator::findLargestContrib(vector seen) { } catch(exception& e) { - errorOut(e, "DeCalculator", "findLargestContribs"); + m->errorOut(e, "DeCalculator", "findLargestContribs"); exit(1); } } @@ -624,7 +637,7 @@ void DeCalculator::removeContrib(int bad, vector& quan) { } catch(exception& e) { - errorOut(e, "DeCalculator", "removeContrib"); + m->errorOut(e, "DeCalculator", "removeContrib"); exit(1); } } @@ -642,7 +655,7 @@ float DeCalculator::findAverage(vector myVector) { } catch(exception& e) { - errorOut(e, "DeCalculator", "findAverage"); + m->errorOut(e, "DeCalculator", "findAverage"); exit(1); } } @@ -662,10 +675,307 @@ 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 db, int& numWanted, vector& indexes) { + 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) { 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)); //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 < db.size(); j++){ + + string dbAligned = db[j]->getAligned(); + string leftDB = dbAligned.substr(firstBaseSpot, (leftSpot-firstBaseSpot+1)); //first 1/3 of the sequence + string rightDB = dbAligned.substr(rightSpot, (lastBaseSpot-rightSpot)); //last 1/3 of the sequence + + Sequence dbLeft(db[j]->getName(), leftDB); + Sequence dbRight(db[j]->getName(), rightDB); + + distcalculator->calcDist(queryLeft, dbLeft); + float distLeft = distcalculator->getDist(); + + distcalculator->calcDist(queryRight, dbRight); + float distRight = distcalculator->getDist(); + + SeqDist subjectLeft; + subjectLeft.seq = db[j]; + subjectLeft.dist = distLeft; + subjectLeft.index = j; + + distsLeft.push_back(subjectLeft); + + SeqDist subjectRight; + subjectRight.seq = db[j]; + 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; + int lasti = 0; + for (int i = 0; i < distsLeft.size(); i++) { + //add left if you havent already + it = seen.find(distsLeft[i].seq->getName()); + if (it == seen.end()) { + dists.push_back(distsLeft[i]); + seen[distsLeft[i].seq->getName()] = distsLeft[i].seq->getName(); + lastLeft = distsLeft[i].dist; + } + + //add right if you havent already + it = seen.find(distsRight[i].seq->getName()); + if (it == seen.end()) { + dists.push_back(distsRight[i]); + seen[distsRight[i].seq->getName()] = distsRight[i].seq->getName(); + lastRight = distsRight[i].dist; + } + + if (dists.size() > numWanted) { lasti = i; break; } //you have enough results + } + + //add in dups + lasti++; + int i = lasti; + while (i < distsLeft.size()) { + if (distsLeft[i].dist == lastLeft) { dists.push_back(distsLeft[i]); numWanted++; } + else { break; } + i++; + } + + i = lasti; + while (i < distsRight.size()) { + if (distsRight[i].dist == lastRight) { dists.push_back(distsRight[i]); numWanted++; } + else { break; } + i++; + } + + if (numWanted > dists.size()) { m->mothurOut("numwanted is larger than the number of template sequences, adjusting numwanted."); m->mothurOutEndLine(); numWanted = dists.size(); } + +//cout << numWanted << endl; + for (int i = 0; i < numWanted; i++) { +//cout << dists[i].seq->getName() << '\t' << dists[i].dist << endl; + Sequence* temp = new Sequence(dists[i].seq->getName(), dists[i].seq->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. + seqsMatches.push_back(temp); + indexes.push_back(dists[i].index); + } + + 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; } + + //check to make sure that is not whole seq + if ((rearPos - frontPos - 1) <= 0) { m->mothurOut("Error, when I trim your sequences, the entire sequence is trimmed."); m->mothurOutEndLine(); exit(1); } +//cout << query->getName() << " front = " << frontPos << " rear = " << rearPos << endl; + //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); + } + + map trimmedPos; + + for (int i = 0; i < newAligned.length(); i++) { + trimmedPos[i] = i+frontPos; + } + + return trimmedPos; + } + catch(exception& e) { + m->errorOut(e, "DeCalculator", "trimSequences"); + exit(1); + } + +} +//*************************************************************************************************************** +