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
 109         if (abort == false) {
 110                 for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 111                 delete templateDB;
 112                 delete alignment;
 113         }
 114 }
 115
 116 //**********************************************************************************************************************
 117
 118 void AlignCommand::help(){
 119         try {
 120                 cout << "The align.seqs command reads a file containing sequences and creates an alignment file and a report file." << "\n";
 121                 cout << "The align.seqs command parameters are template, candidate, search, ksize, align, match, mismatch, gapopen and gapextend.  " << "\n";
 122                 cout << "The template and candidate parameters are required." << "\n";
 123                 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";
 124                 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";
 125                 cout << "The ksize parameter allows you to specify the kmer size for finding most similar template to candidate.  The default is 7." << "\n";
 126                 cout << "The match parameter allows you to specify the bonus for having the same base. The default is 1.0." << "\n";
 127                 cout << "The mistmatch parameter allows you to specify the penalty for having different bases.  The default is -1.0." << "\n";
 128                 cout << "The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -1.0." << "\n";
 129                 cout << "The gapextend parameter allows you to specify the penalty for extending a gap in an alignment.  The default is -2.0." << "\n";
 130                 cout << "The align.seqs command should be in the following format: " << "\n";
 131                 cout << "align.seqs(template=yourTemplateFile, candidate=yourCandidateFile, align=yourAlignmentMethod, search=yourSearchmethod, ksize=yourKmerSize, match=yourMatchBonus, mismatch=yourMismatchpenalty, gapopen=yourGapopenPenalty, gapextend=yourGapExtendPenalty) " << "\n";
 132                 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";
 133                 cout << "Note: No spaces between parameter labels (i.e. candidate), '=' and parameters (i.e.yourFastaFile)." << "\n" << "\n";
 134         }
 135         catch(exception& e) {
 136                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 137                 exit(1);
 138         }
 139         catch(...) {
 140                 cout << "An unknown error has occurred in the AlignCommand class function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 141                 exit(1);
 142         }
 143 }
 144
 145
 146 //**********************************************************************************************************************
 147
 148 int AlignCommand::execute(){
 149         try {
 150                 if (abort == true) {    return 0;       }
 151
 152                 if(search == "kmer")                    {       templateDB = new KmerDB(templateFileName, kmerSize);    }
 153                 else if(search == "suffix")             {       templateDB = new SuffixDB(templateFileName);                    }
 154                 else if(search == "blast")              {       templateDB = new BlastDB(templateFileName, gapOpen, gapExtend, match, misMatch);        }
 155                 else {
 156                         cout << search << " is not a valid search option. I will run the command using kmer, ksize=8." << endl; kmerSize = 8;
 157                         templateDB = new KmerDB(templateFileName, kmerSize);
 158                 }
 159
 160                 if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
 161                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
 162                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 163                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 164                 else {
 165                         cout << align << " is not a valid alignment option. I will run the command using needleman." << endl;
 166                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);
 167                 }
 168
 169                 string alignFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align";
 170                 string reportFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align.report";
 171
 172                 int numFastaSeqs = 0;
 173                 int start = time(NULL);
 174
 175 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 176                 if(processors == 1){
 177                         ifstream inFASTA;
 178                         openInputFile(candidateFileName, inFASTA);
 179                         numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 180                         inFASTA.close();
 181
 182                         lines.push_back(new linePair(0, numFastaSeqs));
 183
 184                         driver(lines[0], alignFileName, reportFileName);
 185
 186                 }
 187                 else{
 188                         vector<int> positions;
 189                         processIDS.resize(0);
 190
 191                         ifstream inFASTA;
 192                         openInputFile(candidateFileName, inFASTA);
 193
 194                         while(!inFASTA.eof()){
 195                                 char c = inFASTA.get();
 196                                 if(c == '>'){   positions.push_back(inFASTA.tellg());   }
 197                                 while (!inFASTA.eof())  {       c = inFASTA.get(); if (c == 10 || c == 13){     break;  }       } // get rest of line if there's any crap there
 198                         }
 199                         inFASTA.close();
 200
 201                         numFastaSeqs = positions.size();
 202
 203                         int numSeqsPerProcessor = numFastaSeqs / processors;
 204
 205                         for (int i = 0; i < processors; i++) {
 206                                 int startPos = positions[ i * numSeqsPerProcessor ];
 207                                 if(i == processors - 1){
 208                                         numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
 209                                 }
 210                                 lines.push_back(new linePair(startPos, numSeqsPerProcessor));
 211                         }
 212                         createProcesses(alignFileName, reportFileName);
 213
 214                         rename((alignFileName + toString(processIDS[0]) + ".temp").c_str(), alignFileName.c_str());
 215                         rename((reportFileName + toString(processIDS[0]) + ".temp").c_str(), reportFileName.c_str());
 216
 217                         for(int i=1;i<processors;i++){
 218                                 appendAlignFiles((alignFileName + toString(processIDS[i]) + ".temp"), alignFileName);
 219                                 remove((alignFileName + toString(processIDS[i]) + ".temp").c_str());
 220
 221                                 appendReportFiles((reportFileName + toString(processIDS[i]) + ".temp"), reportFileName);
 222                                 remove((reportFileName + toString(processIDS[i]) + ".temp").c_str());
 223                         }
 224
 225                 }
 226 #else
 227                 ifstream inFASTA;
 228                 openInputFile(candidateFileName, inFASTA);
 229                 numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 230                 inFASTA.close();
 231
 232                 lines.push_back(new linePair(0, numFastaSeqs));
 233
 234                 driver(lines[0], alignFileName, reportFileName);
 235 #endif
 236
 237                 cout << "It took " << time(NULL) - start << " secs to align " << numFastaSeqs << " sequences" << endl;
 238                 cout << endl;
 239
 240
 241
 242                 return 0;
 243         }
 244         catch(exception& e) {
 245                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 246                 exit(1);
 247         }
 248         catch(...) {
 249                 cout << "An unknown error has occurred in the AlignCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 250                 exit(1);
 251         }
 252 }
 253
 254 //**********************************************************************************************************************
 255
 256 int AlignCommand::driver(linePair* line, string alignFName, string reportFName){
 257         try {
 258                 ofstream alignmentFile;
 259                 openOutputFile(alignFName, alignmentFile);
 260                 NastReport report(reportFName);
 261
 262                 ifstream inFASTA;
 263                 openInputFile(candidateFileName, inFASTA);
 264                 inFASTA.seekg(line->start);
 265
 266                 for(int i=0;i<line->numSeqs;i++){
 267
 268                         Sequence* candidateSeq = new Sequence(inFASTA);
 269                         report.setCandidate(candidateSeq);
 270
 271                         Sequence temp = templateDB->findClosestSequence(candidateSeq);
 272                         Sequence* templateSeq = &temp;
 273
 274                         report.setTemplate(templateSeq);
 275                         report.setSearchParameters(search, templateDB->getSearchScore());
 276
 277                         Nast nast(alignment, candidateSeq, templateSeq);
 278                         report.setAlignmentParameters(align, alignment);
 279                         report.setNastParameters(nast);
 280
 281                         alignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
 282
 283                         report.print();
 284
 285                         delete candidateSeq;
 286                 }
 287
 288                 alignmentFile.close();
 289                 inFASTA.close();
 290
 291                 return 1;
 292         }
 293         catch(exception& e) {
 294                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 295                 exit(1);
 296         }
 297         catch(...) {
 298                 cout << "An unknown error has occurred in the AlignCommand class function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 299                 exit(1);
 300         }
 301 }
 302
 303 /**************************************************************************************************/
 304
 305 void AlignCommand::createProcesses(string alignFileName, string reportFileName) {
 306         try {
 307 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 308                 int process = 0;
 309                 //              processIDS.resize(0);
 310
 311                 //loop through and create all the processes you want
 312                 while (process != processors) {
 313                         int pid = fork();
 314
 315                         if (pid > 0) {
 316                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 317                                 process++;
 318                         }else if (pid == 0){
 319                                 driver(lines[process], alignFileName + toString(getpid()) + ".temp", reportFileName + toString(getpid()) + ".temp");
 320                                 exit(0);
 321                         }else { cout << "unable to spawn the necessary processes." << endl; exit(0); }
 322                 }
 323
 324                 //force parent to wait until all the processes are done
 325                 for (int i=0;i<processors;i++) {
 326                         int temp = processIDS[i];
 327                         wait(&temp);
 328                 }
 329 #endif
 330         }
 331         catch(exception& e) {
 332                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 333                 exit(1);
 334         }
 335         catch(...) {
 336                 cout << "An unknown error has occurred in the AlignCommand class function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 337                 exit(1);
 338         }
 339 }
 340
 341 /**************************************************************************************************/
 342
 343 void AlignCommand::appendAlignFiles(string temp, string filename) {
 344         try{
 345
 346                 //open output file in append mode
 347                 ofstream output;
 348                 openOutputFileAppend(filename, output);
 349
 350                 //open temp file for reading
 351                 ifstream input;
 352                 openInputFile(temp, input);
 353
 354                 string line;
 355                 //read input file and write to output file
 356                 while(input.eof() != true) {
 357                         getline(input, line); //getline removes the newline char
 358                         if (line != "") {
 359                                 output << line << endl;   // Appending back newline char
 360                         }
 361                 }
 362
 363                 input.close();
 364                 output.close();
 365         }
 366         catch(exception& e) {
 367                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 368                 exit(1);
 369         }
 370         catch(...) {
 371                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 372                 exit(1);
 373         }
 374 }
 375
 376 /**************************************************************************************************/
 377
 378 void AlignCommand::appendReportFiles(string temp, string filename) {
 379         try{
 380
 381                 //open output file in append mode
 382                 ofstream output;
 383                 openOutputFileAppend(filename, output);
 384
 385                 //open temp file for reading
 386                 ifstream input;
 387                 openInputFile(temp, input);
 388                 while (!input.eof())    {       char c = input.get(); if (c == 10 || c == 13){  break;  }       } // get header line
 389
 390                 string line;
 391                 //read input file and write to output file
 392                 while(input.eof() != true) {
 393                         getline(input, line); //getline removes the newline char
 394                         if (line != "") {
 395                                 output << line << endl;   // Appending back newline char
 396                         }
 397                 }
 398
 399                 input.close();
 400                 output.close();
 401         }
 402         catch(exception& e) {
 403                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 404                 exit(1);
 405         }
 406         catch(...) {
 407                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 408                 exit(1);
 409         }
 410 }
 411
 412 //**********************************************************************************************************************