X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=decalc.cpp;h=98545f8e24a72dca332bf7571fc0e5e9e6936988;hp=67258714aa3b45cbdeaee2cf8d43576f2de4f446;hb=cf9987b67aa49777a4c91c2d21f96e58bf17aa82;hpb=86c838c428a9e7d26f902f5492738241fa72c4e7 diff --git a/decalc.cpp b/decalc.cpp index 6725871..98545f8 100644 --- a/decalc.cpp +++ b/decalc.cpp @@ -12,7 +12,7 @@ #include "dist.h" #include "eachgapdist.h" #include "ignoregaps.h" -#include "eachgapdistignorens.h" +#include "eachgapdist.h" //*************************************************************************************************************** void DeCalculator::setMask(string ms) { @@ -599,7 +599,7 @@ cout << largest->second << '\t' << largest->first->score << '\t' << largest->fir } } -//*************************************************************************************************************** +*************************************************************************************************************** //used by removeObviousOutliers which was attempt to increase sensitivity of chimera detection...not currently used... int DeCalculator::findLargestContrib(vector seen) { try{ @@ -624,7 +624,7 @@ int DeCalculator::findLargestContrib(vector seen) { exit(1); } } -//*************************************************************************************************************** +*************************************************************************************************************** void DeCalculator::removeContrib(int bad, vector& quan) { try{ @@ -683,23 +683,23 @@ float DeCalculator::getCoef(vector obs, vector qav) { } //*************************************************************************************************************** //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) { +vector DeCalculator::findClosest(Sequence querySeq, vector& thisTemplate, vector& thisFilteredTemplate, int numWanted, int minSim) { try { - indexes.clear(); + //indexes.clear(); - vector seqsMatches; + vector seqsMatches; vector distsLeft; vector distsRight; - Dist* distcalculator = new eachGapDistIgnoreNs(); + Dist* distcalculator = new eachGapDist(); - string queryUnAligned = querySeq->getUnaligned(); + 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(); + string queryAligned = querySeq.getAligned(); //left side bool foundFirstBase = false; @@ -726,7 +726,7 @@ vector DeCalculator::findClosest(Sequence* querySeq, vector= numBases) { rightSpot = i; break; } //last 1/3 + if (baseCount > numBases + 1) { rightSpot = i; break; } //last 1/3 } //trim end @@ -738,27 +738,28 @@ vector DeCalculator::findClosest(Sequence* querySeq, vectorgetName(), 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++){ + for(int j = 0; j < thisFilteredTemplate.size(); j++){ - string dbAligned = db[j]->getAligned(); + 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)); //last 1/3 of the sequence - - Sequence dbLeft(db[j]->getName(), leftDB); - Sequence dbRight(db[j]->getName(), rightDB); + 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; @@ -780,14 +781,7 @@ vector DeCalculator::findClosest(Sequence* querySeq, vectorgetName() << '\t' << distsLeft[i].dist << endl; -// } -// for(int i=0;i<15;i++){ -// cout << "right\t" << db[distsLeft[i].index]->getName() << '\t' << distsRight[i].dist << endl; -// } - + //merge results map seen; @@ -796,83 +790,63 @@ vector DeCalculator::findClosest(Sequence* querySeq, vector dists; float lastRight = distsRight[0].dist; float lastLeft = distsLeft[0].dist; - int lasti = 0; - for (int i = 0; i < distsLeft.size(); i++) { + + 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(db[distsLeft[i].index]->getName()); - if (it == seen.end()) { + it = seen.find(thisTemplate[distsLeft[i].index]->getName()); + if (it == seen.end() && distsLeft[i].dist <= maxDist) { dists.push_back(distsLeft[i]); - seen[db[distsLeft[i].index]->getName()] = db[distsLeft[i].index]->getName(); + 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(db[distsRight[i].index]->getName()); - if (it == seen.end()) { + it = seen.find(thisTemplate[distsRight[i].index]->getName()); + if (it == seen.end() && distsRight[i].dist <= maxDist) { dists.push_back(distsRight[i]); - seen[db[distsRight[i].index]->getName()] = db[distsRight[i].index]->getName(); + 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 (dists.size() > numWanted) { lasti = i; break; } //you have enough results + if (i == numWanted) { break; } + } -// cout << "lastLeft\t" << lastLeft << endl; - - //add in sequences with same distance as last sequence added - lasti++; - int i = lasti; - while (i < distsLeft.size()) { - if (distsLeft[i].dist == lastLeft) { - it = seen.find(db[distsLeft[i].index]->getName()); - - if (it == seen.end()) { -// cout << "newLoop-left\t" << db[distsLeft[i].index]->getName() << '\t' << distsLeft[i].dist << endl; - dists.push_back(distsLeft[i]); - seen[db[distsRight[i].index]->getName()] = db[distsLeft[i].index]->getName(); -// numWanted++; - } + //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++; } - else { break; } - i++; - } - -// cout << "lastRight\t" << lastRight << endl; - //add in sequences with same distance as last sequence added - i = lasti; - while (i < distsRight.size()) { - if (distsRight[i].dist == lastRight) { - it = seen.find(db[distsRight[i].index]->getName()); - - if (it == seen.end()) { -// cout << "newLoop-right\t" << db[distsRight[i].index]->getName() << '\t' << distsRight[i].dist << endl; - dists.push_back(distsRight[i]); - seen[db[distsRight[i].index]->getName()] = db[distsRight[i].index]->getName(); -// numWanted++; - } + + //add in ties from right + i = numWanted; + while (i < distsRight.size()) { + if (distsRight[i].dist == lastRight) { dists.push_back(distsRight[i]); } + else { break; } + i++; } - else { break; } - i++; } - - numWanted = seen.size(); - 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 << numWanted << endl; + for (int i = 0; i < dists.size(); i++) { // cout << db[dists[i].index]->getName() << '\t' << dists[i].dist << endl; - if (db[dists[i].index]->getName() != querySeq->getName()) { - Sequence* temp = new Sequence(db[dists[i].index]->getName(), db[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. - + + 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); - indexes.push_back(dists[i].index); } + } return seqsMatches; @@ -888,7 +862,7 @@ Sequence* DeCalculator::findClosest(Sequence* querySeq, vector db) { Sequence* seqsMatch; - Dist* distcalculator = new eachGapDistIgnoreNs(); + Dist* distcalculator = new eachGapDist(); int index = 0; int smallest = 1000000; @@ -915,17 +889,17 @@ Sequence* DeCalculator::findClosest(Sequence* querySeq, vector db) { } } /***************************************************************************************************************/ -map DeCalculator::trimSeqs(Sequence* query, vector topMatches) { +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(); + 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(); + string aligned = topMatches[i].getAligned(); int pos = 0; //find first spot in this seq @@ -941,7 +915,7 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch } - string aligned = query->getAligned(); + string aligned = query.getAligned(); int pos = 0; //find first position in query that is a non gap character @@ -959,7 +933,7 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch //********find last position in topMatches that is a non gap character***********// for (int i = 0; i < topMatches.size(); i++) { - string aligned = topMatches[i]->getAligned(); + string aligned = topMatches[i].getAligned(); int pos = aligned.length(); //find first spot in this seq @@ -975,7 +949,7 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch } - aligned = query->getAligned(); + aligned = query.getAligned(); pos = aligned.length(); //find last position in query that is a non gap character @@ -992,24 +966,24 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch map trimmedPos; //check to make sure that is not whole seq if ((rearPos - frontPos - 1) <= 0) { - query->setAligned(""); + query.setAligned(""); //trim topMatches for (int i = 0; i < topMatches.size(); i++) { - topMatches[i]->setAligned(""); + topMatches[i].setAligned(""); } }else { //trim query - string newAligned = query->getAligned(); + string newAligned = query.getAligned(); newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1)); - query->setAligned(newAligned); + query.setAligned(newAligned); //trim topMatches for (int i = 0; i < topMatches.size(); i++) { - newAligned = topMatches[i]->getAligned(); + newAligned = topMatches[i].getAligned(); newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1)); - topMatches[i]->setAligned(newAligned); + topMatches[i].setAligned(newAligned); } for (int i = 0; i < newAligned.length(); i++) {