]> git.donarmstrong.com Git - mothur.git/blob - aligncommand.cpp
added logfile feature
[mothur.git] / 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") { 
62                                 mothurOut("template is a required parameter for the align.seqs command."); 
63                                 mothurOutEndLine();
64                                 abort = true; 
65                         }
66                         else if (templateFileName == "not open") { abort = true; }      
67                         
68                         candidateFileName = validParameter.validFile(parameters, "candidate", true);
69                         if (candidateFileName == "not found") { 
70                                 mothurOut("candidate is a required parameter for the align.seqs command."); 
71                                 mothurOutEndLine();
72                                 abort = true; 
73                         }
74                         else if (candidateFileName == "not open") { abort = true; }     
75                         
76                         //check for optional parameter and set defaults
77                         // ...at some point should added some additional type checking...
78                         string temp;
79                         temp = validParameter.validFile(parameters, "ksize", false);            if (temp == "not found"){       temp = "8";                             }
80                         convert(temp, kmerSize); 
81                         
82                         temp = validParameter.validFile(parameters, "match", false);            if (temp == "not found"){       temp = "1.0";                   }
83                         convert(temp, match);  
84                         
85                         temp = validParameter.validFile(parameters, "mismatch", false);         if (temp == "not found"){       temp = "-1.0";                  }
86                         convert(temp, misMatch);  
87                         
88                         temp = validParameter.validFile(parameters, "gapopen", false);          if (temp == "not found"){       temp = "-2.0";                  }
89                         convert(temp, gapOpen);  
90                         
91                         temp = validParameter.validFile(parameters, "gapextend", false);        if (temp == "not found"){       temp = "-1.0";                  }
92                         convert(temp, gapExtend); 
93                         
94                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
95                         convert(temp, processors); 
96                         
97                         search = validParameter.validFile(parameters, "search", false);         if (search == "not found"){     search = "kmer";                }
98                         
99                         align = validParameter.validFile(parameters, "align", false);           if (align == "not found"){      align = "needleman";    }
100                 }
101                 
102         }
103         catch(exception& e) {
104                 errorOut(e, "AlignCommand", "AlignCommand");
105                 exit(1);
106         }
107 }
108
109 //**********************************************************************************************************************
110
111 AlignCommand::~AlignCommand(){  
112
113         if (abort == false) {
114                 for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
115                 delete templateDB;
116                 delete alignment;
117         }
118 }
119
120 //**********************************************************************************************************************
121
122 void AlignCommand::help(){
123         try {
124                 mothurOut("The align.seqs command reads a file containing sequences and creates an alignment file and a report file.\n");
125                 mothurOut("The align.seqs command parameters are template, candidate, search, ksize, align, match, mismatch, gapopen and gapextend.\n");
126                 mothurOut("The template and candidate parameters are required.\n");
127                 mothurOut("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");
128                 mothurOut("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");
129                 mothurOut("The ksize parameter allows you to specify the kmer size for finding most similar template to candidate.  The default is 7.\n");
130                 mothurOut("The match parameter allows you to specify the bonus for having the same base. The default is 1.0.\n");
131                 mothurOut("The mistmatch parameter allows you to specify the penalty for having different bases.  The default is -1.0.\n");
132                 mothurOut("The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -1.0.\n");
133                 mothurOut("The gapextend parameter allows you to specify the penalty for extending a gap in an alignment.  The default is -2.0.\n");
134                 mothurOut("The align.seqs command should be in the following format: \n");
135                 mothurOut("align.seqs(template=yourTemplateFile, candidate=yourCandidateFile, align=yourAlignmentMethod, search=yourSearchmethod, ksize=yourKmerSize, match=yourMatchBonus, mismatch=yourMismatchpenalty, gapopen=yourGapopenPenalty, gapextend=yourGapExtendPenalty) \n");
136                 mothurOut("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");
137                 mothurOut("Note: No spaces between parameter labels (i.e. candidate), '=' and parameters (i.e.yourFastaFile).\n\n");
138         }
139         catch(exception& e) {
140                 errorOut(e, "AlignCommand", "help");
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                         mothurOut(search + " is not a valid search option. I will run the command using kmer, ksize=8.");
157                         mothurOutEndLine();
158                         kmerSize = 8;
159                         
160                         templateDB = new KmerDB(templateFileName, kmerSize);
161                 }
162                 
163                 if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000);        }
164                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);                       }
165                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
166                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
167                 else {
168                         mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
169                         mothurOutEndLine();
170                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000);
171                 }
172                 
173                 string alignFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align";
174                 string reportFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align.report";
175                 
176                 int numFastaSeqs = 0;
177                 int start = time(NULL);
178                 
179 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
180                 if(processors == 1){
181                         ifstream inFASTA;
182                         openInputFile(candidateFileName, inFASTA);
183                         numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
184                         inFASTA.close();
185                         
186                         lines.push_back(new linePair(0, numFastaSeqs));
187                         
188                         driver(lines[0], alignFileName, reportFileName);
189                         
190                 }
191                 else{
192                         vector<int> positions;
193                         processIDS.resize(0);
194                         
195                         ifstream inFASTA;
196                         openInputFile(candidateFileName, inFASTA);
197                         
198                         while(!inFASTA.eof()){
199                                 char c = inFASTA.get();
200                                 if(c == '>'){   int pos = inFASTA.tellg(); positions.push_back(pos-1);  }
201                                 while (!inFASTA.eof())  {       c = inFASTA.get(); if (c == 10 || c == 13){     break;  }       } // get rest of line if there's any crap there
202                         }
203                         inFASTA.close();
204                         
205                         numFastaSeqs = positions.size();
206                         
207                         int numSeqsPerProcessor = numFastaSeqs / processors;
208                         
209                         for (int i = 0; i < processors; i++) {
210                                 int startPos = positions[ i * numSeqsPerProcessor ];
211                                 if(i == processors - 1){
212                                         numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
213                                 }
214                                 lines.push_back(new linePair(startPos, numSeqsPerProcessor));
215                         }
216                         createProcesses(alignFileName, reportFileName); 
217                         
218                         rename((alignFileName + toString(processIDS[0]) + ".temp").c_str(), alignFileName.c_str());
219                         rename((reportFileName + toString(processIDS[0]) + ".temp").c_str(), reportFileName.c_str());
220                         
221                         for(int i=1;i<processors;i++){
222                                 appendAlignFiles((alignFileName + toString(processIDS[i]) + ".temp"), alignFileName);
223                                 remove((alignFileName + toString(processIDS[i]) + ".temp").c_str());
224                                 
225                                 appendReportFiles((reportFileName + toString(processIDS[i]) + ".temp"), reportFileName);
226                                 remove((reportFileName + toString(processIDS[i]) + ".temp").c_str());
227                         }
228                         
229                 }
230 #else
231                 ifstream inFASTA;
232                 openInputFile(candidateFileName, inFASTA);
233                 numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
234                 inFASTA.close();
235                 
236                 lines.push_back(new linePair(0, numFastaSeqs));
237                 
238                 driver(lines[0], alignFileName, reportFileName);
239 #endif
240                 
241                 mothurOut("It took " + toString(time(NULL) - start) + " secs to align " + toString(numFastaSeqs) + " sequences.");
242                 mothurOutEndLine();
243                 mothurOutEndLine();
244                 
245                 return 0;
246         }
247         catch(exception& e) {
248                 errorOut(e, "AlignCommand", "execute");
249                 exit(1);
250         }
251 }
252
253 //**********************************************************************************************************************
254
255 int AlignCommand::driver(linePair* line, string alignFName, string reportFName){
256         try {
257                 ofstream alignmentFile;
258                 openOutputFile(alignFName, alignmentFile);
259                 NastReport report(reportFName);
260                 
261                 ifstream inFASTA;
262                 openInputFile(candidateFileName, inFASTA);
263                 inFASTA.seekg(line->start);
264                 
265                 for(int i=0;i<line->numSeqs;i++){
266                         
267                         Sequence* candidateSeq = new Sequence(inFASTA);
268                         report.setCandidate(candidateSeq);
269                         
270                         Sequence temp = templateDB->findClosestSequence(candidateSeq);
271                         Sequence* templateSeq = &temp;
272                         
273                         report.setTemplate(templateSeq);
274                         report.setSearchParameters(search, templateDB->getSearchScore());
275                         
276                         Nast nast(alignment, candidateSeq, templateSeq);
277                         report.setAlignmentParameters(align, alignment);
278                         report.setNastParameters(nast);
279                         
280                         alignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
281                         
282                         report.print();
283                         
284                         delete candidateSeq;            
285                 }
286                 
287                 alignmentFile.close();
288                 inFASTA.close();
289                 
290                 return 1;
291         }
292         catch(exception& e) {
293                 errorOut(e, "AlignCommand", "driver");
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.resize(0);
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.push_back(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 { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); }
317                 }
318                 
319                 //force parent to wait until all the processes are done
320                 for (int i=0;i<processors;i++) { 
321                         int temp = processIDS[i];
322                         wait(&temp);
323                 }
324 #endif          
325         }
326         catch(exception& e) {
327                 errorOut(e, "AlignCommand", "createProcesses");
328                 exit(1);
329         }
330 }
331
332 /**************************************************************************************************/
333
334 void AlignCommand::appendAlignFiles(string temp, string filename) {
335         try{
336                 
337                 ofstream output;
338                 ifstream input;
339                 openOutputFileAppend(filename, output);
340                 openInputFile(temp, input);
341                 
342                 while(char c = input.get()){
343                         if(input.eof())         {       break;                  }
344                         else                            {       output << c;    }
345                 }
346                 
347                 input.close();
348                 output.close();
349         }
350         catch(exception& e) {
351                 errorOut(e, "AlignCommand", "appendAlignFiles");
352                 exit(1);
353         }
354 }
355
356 /**************************************************************************************************/
357
358 void AlignCommand::appendReportFiles(string temp, string filename) {
359         try{
360                 
361                 ofstream output;
362                 ifstream input;
363                 openOutputFileAppend(filename, output);
364                 openInputFile(temp, input);
365
366                 while (!input.eof())    {       char c = input.get(); if (c == 10 || c == 13){  break;  }       } // get header line
367                                 
368                 while(char c = input.get()){
369                         if(input.eof())         {       break;                  }
370                         else                            {       output << c;    }
371                 }
372                 
373                 input.close();
374                 output.close();
375         }
376         catch(exception& e) {
377                 errorOut(e, "AlignCommand", "appendReportFiles");
378                 exit(1);
379         }
380 }
381
382 //**********************************************************************************************************************