5 * Created by Sarah Westcott on 5/15/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
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).
18 * for basic instructions on how to run nastPlus type:
22 #include "aligncommand.h"
23 #include "sequence.hpp"
25 #include "alignment.hpp"
26 #include "gotohoverlap.hpp"
27 #include "needlemanoverlap.hpp"
28 #include "blastalign.hpp"
29 #include "noalign.hpp"
31 #include "database.hpp"
33 #include "suffixdb.hpp"
34 #include "blastdb.hpp"
37 #include "nastreport.hpp"
40 //**********************************************************************************************************************
41 AlignCommand::AlignCommand(){
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);
54 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
58 cout << "An unknown error has occurred in the AlignCommand class function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
63 //**********************************************************************************************************************
65 AlignCommand::~AlignCommand(){
69 //**********************************************************************************************************************
71 int AlignCommand::execute(){
73 srand( (unsigned)time( NULL ) ); //needed to assign names to temporary files
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); }
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); }
90 int numFastaSeqs=count(istreambuf_iterator<char>(in),istreambuf_iterator<char>(), '>');
93 string candidateAligngmentFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.align";
94 ofstream candidateAlignmentFile;
95 openOutputFile(candidateAligngmentFName, candidateAlignmentFile);
97 string candidateReportFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.report";
98 NastReport report(candidateReportFName);
100 cout << "We are going to align the " << numFastaSeqs << " sequences in " << candidateFileName << "..." << endl;
103 int start = time(NULL);
105 for(int i=0;i<numFastaSeqs;i++){
107 Sequence* candidateSeq = new Sequence(in);
108 report.setCandidate(candidateSeq);
110 Sequence* templateSeq = templateDB->findClosestSequence(candidateSeq);
111 report.setTemplate(templateSeq);
112 report.setSearchParameters(globaldata->getSearch(), templateDB->getSearchScore());
114 Nast nast(alignment, candidateSeq, templateSeq);
115 report.setAlignmentParameters(globaldata->getAlign(), alignment);
116 report.setNastParameters(nast);
118 candidateAlignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
119 candidateAlignmentFile.flush();
121 if((i+1) % 100 == 0){
122 cout << "It has taken " << time(NULL) - start << " secs to align " << i+1 << " sequences" << endl;
129 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
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";
143 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
148 //**********************************************************************************************************************