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 (m->control_pressed) {  templateSequences.clear(); break;  }
  32
  33                         if (temp.getName() != "") {
  34                                 templateSequences.push_back(temp);
  35                                 //save longest base
  36                                 if (temp.getUnaligned().length() > longest)  { longest = temp.getUnaligned().length()+1; }
  37                         }
  38                 }
  39
  40                 numSeqs = templateSequences.size();
  41
  42                 fastaFile.close();
  43                 //all of this is elsewhere already!
  44
  45                 m->mothurOut("DONE.");
  46                 m->mothurOutEndLine();  cout.flush();
  47
  48                 //in case you delete the seqs and then ask for them
  49                 emptySequence = Sequence();
  50                 emptySequence.setName("no_match");
  51                 emptySequence.setUnaligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
  52                 emptySequence.setAligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
  53
  54                 bool needToGenerate = true;
  55                 string kmerDBName;
  56                 if(method == "kmer")                    {
  57                         search = new KmerDB(fastaFileName, kmerSize);
  58
  59                         kmerDBName = fastaFileName.substr(0,fastaFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
  60                         ifstream kmerFileTest(kmerDBName.c_str());
  61
  62                         if(kmerFileTest){       needToGenerate = false;         }
  63                 }
  64                 else if(method == "suffix")             {       search = new SuffixDB(numSeqs);                                                         }
  65                 else if(method == "blast")              {       search = new BlastDB(gapOpen, gapExtend, match, misMatch);      }
  66                 else {
  67                         m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
  68                         m->mothurOutEndLine();
  69                         search = new KmerDB(fastaFileName, 8);
  70                 }
  71
  72                 if (!(m->control_pressed)) {
  73                         if (needToGenerate) {
  74                                 //add sequences to search
  75                                 for (int i = 0; i < templateSequences.size(); i++) {
  76                                         search->addSequence(templateSequences[i]);
  77                                 }
  78                                 search->generateDB();
  79
  80                         }else if ((method == "kmer") && (!needToGenerate)) {
  81                                 ifstream kmerFileTest(kmerDBName.c_str());
  82                                 search->readKmerDB(kmerFileTest);
  83                         }
  84
  85                         search->setNumSeqs(numSeqs);
  86                 }
  87         }
  88         catch(exception& e) {
  89                 m->errorOut(e, "AlignmentDB", "AlignmentDB");
  90                 exit(1);
  91         }
  92 }
  93 /**************************************************************************************************/
  94 AlignmentDB::~AlignmentDB() {  delete search;   }
  95 /**************************************************************************************************/
  96 Sequence AlignmentDB::findClosestSequence(Sequence* seq) {
  97         try{
  98
  99                 vector<int> spot = search->findClosestSequences(seq, 1);
 100
 101                 if (spot.size() != 0)   {               return templateSequences[spot[0]];      }
 102                 else                                    {               return emptySequence;                           }
 103
 104         }
 105         catch(exception& e) {
 106                 m->errorOut(e, "AlignmentDB", "findClosestSequence");
 107                 exit(1);
 108         }
 109 }
 110 /**************************************************************************************************/
 111
 112
 113
 114