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