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 #include "distancedb.hpp"
  36
  37 #include "nast.hpp"
  38 #include "nastreport.hpp"
  39
  40
  41 //**********************************************************************************************************************
  42 AlignCommand::AlignCommand(){
  43         try {
  44                 globaldata = GlobalData::getInstance();
  45                 candidateFileName = globaldata->inputFileName;
  46                 templateFileName = globaldata->getTemplateFile();
  47                 openInputFile(candidateFileName, in);
  48                 convert(globaldata->getKSize(), kmerSize);
  49                 convert(globaldata->getMatch(), match);
  50                 convert(globaldata->getMismatch(), misMatch);
  51                 convert(globaldata->getGapopen(), gapOpen);
  52                 convert(globaldata->getGapextend(), gapExtend);
  53                 distanceFileName = "????";
  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(templateFileName, kmerSize);                                                            }
  79                 else if(globaldata->getSearch() == "suffix")    {       templateDB = new SuffixDB(templateFileName);                                                                            }
  80                 else if(globaldata->getSearch() == "blast")             {       templateDB = new BlastDB(templateFileName, gapOpen, gapExtend, match, misMatch);        }
  81                 else if(globaldata->getSearch() == "distance")  {       templateDB = new DistanceDB(templateFileName, distanceFileName);                                        }
  82                 else {  cout << globaldata->getSearch() << " is not a valid search option. I will run the command using suffix." << endl;
  83                                 templateDB = new SuffixDB(templateFileName);            }
  84
  85                 Alignment* alignment;
  86                 if(globaldata->getAlign() == "gotoh")                                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
  87                 else if(globaldata->getAlign() == "needleman")                  {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
  88                 else if(globaldata->getAlign() == "blast")                              {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
  89                 else if(globaldata->getAlign() == "noalign")                    {       alignment = new NoAlign();                                                                                                      }
  90                 else {  cout << globaldata->getAlign() << " is not a valid alignment option. I will run the command using blast." << endl;
  91                                 alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);    }
  92
  93                 int numFastaSeqs=count(istreambuf_iterator<char>(in),istreambuf_iterator<char>(), '>');
  94                 in.seekg(0);
  95
  96                 string candidateAligngmentFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.align";
  97                 ofstream candidateAlignmentFile;
  98                 openOutputFile(candidateAligngmentFName, candidateAlignmentFile);
  99
 100                 string candidateReportFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.report";
 101                 NastReport report(candidateReportFName);
 102
 103                 cout << "We are going to align the " << numFastaSeqs << " sequences in " << candidateFileName << "..." << endl;
 104                 cout.flush();
 105
 106                 int start = time(NULL);
 107
 108                 for(int i=0;i<numFastaSeqs;i++){
 109
 110                         Sequence* candidateSeq = new Sequence(in);
 111                         report.setCandidate(candidateSeq);
 112
 113                         Sequence* templateSeq = templateDB->findClosestSequence(candidateSeq);
 114                         report.setTemplate(templateSeq);
 115                         report.setSearchParameters(globaldata->getSearch(), templateDB->getSearchScore());
 116
 117                         Nast nast(alignment, candidateSeq, templateSeq);
 118                         report.setAlignmentParameters(globaldata->getAlign(), alignment);
 119                         report.setNastParameters(nast);
 120
 121                         candidateAlignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
 122                         candidateAlignmentFile.flush();
 123
 124                         if((i+1) % 100 == 0){
 125                                 cout << "It has taken " << time(NULL) - start << " secs to align " << i+1 << " sequences" << endl;
 126                         }
 127                         report.print();
 128
 129                         delete candidateSeq;
 130                 }
 131
 132                 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
 133                 cout << endl;
 134
 135                 delete templateDB;
 136                 delete alignment;
 137
 138
 139                 return 0;
 140         }
 141         catch(exception& e) {
 142                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 143                 exit(1);
 144         }
 145         catch(...) {
 146                 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 147                 exit(1);
 148         }
 149 }
 150
 151 //**********************************************************************************************************************