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 = "-1.0";                  }
  81                         convert(temp, gapOpen);
  82
  83                         temp = validParameter.validFile(parameters, "gapextend", false);        if (temp == "not found"){       temp = "-2.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 fasta, candidate, search, ksize, align, match, mismatch, gapopen and gapextend.  " << "\n";
 118                 cout << "The fasta 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(fasta=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, fasta=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. fasta), '=' 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                 srand( (unsigned)time( NULL ) );  //needed to assign names to temporary files
 149
 150
 151                 if(search == "kmer")                    {       templateDB = new KmerDB(templateFileName, kmerSize);    }
 152                 else if(search == "suffix")             {       templateDB = new SuffixDB(templateFileName);                    }
 153                 else if(search == "blast")              {       templateDB = new BlastDB(templateFileName, gapOpen, gapExtend, match, misMatch);        }
 154                 else {
 155                         cout << search << " is not a valid search option. I will run the command using kmer, ksize=8." << endl; kmerSize = 8;
 156                         templateDB = new KmerDB(templateFileName, kmerSize);
 157                 }
 158
 159                 if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
 160                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
 161                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 162                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 163                 else {
 164                         cout << align << " is not a valid alignment option. I will run the command using needleman." << endl;
 165                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);
 166                 }
 167
 168                 string alignFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align";
 169                 string reportFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align.report";
 170
 171                 int numFastaSeqs = 0;
 172                 int start = time(NULL);
 173
 174 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 175                 if(processors == 1){
 176                         ifstream inFASTA;
 177                         openInputFile(candidateFileName, inFASTA);
 178                         numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 179                         inFASTA.close();
 180
 181                         lines.push_back(new linePair(0, numFastaSeqs));
 182
 183                         driver(lines[0], alignFileName, reportFileName);
 184                 }
 185                 else{
 186                         vector<int> positions;
 187
 188                         ifstream inFASTA;
 189                         openInputFile(candidateFileName, inFASTA);
 190
 191                         while(!inFASTA.eof()){
 192                                 char c = inFASTA.get();
 193                                 if(c == '>'){   positions.push_back(inFASTA.tellg());   }
 194                                 while (!inFASTA.eof())  {       c = inFASTA.get(); if (c == 10 || c == 13){     break;  }       } // get rest of line if there's any crap there
 195                         }
 196                         inFASTA.close();
 197
 198                         numFastaSeqs = positions.size();
 199
 200                         int numSeqsPerProcessor = numFastaSeqs / processors;
 201
 202                         for (int i = 0; i < processors; i++) {
 203                                 int startPos = positions[ i * numSeqsPerProcessor ];
 204                                 if(i == processors - 1){
 205                                         numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
 206                                 }
 207                                 lines.push_back(new linePair(startPos, numSeqsPerProcessor));
 208                         }
 209                         createProcesses(alignFileName, reportFileName);
 210
 211
 212                         //append and remove temp files
 213                         map<int, int>::iterator it = processIDS.begin();
 214                         rename((alignFileName + toString(it->second) + ".temp").c_str(), alignFileName.c_str());
 215                         rename((reportFileName + toString(it->second) + ".temp").c_str(), alignFileName.c_str());
 216                         it++;
 217
 218                         for (; it != processIDS.end(); it++) {
 219                                 appendAlignFiles((alignFileName + toString(it->second) + ".temp"), alignFileName);
 220                                 remove((alignFileName + toString(it->second) + ".temp").c_str());
 221
 222                                 appendReportFiles((reportFileName + toString(it->second) + ".temp"), reportFileName);
 223                                 remove((reportFileName + toString(it->second) + ".temp").c_str());
 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* templateSeq = templateDB->findClosestSequence(candidateSeq);
 272                         report.setTemplate(templateSeq);
 273                         report.setSearchParameters(search, templateDB->getSearchScore());
 274
 275                         Nast nast(alignment, candidateSeq, templateSeq);
 276                         report.setAlignmentParameters(align, alignment);
 277                         report.setNastParameters(nast);
 278
 279                         alignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
 280
 281                         report.print();
 282
 283                         delete candidateSeq;
 284                 }
 285
 286                 return 1;
 287         }
 288         catch(exception& e) {
 289                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 290                 exit(1);
 291         }
 292         catch(...) {
 293                 cout << "An unknown error has occurred in the AlignCommand class function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 294                 exit(1);
 295         }
 296 }
 297
 298 /**************************************************************************************************/
 299
 300 void AlignCommand::createProcesses(string alignFileName, string reportFileName) {
 301         try {
 302 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 303                 int process = 0;
 304                 processIDS.clear();
 305
 306                 //loop through and create all the processes you want
 307                 while (process != processors) {
 308                         int pid = fork();
 309
 310                         if (pid > 0) {
 311                                 processIDS[lines[process]->numSeqs] = pid;  //create map from line number to pid so you can append files in correct order later
 312                                 process++;
 313                         }else if (pid == 0){
 314                                 driver(lines[process], alignFileName + toString(getpid()) + ".temp", reportFileName + toString(getpid()) + ".temp");
 315                                 exit(0);
 316                         }else { cout << "unable to spawn the necessary processes." << endl; exit(0); }
 317                 }
 318
 319                 //force parent to wait until all the processes are done
 320                 for (map<int, int>::iterator it = processIDS.begin(); it != processIDS.end(); it++) {
 321                         int temp = it->second;
 322                         wait(&temp);
 323                 }
 324 #endif
 325         }
 326         catch(exception& e) {
 327                 cout << "Standard Error: " << e.what() << " has occurred in the AlignCommand class Function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 328                 exit(1);
 329         }
 330         catch(...) {
 331                 cout << "An unknown error has occurred in the AlignCommand class function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 332                 exit(1);
 333         }
 334 }
 335
 336 /**************************************************************************************************/
 337
 338 void AlignCommand::appendAlignFiles(string temp, string filename) {
 339         try{
 340
 341                 //open output file in append mode
 342                 ofstream output;
 343                 openOutputFileAppend(filename, output);
 344
 345                 //open temp file for reading
 346                 ifstream input;
 347                 openInputFile(temp, input);
 348
 349                 string line;
 350                 //read input file and write to output file
 351                 while(input.eof() != true) {
 352                         getline(input, line); //getline removes the newline char
 353                         if (line != "") {
 354                                 output << line << endl;   // Appending back newline char
 355                         }
 356                 }
 357
 358                 input.close();
 359                 output.close();
 360         }
 361         catch(exception& e) {
 362                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 363                 exit(1);
 364         }
 365         catch(...) {
 366                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 367                 exit(1);
 368         }
 369 }
 370
 371 /**************************************************************************************************/
 372
 373 void AlignCommand::appendReportFiles(string temp, string filename) {
 374         try{
 375
 376                 //open output file in append mode
 377                 ofstream output;
 378                 openOutputFileAppend(filename, output);
 379
 380                 //open temp file for reading
 381                 ifstream input;
 382                 openInputFile(temp, input);
 383                 while (!input.eof())    {       char c = input.get(); if (c == 10 || c == 13){  break;  }       } // get header line
 384
 385                 string line;
 386                 //read input file and write to output file
 387                 while(input.eof() != true) {
 388                         getline(input, line); //getline removes the newline char
 389                         if (line != "") {
 390                                 output << line << endl;   // Appending back newline char
 391                         }
 392                 }
 393
 394                 input.close();
 395                 output.close();
 396         }
 397         catch(exception& e) {
 398                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 399                 exit(1);
 400         }
 401         catch(...) {
 402                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 403                 exit(1);
 404         }
 405 }
 406
 407 //**********************************************************************************************************************