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                 if(globaldata->getFastaFile() == ""){
  45                         cout << "you forgot a template file" << endl;
  46                 }
  47                 openInputFile(globaldata->getCandidateFile(), in);
  48
  49                 convert(globaldata->getKSize(), kmerSize);
  50                 convert(globaldata->getMatch(), match);
  51                 convert(globaldata->getMismatch(), misMatch);
  52                 convert(globaldata->getGapopen(), gapOpen);
  53                 convert(globaldata->getGapextend(), gapExtend);
  54         }
  55         catch(exception& e) {
  56                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  57                 exit(1);
  58         }
  59         catch(...) {
  60                 cout << "An unknown error has occurred in the AlignCommand class function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  61                 exit(1);
  62         }
  63 }
  64
  65 //**********************************************************************************************************************
  66
  67 AlignCommand::~AlignCommand(){
  68
  69 }
  70
  71 //**********************************************************************************************************************
  72
  73 int AlignCommand::execute(){
  74         try {
  75                 srand( (unsigned)time( NULL ) );  //needed to assign names to temporary files
  76
  77                 Database* templateDB;
  78                 if(globaldata->getSearch() == "kmer")                   {       templateDB = new KmerDB(globaldata->getFastaFile() , kmerSize); }
  79                 else if(globaldata->getSearch() == "suffix")    {       templateDB = new SuffixDB(globaldata->getFastaFile());                  }
  80                 else if(globaldata->getSearch() == "blast")             {       templateDB = new BlastDB(globaldata->getFastaFile(), gapOpen, gapExtend, match, misMatch);      }
  81                 else {
  82                         cout << globaldata->getSearch() << " is not a valid search option. I will run the command using kmer, ksize=8." << endl;
  83                         templateDB = new KmerDB(globaldata->getFastaFile(), kmerSize);
  84                 }
  85
  86                 Alignment* alignment;
  87                 if(globaldata->getAlign() == "gotoh")                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
  88                 else if(globaldata->getAlign() == "needleman")  {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
  89                 else if(globaldata->getAlign() == "blast")              {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
  90                 else if(globaldata->getAlign() == "noalign")    {       alignment = new NoAlign();                                                                                                      }
  91                 else {
  92                         cout << globaldata->getAlign() << " is not a valid alignment option. I will run the command using needleman." << endl;
  93                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);
  94                 }
  95
  96                 int numFastaSeqs=count(istreambuf_iterator<char>(in),istreambuf_iterator<char>(), '>');
  97                 in.seekg(0);
  98
  99                 candidateFileName = globaldata->getCandidateFile();
 100                 string candidateAligngmentFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align";
 101                 ofstream candidateAlignmentFile;
 102                 openOutputFile(candidateAligngmentFName, candidateAlignmentFile);
 103
 104                 string candidateReportFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align.report";
 105                 NastReport report(candidateReportFName);
 106
 107                 cout << "We are going to align the " << numFastaSeqs << " sequences in " << candidateFileName << "..." << endl;
 108                 cout.flush();
 109
 110                 int start = time(NULL);
 111
 112                 for(int i=0;i<numFastaSeqs;i++){
 113
 114                         Sequence* candidateSeq = new Sequence(in);
 115                         report.setCandidate(candidateSeq);
 116
 117                         Sequence* templateSeq = templateDB->findClosestSequence(candidateSeq);
 118                         report.setTemplate(templateSeq);
 119                         report.setSearchParameters(globaldata->getSearch(), templateDB->getSearchScore());
 120
 121
 122                         Nast nast(alignment, candidateSeq, templateSeq);
 123                         report.setAlignmentParameters(globaldata->getAlign(), alignment);
 124                         report.setNastParameters(nast);
 125
 126                         candidateAlignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
 127                         candidateAlignmentFile.flush();
 128
 129                         if((i+1) % 100 == 0){
 130                                 cout << "It has taken " << time(NULL) - start << " secs to align " << i+1 << " sequences" << endl;
 131                         }
 132                         report.print();
 133
 134                         delete candidateSeq;
 135                 }
 136
 137                 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
 138                 cout << endl;
 139
 140                 delete templateDB;
 141                 delete alignment;
 142
 143
 144                 return 0;
 145         }
 146         catch(exception& e) {
 147                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 148                 exit(1);
 149         }
 150         catch(...) {
 151                 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 152                 exit(1);
 153         }
 154 }
 155
 156 //**********************************************************************************************************************