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  */
  16
  17 #include "aligncommand.h"
  18 #include "sequence.hpp"
  19
  20 #include "gotohoverlap.hpp"
  21 #include "needlemanoverlap.hpp"
  22 #include "blastalign.hpp"
  23 #include "noalign.hpp"
  24
  25 #include "kmerdb.hpp"
  26 #include "suffixdb.hpp"
  27 #include "blastdb.hpp"
  28
  29 #include "nast.hpp"
  30 #include "nastreport.hpp"
  31
  32
  33 //**********************************************************************************************************************
  34
  35 AlignCommand::AlignCommand(string option){
  36         try {
  37                 //              globaldata = GlobalData::getInstance();
  38                 abort = false;
  39
  40                 //allow user to run help
  41                 if(option == "help") { help(); abort = true; }
  42
  43                 else {
  44
  45                         //valid paramters for this command
  46                         string AlignArray[] =  {"template","candidate","search","ksize","align","match","mismatch","gapopen","gapextend", "processors"};
  47                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  48
  49                         OptionParser parser(option);
  50                         map<string, string> parameters = parser.getParameters();
  51
  52                         ValidParameters validParameter;
  53
  54                         //check to make sure all parameters are valid for command
  55                         for (map<string, string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  56                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  57                         }
  58
  59                         //check for required parameters
  60                         templateFileName = validParameter.validFile(parameters, "template", true);
  61                         if (templateFileName == "not found") { cout << "template is a required parameter for the align.seqs command." << endl; abort = true; }
  62                         else if (templateFileName == "not open") { abort = true; }
  63
  64                         candidateFileName = validParameter.validFile(parameters, "candidate", true);
  65                         if (candidateFileName == "not found") { cout << "candidate is a required parameter for the align.seqs command." << endl; abort = true; }
  66                         else if (candidateFileName == "not open") { abort = true; }
  67
  68                         //check for optional parameter and set defaults
  69                         // ...at some point should added some additional type checking...
  70                         string temp;
  71                         temp = validParameter.validFile(parameters, "ksize", false);            if (temp == "not found"){       temp = "8";                             }
  72                         convert(temp, kmerSize);
  73
  74                         temp = validParameter.validFile(parameters, "match", false);            if (temp == "not found"){       temp = "1.0";                   }
  75                         convert(temp, match);
  76
  77                         temp = validParameter.validFile(parameters, "mismatch", false);         if (temp == "not found"){       temp = "-1.0";                  }
  78                         convert(temp, misMatch);
  79
  80                         temp = validParameter.validFile(parameters, "gapopen", false);          if (temp == "not found"){       temp = "-2.0";                  }
  81                         convert(temp, gapOpen);
  82
  83                         temp = validParameter.validFile(parameters, "gapextend", false);        if (temp == "not found"){       temp = "-1.0";                  }
  84                         convert(temp, gapExtend);
  85
  86                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
  87                         convert(temp, processors);
  88
  89                         search = validParameter.validFile(parameters, "search", false);         if (search == "not found"){     search = "kmer";                }
  90
  91                         align = validParameter.validFile(parameters, "align", false);           if (align == "not found"){      align = "needleman";    }
  92                 }
  93
  94         }
  95         catch(exception& e) {
  96                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  97                 exit(1);
  98         }
  99         catch(...) {
 100                 cout << "An unknown error has occurred in the AlignCommand class function AlignCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 101                 exit(1);
 102         }
 103 }
 104
 105 //**********************************************************************************************************************
 106
 107 AlignCommand::~AlignCommand(){
 108         delete templateDB;
 109         delete alignment;
 110 }
 111
 112 //**********************************************************************************************************************
 113
 114 void AlignCommand::help(){
 115         try {
 116                 cout << "The align.seqs command reads a file containing sequences and creates an alignment file and a report file." << "\n";
 117                 cout << "The align.seqs command parameters are template, candidate, search, ksize, align, match, mismatch, gapopen and gapextend.  " << "\n";
 118                 cout << "The template and candidate parameters are required." << "\n";
 119                 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";
 120                 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";
 121                 cout << "The ksize parameter allows you to specify the kmer size for finding most similar template to candidate.  The default is 7." << "\n";
 122                 cout << "The match parameter allows you to specify the bonus for having the same base. The default is 1.0." << "\n";
 123                 cout << "The mistmatch parameter allows you to specify the penalty for having different bases.  The default is -1.0." << "\n";
 124                 cout << "The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -1.0." << "\n";
 125                 cout << "The gapextend parameter allows you to specify the penalty for extending a gap in an alignment.  The default is -2.0." << "\n";
 126                 cout << "The align.seqs command should be in the following format: " << "\n";
 127                 cout << "align.seqs(template=yourTemplateFile, candidate=yourCandidateFile, align=yourAlignmentMethod, search=yourSearchmethod, ksize=yourKmerSize, match=yourMatchBonus, mismatch=yourMismatchpenalty, gapopen=yourGapopenPenalty, gapextend=yourGapExtendPenalty) " << "\n";
 128                 cout << "Example align.seqs(candidate=candidate.fasta, template=core.filtered, align=kmer, search=gotoh, ksize=8, match=2.0, mismatch=3.0, gapopen=-2.0, gapextend=-1.0)" << "\n";
 129                 cout << "Note: No spaces between parameter labels (i.e. candidate), '=' and parameters (i.e.yourFastaFile)." << "\n" << "\n";
 130         }
 131         catch(exception& e) {
 132                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 133                 exit(1);
 134         }
 135         catch(...) {
 136                 cout << "An unknown error has occurred in the AlignCommand class function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 137                 exit(1);
 138         }
 139 }
 140
 141
 142 //**********************************************************************************************************************
 143
 144 int AlignCommand::execute(){
 145         try {
 146                 if (abort == true) {    return 0;       }
 147
 148                 if(search == "kmer")                    {       templateDB = new KmerDB(templateFileName, kmerSize);    }
 149                 else if(search == "suffix")             {       templateDB = new SuffixDB(templateFileName);                    }
 150                 else if(search == "blast")              {       templateDB = new BlastDB(templateFileName, gapOpen, gapExtend, match, misMatch);        }
 151                 else {
 152                         cout << search << " is not a valid search option. I will run the command using kmer, ksize=8." << endl; kmerSize = 8;
 153                         templateDB = new KmerDB(templateFileName, kmerSize);
 154                 }
 155
 156                 if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
 157                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
 158                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 159                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 160                 else {
 161                         cout << align << " is not a valid alignment option. I will run the command using needleman." << endl;
 162                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);
 163                 }
 164
 165                 string alignFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align";
 166                 string reportFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align.report";
 167
 168                 int numFastaSeqs = 0;
 169                 int start = time(NULL);
 170
 171 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 172                 if(processors == 1){
 173                         ifstream inFASTA;
 174                         openInputFile(candidateFileName, inFASTA);
 175                         numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 176                         inFASTA.close();
 177
 178                         lines.push_back(new linePair(0, numFastaSeqs));
 179
 180                         driver(lines[0], alignFileName, reportFileName);
 181                 }
 182                 else{
 183                         vector<int> positions;
 184                         processIDS.resize(0);
 185
 186                         ifstream inFASTA;
 187                         openInputFile(candidateFileName, inFASTA);
 188
 189                         while(!inFASTA.eof()){
 190                                 char c = inFASTA.get();
 191                                 if(c == '>'){   positions.push_back(inFASTA.tellg());   }
 192                                 while (!inFASTA.eof())  {       c = inFASTA.get(); if (c == 10 || c == 13){     break;  }       } // get rest of line if there's any crap there
 193                         }
 194                         inFASTA.close();
 195
 196                         numFastaSeqs = positions.size();
 197
 198                         int numSeqsPerProcessor = numFastaSeqs / processors;
 199
 200                         for (int i = 0; i < processors; i++) {
 201                                 int startPos = positions[ i * numSeqsPerProcessor ];
 202                                 if(i == processors - 1){
 203                                         numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
 204                                 }
 205                                 lines.push_back(new linePair(startPos, numSeqsPerProcessor));
 206                         }
 207                         createProcesses(alignFileName, reportFileName);
 208
 209                         rename((alignFileName + toString(processIDS[0]) + ".temp").c_str(), alignFileName.c_str());
 210                         rename((reportFileName + toString(processIDS[0]) + ".temp").c_str(), reportFileName.c_str());
 211
 212                         for(int i=1;i<processors;i++){
 213                                 appendAlignFiles((alignFileName + toString(processIDS[i]) + ".temp"), alignFileName);
 214                                 remove((alignFileName + toString(processIDS[i]) + ".temp").c_str());
 215
 216                                 appendReportFiles((reportFileName + toString(processIDS[i]) + ".temp"), reportFileName);
 217                                 remove((reportFileName + toString(processIDS[i]) + ".temp").c_str());
 218                         }
 219                 }
 220 #else
 221                 ifstream inFASTA;
 222                 openInputFile(candidateFileName, inFASTA);
 223                 numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 224                 inFASTA.close();
 225
 226                 lines.push_back(new linePair(0, numFastaSeqs));
 227
 228                 driver(lines[0], alignFileName, reportFileName);
 229 #endif
 230
 231                 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
 232                 cout << endl;
 233
 234
 235
 236                 return 0;
 237         }
 238         catch(exception& e) {
 239                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 240                 exit(1);
 241         }
 242         catch(...) {
 243                 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 244                 exit(1);
 245         }
 246 }
 247
 248 //**********************************************************************************************************************
 249
 250 int AlignCommand::driver(linePair* line, string alignFName, string reportFName){
 251         try {
 252                 ofstream alignmentFile;
 253                 openOutputFile(alignFName, alignmentFile);
 254                 NastReport report(reportFName);
 255
 256                 ifstream inFASTA;
 257                 openInputFile(candidateFileName, inFASTA);
 258                 inFASTA.seekg(line->start);
 259
 260                 for(int i=0;i<line->numSeqs;i++){
 261
 262                         Sequence* candidateSeq = new Sequence(inFASTA);
 263                         report.setCandidate(candidateSeq);
 264
 265                         Sequence* templateSeq = templateDB->findClosestSequence(candidateSeq);
 266                         report.setTemplate(templateSeq);
 267                         report.setSearchParameters(search, templateDB->getSearchScore());
 268
 269                         Nast nast(alignment, candidateSeq, templateSeq);
 270                         report.setAlignmentParameters(align, alignment);
 271                         report.setNastParameters(nast);
 272
 273                         alignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
 274
 275                         report.print();
 276
 277                         delete candidateSeq;
 278                 }
 279
 280                 alignmentFile.close();
 281                 inFASTA.close();
 282
 283                 return 1;
 284         }
 285         catch(exception& e) {
 286                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 287                 exit(1);
 288         }
 289         catch(...) {
 290                 cout << "An unknown error has occurred in the AlignCommand class function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 291                 exit(1);
 292         }
 293 }
 294
 295 /**************************************************************************************************/
 296
 297 void AlignCommand::createProcesses(string alignFileName, string reportFileName) {
 298         try {
 299 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 300                 int process = 0;
 301                 //              processIDS.resize(0);
 302
 303                 //loop through and create all the processes you want
 304                 while (process != processors) {
 305                         int pid = fork();
 306
 307                         if (pid > 0) {
 308                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 309                                 process++;
 310                         }else if (pid == 0){
 311                                 driver(lines[process], alignFileName + toString(getpid()) + ".temp", reportFileName + toString(getpid()) + ".temp");
 312                                 exit(0);
 313                         }else { cout << "unable to spawn the necessary processes." << endl; exit(0); }
 314                 }
 315
 316                 //force parent to wait until all the processes are done
 317                 for (int i=0;i<processors;i++) {
 318                         int temp = processIDS[i];
 319                         wait(&temp);
 320                 }
 321 #endif
 322         }
 323         catch(exception& e) {
 324                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 325                 exit(1);
 326         }
 327         catch(...) {
 328                 cout << "An unknown error has occurred in the AlignCommand class function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 329                 exit(1);
 330         }
 331 }
 332
 333 /**************************************************************************************************/
 334
 335 void AlignCommand::appendAlignFiles(string temp, string filename) {
 336         try{
 337
 338                 //open output file in append mode
 339                 ofstream output;
 340                 openOutputFileAppend(filename, output);
 341
 342                 //open temp file for reading
 343                 ifstream input;
 344                 openInputFile(temp, input);
 345
 346                 string line;
 347                 //read input file and write to output file
 348                 while(input.eof() != true) {
 349                         getline(input, line); //getline removes the newline char
 350                         if (line != "") {
 351                                 output << line << endl;   // Appending back newline char
 352                         }
 353                 }
 354
 355                 input.close();
 356                 output.close();
 357         }
 358         catch(exception& e) {
 359                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 360                 exit(1);
 361         }
 362         catch(...) {
 363                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 364                 exit(1);
 365         }
 366 }
 367
 368 /**************************************************************************************************/
 369
 370 void AlignCommand::appendReportFiles(string temp, string filename) {
 371         try{
 372
 373                 //open output file in append mode
 374                 ofstream output;
 375                 openOutputFileAppend(filename, output);
 376
 377                 //open temp file for reading
 378                 ifstream input;
 379                 openInputFile(temp, input);
 380                 while (!input.eof())    {       char c = input.get(); if (c == 10 || c == 13){  break;  }       } // get header line
 381
 382                 string line;
 383                 //read input file and write to output file
 384                 while(input.eof() != true) {
 385                         getline(input, line); //getline removes the newline char
 386                         if (line != "") {
 387                                 output << line << endl;   // Appending back newline char
 388                         }
 389                 }
 390
 391                 input.close();
 392                 output.close();
 393         }
 394         catch(exception& e) {
 395                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 396                 exit(1);
 397         }
 398         catch(...) {
 399                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 400                 exit(1);
 401         }
 402 }
 403
 404 //**********************************************************************************************************************