5 * Created by Pat Schloss on 12/29/08.
6 * Copyright 2008 Patrick D. Schloss. All rights reserved.
11 #include "database.hpp"
12 #include "sequence.hpp"
13 #include "distancedb.hpp"
14 #include "onegapignore.h"
16 /**************************************************************************************************/
17 DistanceDB::DistanceDB() {
19 templateAligned = true;
20 templateSeqsLength = 0;
21 distCalculator = new oneGapIgnoreTermGapDist();
24 m->errorOut(e, "DistanceDB", "DistanceDB");
28 /**************************************************************************************************/
29 void DistanceDB::addSequence(Sequence seq) {
31 //are the template sequences aligned
32 if (!isAligned(seq.getAligned())) {
33 templateAligned = false;
34 m->mothurOut(seq.getName() + " is not aligned. Sequences must be aligned to use the distance method.");
35 m->mothurOutEndLine();
38 if (templateSeqsLength == 0) { templateSeqsLength = seq.getAligned().length(); }
43 m->errorOut(e, "DistanceDB", "addSequence");
47 /**************************************************************************************************/
48 //returns indexes to top matches
49 vector<int> DistanceDB::findClosestSequences(Sequence* query, int numWanted){
51 vector<int> topMatches;
52 bool templateSameLength = true;
53 string sequence = query->getAligned();
54 vector<seqDist> dists;
58 if (numWanted > data.size()){
59 m->mothurOut("numwanted is larger than the number of template sequences, using "+ toString(data.size()) + ".");
60 m->mothurOutEndLine();
61 numWanted = data.size();
64 if (sequence.length() != templateSeqsLength) { templateSameLength = false; }
66 if (templateSameLength && templateAligned) {
69 dists.resize(data.size());
71 //calc distance from this sequence to every sequence in the template
72 for (int i = 0; i < data.size(); i++) {
73 distCalculator->calcDist(*query, data[i]);
74 float dist = distCalculator->getDist();
76 //save distance to each template sequence
82 sort(dists.begin(), dists.end(), compareSequenceDistance); //sorts by distance lowest to highest
84 //save distance of best match
85 searchScore = dists[0].dist;
87 //fill topmatches with numwanted closest sequences indexes
88 for (int i = 0; i < numWanted; i++) {
89 topMatches.push_back(dists[i].seq2);
93 float smallDist = 100000;
94 for (int i = 0; i < data.size(); i++) {
95 distCalculator->calcDist(*query, data[i]);
96 float dist = distCalculator->getDist();
99 if (dist < smallDist) {
104 searchScore = smallDist;
105 topMatches.push_back(bestIndex);
109 m->mothurOut("cannot find closest matches using distance method for " + query->getName() + " without aligned template sequences of the same length.");
110 m->mothurOutEndLine();
116 catch(exception& e) {
117 m->errorOut(e, "DistanceDB", "findClosestSequence");
121 /**************************************************************************************************/
122 bool DistanceDB::isAligned(string seq){
126 int pos = seq.find_first_of(".-");
128 if (pos != seq.npos) {
130 }else { aligned = false; }
135 catch(exception& e) {
136 m->errorOut(e, "DistanceDB", "isAligned");
141 /**************************************************************************************************/