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