#include "dist.h"
#include "eachgapdist.h"
#include "ignoregaps.h"
-#include "eachgapdistignorens.h"
+#include "eachgapdist.h"
//***************************************************************************************************************
void DeCalculator::setMask(string ms) {
}
//***************************************************************************************************************
//gets closest matches to each end, since chimeras will most likely have different parents on each end
-vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*> db, int& numWanted, vector<int>& indexes) {
+vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*>& thisTemplate, vector<Sequence*>& thisFilteredTemplate, int numWanted, int minSim) {
try {
- indexes.clear();
+ //indexes.clear();
vector<Sequence*> seqsMatches;
vector<SeqDist> distsLeft;
vector<SeqDist> distsRight;
- Dist* distcalculator = new eachGapDistIgnoreNs();
+ Dist* distcalculator = new eachGapDist();
string queryUnAligned = querySeq->getUnaligned();
int numBases = int(queryUnAligned.length() * 0.33);
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);
+ Sequence dbLeft(thisFilteredTemplate[j]->getName(), leftDB);
+ Sequence dbRight(thisFilteredTemplate[j]->getName(), rightDB);
distcalculator->calcDist(queryLeft, dbLeft);
float distLeft = distcalculator->getDist();
//sort by smallest distance
sort(distsRight.begin(), distsRight.end(), compareSeqDist);
sort(distsLeft.begin(), distsLeft.end(), compareSeqDist);
-// cout << distsLeft.size() << '\t' << distsRight.size() << endl;
-// for(int i=0;i<15;i++){
-// cout << "left\t" << db[distsLeft[i].index]->getName() << '\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<string, string> seen;
vector<SeqDist> dists;
float lastRight = distsRight[0].dist;
float lastLeft = distsLeft[0].dist;
- int lasti = 0;
- for (int i = 0; i < distsLeft.size(); i++) {
+ //int lasti = 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());
+ it = seen.find(thisTemplate[distsLeft[i].index]->getName());
if (it == seen.end()) {
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());
+ it = seen.find(thisTemplate[distsRight[i].index]->getName());
if (it == seen.end()) {
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
}
-// cout << "lastLeft\t" << lastLeft << endl;
-
- //add in dups
- 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;
-
- 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;
- 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.
-
- seqsMatches.push_back(temp);
- indexes.push_back(dists[i].index);
+ if ((thisTemplate[dists[i].index]->getName() != querySeq->getName()) && (((1.0-dists[i].dist)*100) >= minSim)) {
+ Sequence* temp = new Sequence(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;
Sequence* seqsMatch;
- Dist* distcalculator = new eachGapDistIgnoreNs();
+ Dist* distcalculator = new eachGapDist();
int index = 0;
int smallest = 1000000;
projects / mothur.git / blobdiff