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"
35 #include "distancedb.hpp"
38 #include "nastreport.hpp"
41 //**********************************************************************************************************************
42 AlignCommand::AlignCommand(){
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 = "????";
56 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
60 cout << "An unknown error has occurred in the AlignCommand class function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
65 //**********************************************************************************************************************
67 AlignCommand::~AlignCommand(){
71 //**********************************************************************************************************************
73 int AlignCommand::execute(){
75 srand( (unsigned)time( NULL ) ); //needed to assign names to temporary files
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); }
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); }
93 int numFastaSeqs=count(istreambuf_iterator<char>(in),istreambuf_iterator<char>(), '>');
96 string candidateAligngmentFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.align";
97 ofstream candidateAlignmentFile;
98 openOutputFile(candidateAligngmentFName, candidateAlignmentFile);
100 string candidateReportFName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + globaldata->getSearch() + '.' + globaldata->getAlign() + ".nast.report";
101 NastReport report(candidateReportFName);
103 cout << "We are going to align the " << numFastaSeqs << " sequences in " << candidateFileName << "..." << endl;
106 int start = time(NULL);
108 for(int i=0;i<numFastaSeqs;i++){
110 Sequence* candidateSeq = new Sequence(in);
111 report.setCandidate(candidateSeq);
113 Sequence* templateSeq = templateDB->findClosestSequence(candidateSeq);
114 report.setTemplate(templateSeq);
115 report.setSearchParameters(globaldata->getSearch(), templateDB->getSearchScore());
117 Nast nast(alignment, candidateSeq, templateSeq);
118 report.setAlignmentParameters(globaldata->getAlign(), alignment);
119 report.setNastParameters(nast);
121 candidateAlignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
122 candidateAlignmentFile.flush();
124 if((i+1) % 100 == 0){
125 cout << "It has taken " << time(NULL) - start << " secs to align " << i+1 << " sequences" << endl;
132 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
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";
146 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
151 //**********************************************************************************************************************