distancedb.cpp

   1 /*
   2  *  distancedb.cpp
   3  *
   4  *
   5  *  Created by Pat Schloss on 12/29/08.
   6  *  Copyright 2008 Patrick D. Schloss. All rights reserved.
   7  *
   8  */
   9
  10
  11 #include "database.hpp"
  12 #include "sequence.hpp"
  13 #include "distancedb.hpp"
  14 #include "eachgapdist.h"
  15
  16 /**************************************************************************************************/
  17 DistanceDB::DistanceDB() {
  18         try {
  19                 templateAligned = true;
  20                 templateSeqsLength = 0;
  21                 distCalculator = new eachGapDist();
  22         }
  23         catch(exception& e) {
  24                 m->errorOut(e, "DistanceDB", "DistanceDB");
  25                 exit(1);
  26         }
  27 }
  28 /**************************************************************************************************/
  29 void DistanceDB::addSequence(Sequence seq) {
  30         try {
  31                 //are the template sequences aligned
  32                 if (!isAligned(seq.getAligned())) { templateAligned = false; m->mothurOut(seq.getName() + " is not aligned. Sequences must be aligned to use the distance method."); m->mothurOutEndLine(); }
  33
  34                 if (templateSeqsLength == 0) { templateSeqsLength = seq.getAligned().length(); }
  35
  36                 data.push_back(seq);
  37         }
  38         catch(exception& e) {
  39                 m->errorOut(e, "DistanceDB", "addSequence");
  40                 exit(1);
  41         }
  42 }
  43 /**************************************************************************************************/
  44 //returns indexes to top matches
  45 vector<int> DistanceDB::findClosestSequences(Sequence* query, int numWanted){
  46         try {
  47                 vector<int> topMatches;
  48                 bool templateSameLength = true;
  49                 string sequence = query->getAligned();
  50                 vector<seqDist> dists;
  51
  52                 searchScore = -1.0;
  53
  54                 if (numWanted > data.size()) { m->mothurOut("numwanted is larger than the number of template sequences, using "+ toString(data.size()) + "."); m->mothurOutEndLine(); numWanted = data.size(); }
  55
  56                 if (sequence.length() != templateSeqsLength) { templateSameLength = false; }
  57
  58                 if (templateSameLength && templateAligned) {
  59                         //calc distance from this sequence to every sequence in the template
  60                         for (int i = 0; i < data.size(); i++) {
  61                                 distCalculator->calcDist(*query, data[i]);
  62                                 float dist = distCalculator->getDist();
  63
  64                                 //save distance to each template sequence
  65                                 seqDist temp(-1, i, dist);
  66                                 dists.push_back(temp);
  67                         }
  68
  69                         sort(dists.begin(), dists.end(), compareSequenceDistance);  //sorts by distance lowest to highest
  70
  71                         //save distance of best match
  72                         searchScore = dists[0].dist;
  73
  74                         //fill topmatches with numwanted closest sequences indexes
  75                         for (int i = 0; i < numWanted; i++) {
  76                                 topMatches.push_back(dists[i].seq2);
  77                         }
  78
  79                 }else{
  80                         m->mothurOut("cannot find closest matches using distance method for " + query->getName() + " without aligned template sequences of the same length."); m->mothurOutEndLine();
  81                         exit(1);
  82                 }
  83
  84                 return topMatches;
  85         }
  86         catch(exception& e) {
  87                 m->errorOut(e, "DistanceDB", "findClosestSequence");
  88                 exit(1);
  89         }
  90 }
  91 /**************************************************************************************************/
  92 bool DistanceDB::isAligned(string seq){
  93         try {
  94                 bool aligned;
  95
  96                 int pos = seq.find_first_of(".-");
  97
  98                 if (pos != seq.npos) {
  99                         aligned = true;
 100                 }else { aligned = false; }
 101
 102
 103                 return aligned;
 104         }
 105         catch(exception& e) {
 106                 m->errorOut(e, "DistanceDB", "isAligned");
 107                 exit(1);
 108         }
 109 }
 110
 111 /**************************************************************************************************/