aligncommand.cpp

   1 /*
   2  *  aligncommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 5/15/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  *      This version of nast does everything I think that the greengenes nast server does and then some.  I have added the
   9  *      feature of allowing users to define their database, kmer size for searching, alignment penalty values and alignment
  10  *      method.  This latter feature is perhaps most significant.  nastPlus enables a user to use either a Needleman-Wunsch
  11  *      (non-affine gap penalty) or Gotoh (affine gap penalty) pairwise alignment algorithm.  This is significant because it
  12  *      allows for a global alignment and not the local alignment provided by bLAst.  Furthermore, it has the potential to
  13  *      provide a better alignment because of the banding method employed by blast (I'm not sure about this).
  14  *
  15  *      to compile type:
  16  *              make
  17  *
  18  *      for basic instructions on how to run nastPlus type:
  19  *              ./nastPlus
  20  */
  21
  22 #include "aligncommand.h"
  23 #include "sequence.hpp"
  24
  25 #include "alignment.hpp"
  26 #include "gotohoverlap.hpp"
  27 #include "needlemanoverlap.hpp"
  28 #include "blastalign.hpp"
  29 #include "noalign.hpp"
  30
  31 #include "database.hpp"
  32 #include "kmerdb.hpp"
  33 #include "suffixdb.hpp"
  34 #include "blastdb.hpp"
  35
  36 #include "nast.hpp"
  37 #include "nastreport.hpp"
  38
  39
  40 //**********************************************************************************************************************
  41 AlignCommand::AlignCommand(){
  42         try {
  43                 globaldata = GlobalData::getInstance();
  44                 candidateFileName = globaldata->inputFileName;
  45                 templateFileName = globaldata->getTemplateFile();
  46                 openInputFile(candidateFileName, in);
  47                 convert(globaldata->getKSize(), kmerSize);
  48                 convert(globaldata->getMatch(), match);
  49                 convert(globaldata->getMismatch(), misMatch);
  50                 convert(globaldata->getGapopen(), gapOpen);
  51                 convert(globaldata->getGapextend(), gapExtend);
  52         }
  53         catch(exception& e) {
  54                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  55                 exit(1);
  56         }
  57         catch(...) {
  58                 cout << "An unknown error has occurred in the AlignCommand class function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  59                 exit(1);
  60         }
  61 }
  62
  63 //**********************************************************************************************************************
  64
  65 AlignCommand::~AlignCommand(){
  66
  67 }
  68
  69 //**********************************************************************************************************************
  70
  71 int AlignCommand::execute(){
  72         try {
  73                 srand( (unsigned)time( NULL ) );  //needed to assign names to temporary files
  74
  75                 Database* templateDB;
  76                 if(globaldata->getSearch() == "kmer")                   {       templateDB = new KmerDB(templateFileName, kmerSize);                                                            }
  77                 else if(globaldata->getSearch() == "suffix")    {       templateDB = new SuffixDB(templateFileName);                                                                            }
  78                 else if(globaldata->getSearch() == "blast")             {       templateDB = new BlastDB(templateFileName, gapOpen, gapExtend, match, misMatch);        }
  79                 else {  cout << globaldata->getSearch() << " is not a valid search option. I will run the command using suffix." << endl;
  80                                 templateDB = new SuffixDB(templateFileName);            }
  81
  82                 Alignment* alignment;
  83                 if(globaldata->getAlign() == "gotoh")                                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
  84                 else if(globaldata->getAlign() == "needleman")                  {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
  85                 else if(globaldata->getAlign() == "blast")                              {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
  86                 else if(globaldata->getAlign() == "noalign")                    {       alignment = new NoAlign();                                                                                                      }
  87                 else {  cout << globaldata->getAlign() << " is not a valid alignment option. I will run the command using blast." << endl;
  88                                 alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);    }
  89
  90                 int numFastaSeqs=count(istreambuf_iterator<char>(in),istreambuf_iterator<char>(), '>');
  91                 in.seekg(0);
  92
  93                 string candidateAligngmentFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.align";
  94                 ofstream candidateAlignmentFile;
  95                 openOutputFile(candidateAligngmentFName, candidateAlignmentFile);
  96
  97                 string candidateReportFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.report";
  98                 NastReport report(candidateReportFName);
  99
 100                 cout << "We are going to align the " << numFastaSeqs << " sequences in " << candidateFileName << "..." << endl;
 101                 cout.flush();
 102
 103                 int start = time(NULL);
 104
 105                 for(int i=0;i<numFastaSeqs;i++){
 106
 107                         Sequence* candidateSeq = new Sequence(in);
 108                         report.setCandidate(candidateSeq);
 109
 110                         Sequence* templateSeq = templateDB->findClosestSequence(candidateSeq);
 111                         report.setTemplate(templateSeq);
 112                         report.setSearchParameters(globaldata->getSearch(), templateDB->getSearchScore());
 113
 114                         Nast nast(alignment, candidateSeq, templateSeq);
 115                         report.setAlignmentParameters(globaldata->getAlign(), alignment);
 116                         report.setNastParameters(nast);
 117
 118                         candidateAlignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
 119                         candidateAlignmentFile.flush();
 120
 121                         if((i+1) % 100 == 0){
 122                                 cout << "It has taken " << time(NULL) - start << " secs to align " << i+1 << " sequences" << endl;
 123                         }
 124                         report.print();
 125
 126                         delete candidateSeq;
 127                 }
 128
 129                 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
 130                 cout << endl;
 131
 132                 delete templateDB;
 133                 delete alignment;
 134
 135
 136                 return 0;
 137         }
 138         catch(exception& e) {
 139                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 140                 exit(1);
 141         }
 142         catch(...) {
 143                 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 144                 exit(1);
 145         }
 146 }
 147
 148 //**********************************************************************************************************************