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(string option){
  42         try {
  43                 globaldata = GlobalData::getInstance();
  44                 abort = false;
  45
  46                 //allow user to run help
  47                 if(option == "help") { help(); abort = true; }
  48
  49                 else {
  50                         //valid paramters for this command
  51                         string AlignArray[] =  {"fasta","candidate","search","ksize","align","match","mismatch","gapopen","gapextend"};
  52                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  53
  54                         parser = new OptionParser();
  55
  56                         parser->parse(option, parameters);      delete parser;
  57
  58                         ValidParameters* validParameter = new ValidParameters();
  59
  60                         //check to make sure all parameters are valid for command
  61                         for (it = parameters.begin(); it != parameters.end(); it++) {
  62                                 if (validParameter->isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  63                         }
  64
  65                         //check for required parameters
  66                         templateFileName = validParameter->validFile(parameters, "fasta", true);
  67                         if (templateFileName == "not found") { cout << "fasta is a required parameter for the align.seqs command." << endl; abort = true; }
  68                         else if (templateFileName == "not open") { abort = true; }
  69                         else { globaldata->setFastaFile(templateFileName); }
  70
  71                         candidateFileName = validParameter->validFile(parameters, "candidate", true);
  72                         if (candidateFileName == "not found") { cout << "candidate is a required parameter for the align.seqs command." << endl; abort = true; }
  73                         else if (candidateFileName == "not open") { abort = true; }
  74                         else {
  75                                 globaldata->setCandidateFile(candidateFileName);
  76                                 openInputFile(candidateFileName, in);
  77                         }
  78
  79
  80                         //check for optional parameter and set defaults
  81                         // ...at some point should added some additional type checking...
  82                         string temp;
  83                         temp = validParameter->validFile(parameters, "ksize", false);                   if (temp == "not found") { temp = "8"; }
  84                         convert(temp, kmerSize);
  85
  86                         temp = validParameter->validFile(parameters, "match", false);                   if (temp == "not found") { temp = "1.0"; }
  87                         convert(temp, match);
  88
  89                         temp = validParameter->validFile(parameters, "mismatch", false);                        if (temp == "not found") { temp = "-1.0"; }
  90                         convert(temp, misMatch);
  91
  92                         temp = validParameter->validFile(parameters, "gapopen", false);                 if (temp == "not found") { temp = "-1.0"; }
  93                         convert(temp, gapOpen);
  94
  95                         temp = validParameter->validFile(parameters, "gapextend", false);               if (temp == "not found") { temp = "-2.0"; }
  96                         convert(temp, gapExtend);
  97
  98                         search = validParameter->validFile(parameters, "search", false);                        if (search == "not found")      { search = "kmer";              }
  99                         align = validParameter->validFile(parameters, "align", false);                  if (align == "not found")       { align = "needleman";  }
 100
 101                         delete validParameter;
 102                 }
 103
 104         }
 105         catch(exception& e) {
 106                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 107                 exit(1);
 108         }
 109         catch(...) {
 110                 cout << "An unknown error has occurred in the AlignCommand class function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 111                 exit(1);
 112         }
 113 }
 114
 115 //**********************************************************************************************************************
 116
 117 AlignCommand::~AlignCommand(){
 118
 119 }
 120 //**********************************************************************************************************************
 121
 122 void AlignCommand::help(){
 123         try {
 124                 cout << "The align.seqs command reads a file containing sequences and creates an alignment file and a report file." << "\n";
 125                 cout << "The align.seqs command parameters are fasta, candidate, search, ksize, align, match, mismatch, gapopen and gapextend.  " << "\n";
 126                 cout << "The fasta and candidate parameters are required." << "\n";
 127                 cout << "The search parameter allows you to specify the method to find most similar template.  Your options are: suffix, kmer and blast. The default is kmer." << "\n";
 128                 cout << "The align parameter allows you to specify the alignment method to use.  Your options are: gotoh, needleman, blast and noalign. The default is needleman." << "\n";
 129                 cout << "The ksize parameter allows you to specify the kmer size for finding most similar template to candidate.  The default is 7." << "\n";
 130                 cout << "The match parameter allows you to specify the bonus for having the same base. The default is 1.0." << "\n";
 131                 cout << "The mistmatch parameter allows you to specify the penalty for having different bases.  The default is -1.0." << "\n";
 132                 cout << "The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -1.0." << "\n";
 133                 cout << "The gapextend parameter allows you to specify the penalty for extending a gap in an alignment.  The default is -2.0." << "\n";
 134                 cout << "The align.seqs command should be in the following format: " << "\n";
 135                 cout << "align.seqs(fasta=yourTemplateFile, candidate=yourCandidateFile, align=yourAlignmentMethod, search=yourSearchmethod, ksize=yourKmerSize, match=yourMatchBonus, mismatch=yourMismatchpenalty, gapopen=yourGapopenPenalty, gapextend=yourGapExtendPenalty) " << "\n";
 136                 cout << "Example align.seqs(candidate=candidate.fasta, fasta=core.filtered, align=kmer, search=gotoh, ksize=8, match=2.0, mismatch=3.0, gapopen=-2.0, gapextend=-1.0)" << "\n";
 137                 cout << "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile)." << "\n" << "\n";
 138         }
 139         catch(exception& e) {
 140                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 141                 exit(1);
 142         }
 143         catch(...) {
 144                 cout << "An unknown error has occurred in the AlignCommand class function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 145                 exit(1);
 146         }
 147 }
 148
 149
 150 //**********************************************************************************************************************
 151
 152 int AlignCommand::execute(){
 153         try {
 154                 if (abort == true) {    return 0;       }
 155
 156                 srand( (unsigned)time( NULL ) );  //needed to assign names to temporary files
 157
 158                 Database* templateDB;
 159                 if(search == "kmer")                    {       templateDB = new KmerDB(templateFileName, kmerSize);    }
 160                 else if(search == "suffix")             {       templateDB = new SuffixDB(templateFileName);                    }
 161                 else if(search == "blast")              {       templateDB = new BlastDB(templateFileName, gapOpen, gapExtend, match, misMatch);        }
 162                 else {
 163                         cout << search << " is not a valid search option. I will run the command using kmer, ksize=8." << endl; kmerSize = 8;
 164                         templateDB = new KmerDB(templateFileName, kmerSize);
 165                 }
 166
 167                 Alignment* alignment;
 168                 if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
 169                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
 170                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 171                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 172                 else {
 173                         cout << align << " is not a valid alignment option. I will run the command using needleman." << endl;
 174                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);
 175                 }
 176
 177                 int numFastaSeqs=count(istreambuf_iterator<char>(in),istreambuf_iterator<char>(), '>');
 178                 in.seekg(0);
 179
 180                 string candidateAligngmentFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align";
 181                 ofstream candidateAlignmentFile;
 182                 openOutputFile(candidateAligngmentFName, candidateAlignmentFile);
 183
 184                 string candidateReportFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align.report";
 185                 NastReport report(candidateReportFName);
 186
 187                 cout << "We are going to align the " << numFastaSeqs << " sequences in " << candidateFileName << "..." << endl;
 188                 cout.flush();
 189
 190                 int start = time(NULL);
 191
 192                 for(int i=0;i<numFastaSeqs;i++){
 193
 194                         Sequence* candidateSeq = new Sequence(in);
 195                         report.setCandidate(candidateSeq);
 196
 197                         Sequence* templateSeq = templateDB->findClosestSequence(candidateSeq);
 198                         report.setTemplate(templateSeq);
 199                         report.setSearchParameters(search, templateDB->getSearchScore());
 200
 201
 202                         Nast nast(alignment, candidateSeq, templateSeq);
 203                         report.setAlignmentParameters(align, alignment);
 204                         report.setNastParameters(nast);
 205
 206                         candidateAlignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
 207                         candidateAlignmentFile.flush();
 208
 209                         if((i+1) % 100 == 0){
 210                                 cout << "It has taken " << time(NULL) - start << " secs to align " << i+1 << " sequences" << endl;
 211                         }
 212                         report.print();
 213
 214                         delete candidateSeq;
 215                 }
 216
 217                 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
 218                 cout << endl;
 219
 220                 delete templateDB;
 221                 delete alignment;
 222
 223
 224                 return 0;
 225         }
 226         catch(exception& e) {
 227                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 228                 exit(1);
 229         }
 230         catch(...) {
 231                 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 232                 exit(1);
 233         }
 234 }
 235
 236 //**********************************************************************************************************************