alignmentdb.cpp

   1 /*
   2  *  alignmentdb.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 11/4/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "alignmentdb.h"
  11 #include "kmerdb.hpp"
  12 #include "suffixdb.hpp"
  13 #include "blastdb.hpp"
  14
  15
  16 /**************************************************************************************************/
  17 AlignmentDB::AlignmentDB(string fastaFileName, string method, int kmerSize, float gapOpen, float gapExtend, float match, float misMatch){               //      This assumes that the template database is in fasta format, may
  18         try {                                                                                   //      need to alter this in the future?
  19                 m = MothurOut::getInstance();
  20                 longest = 0;
  21
  22                 ifstream fastaFile;
  23                 openInputFile(fastaFileName, fastaFile);
  24
  25                 m->mothurOutEndLine();
  26                 m->mothurOut("Reading in the " + fastaFileName + " template sequences...\t");   cout.flush();
  27
  28                 while (!fastaFile.eof()) {
  29                         Sequence temp(fastaFile);  gobble(fastaFile);
  30
  31                         if (temp.getName() != "") {
  32                                 templateSequences.push_back(temp);
  33                                 //save longest base
  34                                 if (temp.getUnaligned().length() > longest)  { longest = temp.getUnaligned().length()+1; }
  35                         }
  36                 }
  37
  38                 numSeqs = templateSequences.size();
  39
  40                 fastaFile.close();
  41                 //all of this is elsewhere already!
  42
  43                 m->mothurOut("DONE.");
  44                 m->mothurOutEndLine();  cout.flush();
  45
  46                 //in case you delete the seqs and then ask for them
  47                 emptySequence = Sequence();
  48                 emptySequence.setName("no_match");
  49                 emptySequence.setUnaligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
  50                 emptySequence.setAligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
  51
  52                 bool needToGenerate = true;
  53                 string kmerDBName;
  54                 if(method == "kmer")                    {
  55                         search = new KmerDB(fastaFileName, kmerSize);
  56
  57                         kmerDBName = fastaFileName.substr(0,fastaFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
  58                         ifstream kmerFileTest(kmerDBName.c_str());
  59
  60                         if(kmerFileTest){       needToGenerate = false;         }
  61                 }
  62                 else if(method == "suffix")             {       search = new SuffixDB(numSeqs);                                                         }
  63                 else if(method == "blast")              {       search = new BlastDB(gapOpen, gapExtend, match, misMatch);      }
  64                 else {
  65                         m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
  66                         m->mothurOutEndLine();
  67                         search = new KmerDB(fastaFileName, 8);
  68                 }
  69
  70                 if (needToGenerate) {
  71
  72                         //add sequences to search
  73                         for (int i = 0; i < templateSequences.size(); i++) {
  74                                 search->addSequence(templateSequences[i]);
  75                         }
  76                         search->generateDB();
  77
  78                 }else if ((method == "kmer") && (!needToGenerate)) {
  79                         ifstream kmerFileTest(kmerDBName.c_str());
  80                         search->readKmerDB(kmerFileTest);
  81                 }
  82
  83                 search->setNumSeqs(numSeqs);
  84         }
  85         catch(exception& e) {
  86                 m->errorOut(e, "AlignmentDB", "AlignmentDB");
  87                 exit(1);
  88         }
  89 }
  90 /**************************************************************************************************/
  91 AlignmentDB::~AlignmentDB() {  delete search;   }
  92 /**************************************************************************************************/
  93 Sequence AlignmentDB::findClosestSequence(Sequence* seq) {
  94         try{
  95
  96                 vector<int> spot = search->findClosestSequences(seq, 1);
  97
  98                 if (spot.size() != 0)   {               return templateSequences[spot[0]];      }
  99                 else                                    {               return emptySequence;                           }
 100
 101         }
 102         catch(exception& e) {
 103                 m->errorOut(e, "AlignmentDB", "findClosestSequence");
 104                 exit(1);
 105         }
 106 }
 107 /**************************************************************************************************/
 108
 109
 110
 111