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 #include "referencedb.h"
  15
  16 /**************************************************************************************************/
  17 AlignmentDB::AlignmentDB(string fastaFileName, string s, int kmerSize, float gapOpen, float gapExtend, float match, float misMatch, int tid){           //      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                 method = s;
  22                 bool needToGenerate = true;
  23                 ReferenceDB* rdb = ReferenceDB::getInstance();
  24                 bool silent = false;
  25                 threadID = tid;
  26
  27                 if (fastaFileName == "saved-silent") {
  28                         fastaFileName = "saved"; silent = true;
  29                 }
  30
  31                 if (fastaFileName == "saved") {
  32                         int start = time(NULL);
  33
  34                         if (!silent) { m->mothurOutEndLine();  m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory."); m->mothurOutEndLine(); }
  35
  36                         for (int i = 0; i < rdb->referenceSeqs.size(); i++) {
  37                                 templateSequences.push_back(rdb->referenceSeqs[i]);
  38                                 //save longest base
  39                                 if (rdb->referenceSeqs[i].getUnaligned().length() >= longest)  { longest = (rdb->referenceSeqs[i].getUnaligned().length()+1); }
  40                         }
  41                         fastaFileName = rdb->getSavedReference();
  42
  43                         numSeqs = templateSequences.size();
  44                         if (!silent) { m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();  }
  45
  46                 }else {
  47                         int start = time(NULL);
  48                         m->mothurOutEndLine();
  49                         m->mothurOut("Reading in the " + fastaFileName + " template sequences...\t");   cout.flush();
  50
  51                         #ifdef USE_MPI
  52                                 int pid, processors;
  53                                 vector<unsigned long long> positions;
  54
  55                                 MPI_Status status;
  56                                 MPI_File inMPI;
  57                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
  58                                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
  59                                 int tag = 2001;
  60
  61                                 char inFileName[1024];
  62                                 strcpy(inFileName, fastaFileName.c_str());
  63
  64                                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
  65
  66                                 if (pid == 0) {
  67                                         positions = m->setFilePosFasta(fastaFileName, numSeqs); //fills MPIPos, returns numSeqs
  68
  69                                         //send file positions to all processes
  70                                         for(int i = 1; i < processors; i++) {
  71                                                 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
  72                                                 MPI_Send(&positions[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
  73                                         }
  74                                 }else{
  75                                         MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
  76                                         positions.resize(numSeqs+1);
  77                                         MPI_Recv(&positions[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
  78                                 }
  79
  80                                 //read file
  81                                 for(int i=0;i<numSeqs;i++){
  82
  83                                         if (m->control_pressed) {  templateSequences.clear(); break;  }
  84
  85                                         //read next sequence
  86                                         int length = positions[i+1] - positions[i];
  87                                         char* buf4 = new char[length];
  88
  89                                         MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
  90
  91                                         string tempBuf = buf4;
  92                                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
  93                                         delete buf4;
  94
  95                                         istringstream iss (tempBuf,istringstream::in);
  96
  97                                         Sequence temp(iss);
  98                                         if (temp.getName() != "") {
  99                                                 templateSequences.push_back(temp);
 100
 101                                                 if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
 102
 103                                                 //save longest base
 104                                                 if (temp.getUnaligned().length() >= longest)  { longest = temp.getUnaligned().length()+1; }
 105                                         }
 106                                 }
 107
 108                                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 109
 110                                 MPI_File_close(&inMPI);
 111
 112                 #else
 113                         ifstream fastaFile;
 114                         m->openInputFile(fastaFileName, fastaFile);
 115
 116                         while (!fastaFile.eof()) {
 117                                 Sequence temp(fastaFile);  m->gobble(fastaFile);
 118
 119                                 if (m->control_pressed) {  templateSequences.clear(); break;  }
 120
 121                                 if (temp.getName() != "") {
 122                                         templateSequences.push_back(temp);
 123
 124                                         if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
 125
 126                                         //save longest base
 127                                         if (temp.getUnaligned().length() >= longest)  { longest = (temp.getUnaligned().length()+1); }
 128                                 }
 129                         }
 130                         fastaFile.close();
 131                 #endif
 132
 133                         numSeqs = templateSequences.size();
 134                         //all of this is elsewhere already!
 135
 136                         m->mothurOut("DONE.");
 137                         m->mothurOutEndLine();  cout.flush();
 138                         m->mothurOut("It took " + toString(time(NULL) - start) + " to read  " + toString(templateSequences.size()) + " sequences."); m->mothurOutEndLine();
 139
 140                 }
 141
 142
 143                 //in case you delete the seqs and then ask for them
 144                 emptySequence = Sequence();
 145                 emptySequence.setName("no_match");
 146                 emptySequence.setUnaligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
 147                 emptySequence.setAligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
 148
 149
 150                 string kmerDBName;
 151                 if(method == "kmer")                    {
 152                         search = new KmerDB(fastaFileName, kmerSize);
 153
 154                         #ifdef USE_MPI
 155                         #else
 156                                 kmerDBName = fastaFileName.substr(0,fastaFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
 157                                 ifstream kmerFileTest(kmerDBName.c_str());
 158
 159                                 if(kmerFileTest){
 160                                         bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion());
 161                                         if (GoodFile) {  needToGenerate = false;        }
 162                                 }
 163                         #endif
 164                 }
 165                 else if(method == "suffix")             {       search = new SuffixDB(numSeqs);                                                         }
 166                 else if(method == "blast")              {       search = new BlastDB(fastaFileName.substr(0,fastaFileName.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch, "", threadID);     }
 167                 else {
 168                         method = "kmer";
 169                         m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
 170                         m->mothurOutEndLine();
 171                         search = new KmerDB(fastaFileName, 8);
 172                 }
 173
 174                 if (!(m->control_pressed)) {
 175                         if (needToGenerate) {
 176                                 //add sequences to search
 177                                 for (int i = 0; i < templateSequences.size(); i++) {
 178                                         search->addSequence(templateSequences[i]);
 179
 180                                         if (m->control_pressed) {  templateSequences.clear(); break;  }
 181                                 }
 182
 183                                 if (m->control_pressed) {  templateSequences.clear();  }
 184
 185                                 search->generateDB();
 186
 187                         }else if ((method == "kmer") && (!needToGenerate)) {
 188                                 ifstream kmerFileTest(kmerDBName.c_str());
 189                                 search->readKmerDB(kmerFileTest);
 190                         }
 191
 192                         search->setNumSeqs(numSeqs);
 193                 }
 194
 195         }
 196         catch(exception& e) {
 197                 m->errorOut(e, "AlignmentDB", "AlignmentDB");
 198                 exit(1);
 199         }
 200 }
 201 /**************************************************************************************************/
 202 AlignmentDB::AlignmentDB(string s){
 203         try {
 204                 m = MothurOut::getInstance();
 205                 method = s;
 206
 207                 if(method == "suffix")          {       search = new SuffixDB();        }
 208                 else if(method == "blast")      {       search = new BlastDB("", 0);            }
 209                 else                                            {       search = new KmerDB();          }
 210
 211
 212                 //in case you delete the seqs and then ask for them
 213                 emptySequence = Sequence();
 214                 emptySequence.setName("no_match");
 215                 emptySequence.setUnaligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
 216                 emptySequence.setAligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
 217
 218         }
 219         catch(exception& e) {
 220                 m->errorOut(e, "AlignmentDB", "AlignmentDB");
 221                 exit(1);
 222         }
 223 }
 224 /**************************************************************************************************/
 225 AlignmentDB::~AlignmentDB() {  delete search;   }
 226 /**************************************************************************************************/
 227 Sequence AlignmentDB::findClosestSequence(Sequence* seq) {
 228         try{
 229
 230                 vector<int> spot = search->findClosestSequences(seq, 1);
 231
 232                 if (spot.size() != 0)   {               return templateSequences[spot[0]];      }
 233                 else                                    {               return emptySequence;                           }
 234
 235         }
 236         catch(exception& e) {
 237                 m->errorOut(e, "AlignmentDB", "findClosestSequence");
 238                 exit(1);
 239         }
 240 }
 241 /**************************************************************************************************/
 242
 243
 244
 245
 246
 247