]> git.donarmstrong.com Git - mothur.git/blob - chopseqscommand.cpp
Merge remote-tracking branch 'mothur/master'
[mothur.git] / chopseqscommand.cpp
1 /*
2  *  chopseqscommand.cpp
3  *  Mothur
4  *
5  *  Created by westcott on 5/10/10.
6  *  Copyright 2010 Schloss Lab. All rights reserved.
7  *
8  */
9
10 #include "chopseqscommand.h"
11 #include "sequence.hpp"
12
13 //**********************************************************************************************************************
14 vector<string> ChopSeqsCommand::setParameters(){        
15         try {
16                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
17                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
18         CommandParameter pnumbases("numbases", "Number", "", "0", "", "", "",false,true); parameters.push_back(pnumbases);
19                 CommandParameter pcountgaps("countgaps", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcountgaps);
20                 CommandParameter pshort("short", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pshort);
21                 CommandParameter pkeep("keep", "Multiple", "front-back", "front", "", "", "",false,false); parameters.push_back(pkeep);
22                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
23                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
24                 
25                 vector<string> myArray;
26                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
27                 return myArray;
28         }
29         catch(exception& e) {
30                 m->errorOut(e, "ChopSeqsCommand", "setParameters");
31                 exit(1);
32         }
33 }
34 //**********************************************************************************************************************
35 string ChopSeqsCommand::getHelpString(){        
36         try {
37                 string helpString = "";
38                 helpString += "The chop.seqs command reads a fasta file and outputs a .chop.fasta containing the trimmed sequences. Note: If a sequence is completely 'chopped', an accnos file will be created with the names of the sequences removed. \n";
39                 helpString += "The chop.seqs command parameters are fasta, numbases, countgaps and keep. fasta is required unless you have a valid current fasta file. numbases is required.\n";
40                 helpString += "The chop.seqs command should be in the following format: chop.seqs(fasta=yourFasta, numbases=yourNum, keep=yourKeep).\n";
41                 helpString += "The numbases parameter allows you to specify the number of bases you want to keep.\n";
42                 helpString += "The keep parameter allows you to specify whether you want to keep the front or the back of your sequence, default=front.\n";
43                 helpString += "The countgaps parameter allows you to specify whether you want to count gaps as bases, default=false.\n";
44                 helpString += "The short parameter allows you to specify you want to keep sequences that are too short to chop, default=false.\n";
45                 helpString += "The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n";
46         helpString += "For example, if you ran chop.seqs with numbases=200 and short=t, if a sequence had 100 bases mothur would keep the sequence rather than eliminate it.\n";
47                 helpString += "Example chop.seqs(fasta=amazon.fasta, numbases=200, keep=front).\n";
48                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
49                 return helpString;
50         }
51         catch(exception& e) {
52                 m->errorOut(e, "ChopSeqsCommand", "getHelpString");
53                 exit(1);
54         }
55 }
56 //**********************************************************************************************************************
57 string ChopSeqsCommand::getOutputFileNameTag(string type, string inputName=""){ 
58         try {
59         string outputFileName = "";
60                 map<string, vector<string> >::iterator it;
61         
62         //is this a type this command creates
63         it = outputTypes.find(type);
64         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
65         else {
66             if (type == "fasta") {  outputFileName =  "chop.fasta"; }
67             else if (type == "accnos") {  outputFileName =  "chop.accnos"; }
68             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
69         }
70         return outputFileName;
71         }
72         catch(exception& e) {
73                 m->errorOut(e, "ChopSeqsCommand", "getOutputFileNameTag");
74                 exit(1);
75         }
76 }
77 //**********************************************************************************************************************
78 ChopSeqsCommand::ChopSeqsCommand(){     
79         try {
80                 abort = true; calledHelp = true; 
81                 setParameters();
82                 vector<string> tempOutNames;
83                 outputTypes["fasta"] = tempOutNames;
84                 outputTypes["accnos"] = tempOutNames;
85         }
86         catch(exception& e) {
87                 m->errorOut(e, "ChopSeqsCommand", "ChopSeqsCommand");
88                 exit(1);
89         }
90 }
91 //**********************************************************************************************************************
92 ChopSeqsCommand::ChopSeqsCommand(string option)  {
93         try {
94                 abort = false; calledHelp = false;   
95                 
96                 //allow user to run help
97                 if(option == "help") { help(); abort = true; calledHelp = true; }
98                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
99                 
100                 else {
101                         vector<string> myArray = setParameters();
102                         
103                         OptionParser parser(option);
104                         map<string,string> parameters = parser.getParameters();
105                         
106                         ValidParameters validParameter;
107                         map<string,string>::iterator it;
108                         
109                         //check to make sure all parameters are valid for command
110                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) { 
111                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
112                         }
113                         
114                         //initialize outputTypes
115                         vector<string> tempOutNames;
116                         outputTypes["fasta"] = tempOutNames;
117                         outputTypes["accnos"] = tempOutNames;
118                 
119                         //if the user changes the input directory command factory will send this info to us in the output parameter 
120                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
121                         if (inputDir == "not found"){   inputDir = "";          }
122                         else {
123                                 string path;
124                                 it = parameters.find("fasta");
125                                 //user has given a template file
126                                 if(it != parameters.end()){ 
127                                         path = m->hasPath(it->second);
128                                         //if the user has not given a path then, add inputdir. else leave path alone.
129                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
130                                 }
131                         }
132
133                         //check for required parameters
134                         fastafile = validParameter.validFile(parameters, "fasta", true);
135                         if (fastafile == "not open") { abort = true; }
136                         else if (fastafile == "not found") {                            //if there is a current fasta file, use it
137                                 fastafile = m->getFastaFile(); 
138                                 if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
139                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
140                         }else { m->setFastaFile(fastafile); }   
141                         
142                         //if the user changes the output directory command factory will send this info to us in the output parameter 
143                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(fastafile);      }
144                         
145                         string temp = validParameter.validFile(parameters, "numbases", false);  if (temp == "not found") { temp = "0"; } 
146                         m->mothurConvert(temp, numbases);   
147                         
148             temp = validParameter.validFile(parameters, "processors", false);   if (temp == "not found"){       temp = m->getProcessors();      }
149                         m->setProcessors(temp);
150                         m->mothurConvert(temp, processors);
151             
152                         temp = validParameter.validFile(parameters, "countgaps", false);        if (temp == "not found") { temp = "f"; } 
153                         countGaps = m->isTrue(temp);  
154                         
155                         temp = validParameter.validFile(parameters, "short", false);    if (temp == "not found") { temp = "f"; } 
156                         Short = m->isTrue(temp);   
157                 
158                         keep = validParameter.validFile(parameters, "keep", false);             if (keep == "not found") { keep = "front"; } 
159                                 
160                         if (numbases == 0)  { m->mothurOut("You must provide the number of bases you want to keep for the chops.seqs command."); m->mothurOutEndLine(); abort = true;  }
161                 }
162
163         }
164         catch(exception& e) {
165                 m->errorOut(e, "ChopSeqsCommand", "ChopSeqsCommand");
166                 exit(1);
167         }
168 }
169 //**********************************************************************************************************************
170
171 int ChopSeqsCommand::execute(){
172         try {
173                 
174                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
175                 
176                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + getOutputFileNameTag("fasta");
177                 string outputFileNameAccnos = outputDir + m->getRootName(m->getSimpleName(fastafile)) + getOutputFileNameTag("accnos");
178         
179         
180         vector<unsigned long long> positions; 
181         vector<linePair> lines;
182 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
183         positions = m->divideFile(fastafile, processors);
184         for (int i = 0; i < (positions.size()-1); i++) {        lines.push_back(linePair(positions[i], positions[(i+1)]));      }
185 #else
186         int numSeqs = 0;
187         positions = m->setFilePosFasta(fastafile, numSeqs); 
188         if (positions.size() < processors) { processors = positions.size(); }
189                 
190         //figure out how many sequences you have to process
191         int numSeqsPerProcessor = numSeqs / processors;
192         for (int i = 0; i < processors; i++) {
193             int startIndex =  i * numSeqsPerProcessor;
194             if(i == (processors - 1)){  numSeqsPerProcessor = numSeqs - i * numSeqsPerProcessor;        }
195             lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
196         }
197 #endif
198         
199         bool wroteAccnos = false;
200         if(processors == 1) {   wroteAccnos = driver(lines[0], fastafile, outputFileName, outputFileNameAccnos);        }
201         else                {   wroteAccnos = createProcesses(lines, fastafile, outputFileName, outputFileNameAccnos);  }
202         
203         if (m->control_pressed) {  return 0; }
204                 
205                 m->mothurOutEndLine();
206                 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
207                 m->mothurOut(outputFileName); m->mothurOutEndLine();    outputNames.push_back(outputFileName); outputTypes["fasta"].push_back(outputFileName);
208                 
209                 if (wroteAccnos) { m->mothurOut(outputFileNameAccnos); m->mothurOutEndLine(); outputNames.push_back(outputFileNameAccnos); outputTypes["accnos"].push_back(outputFileNameAccnos); }
210                 else {  m->mothurRemove(outputFileNameAccnos);  }
211                 
212                 m->mothurOutEndLine();
213                 
214                 //set fasta file as new current fastafile
215                 string current = "";
216                 itTypes = outputTypes.find("fasta");
217                 if (itTypes != outputTypes.end()) {
218                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
219                 }
220                 
221                 if (wroteAccnos) { //set accnos file as new current accnosfile
222                         itTypes = outputTypes.find("accnos");
223                         if (itTypes != outputTypes.end()) {
224                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
225                         }
226                 }
227                 
228                 
229                 return 0;               
230         }
231
232         catch(exception& e) {
233                 m->errorOut(e, "ChopSeqsCommand", "execute");
234                 exit(1);
235         }
236 }
237 /**************************************************************************************************/
238 bool ChopSeqsCommand::createProcesses(vector<linePair> lines, string filename, string outFasta, string outAccnos) {
239         try {
240                 int process = 1;
241                 bool wroteAccnos = false;
242                 vector<int> processIDS;
243         vector<string> nonBlankAccnosFiles;
244                 
245 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
246                 
247                 //loop through and create all the processes you want
248                 while (process != processors) {
249                         int pid = fork();
250                         
251                         if (pid > 0) {
252                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
253                                 process++;
254                         }else if (pid == 0){
255                                 wroteAccnos = driver(lines[process], filename, outFasta + toString(getpid()) + ".temp", outAccnos + toString(getpid()) + ".temp");
256                                 
257                                 //pass numSeqs to parent
258                                 ofstream out;
259                                 string tempFile = fastafile + toString(getpid()) + ".bool.temp";
260                                 m->openOutputFile(tempFile, out);
261                                 out << wroteAccnos << endl;                             
262                                 out.close();
263                                 
264                                 exit(0);
265                         }else { 
266                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
267                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
268                                 exit(0);
269                         }
270                 }
271                 
272                 //do your part
273                 wroteAccnos = driver(lines[0], filename, outFasta, outAccnos);
274         
275                 //force parent to wait until all the processes are done
276                 for (int i=0;i<processIDS.size();i++) { 
277                         int temp = processIDS[i];
278                         wait(&temp);
279                 }
280                 
281         
282                 if (wroteAccnos) { nonBlankAccnosFiles.push_back(outAccnos); }
283                 else { m->mothurRemove(outAccnos); } //remove so other files can be renamed to it
284         
285                 //parent reads in and combine Filter info
286                 for (int i = 0; i < processIDS.size(); i++) {
287                         string tempFilename = fastafile + toString(processIDS[i]) + ".bool.temp";
288                         ifstream in;
289                         m->openInputFile(tempFilename, in);
290                         
291                         bool temp;
292                         in >> temp; m->gobble(in); 
293             if (temp) { wroteAccnos = temp; nonBlankAccnosFiles.push_back(outAccnos + toString(processIDS[i]) + ".temp");  }
294                         else { m->mothurRemove((outAccnos + toString(processIDS[i]) + ".temp"));  }
295             
296                         in.close();
297                         m->mothurRemove(tempFilename);
298                 }
299 #else
300                 //////////////////////////////////////////////////////////////////////////////////////////////////////
301                 //Windows version shared memory, so be careful when passing variables through the seqSumData struct. 
302                 //Above fork() will clone, so memory is separate, but that's not the case with windows, 
303                 //Taking advantage of shared memory to allow both threads to add info to vectors.
304                 //////////////////////////////////////////////////////////////////////////////////////////////////////
305                 
306                 vector<chopData*> pDataArray; 
307                 DWORD   dwThreadIdArray[processors-1];
308                 HANDLE  hThreadArray[processors-1]; 
309                 
310                 //Create processor worker threads.
311                 for( int i=0; i<processors-1; i++ ){
312             
313             string extension = "";
314             if (i != 0) { extension = toString(i) + ".temp"; processIDS.push_back(i); }
315                         // Allocate memory for thread data.
316                         chopData* tempChop = new chopData(filename, (outFasta+extension), (outAccnos+extension), m, lines[i].start, lines[i].end, keep, countGaps, numbases, Short);
317                         pDataArray.push_back(tempChop);
318                         
319                         //MyChopThreadFunction is in header. It must be global or static to work with the threads.
320                         //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
321                         hThreadArray[i] = CreateThread(NULL, 0, MyChopThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);   
322                 }
323                 
324         //do your part
325                 wroteAccnos = driver(lines[processors-1], filename, (outFasta + toString(processors-1) + ".temp"), (outAccnos + toString(processors-1) + ".temp"));
326         processIDS.push_back(processors-1);
327         
328                 //Wait until all threads have terminated.
329                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
330                 
331         if (wroteAccnos) { nonBlankAccnosFiles.push_back(outAccnos); }
332                 else { m->mothurRemove(outAccnos); } //remove so other files can be renamed to it
333
334                 //Close all thread handles and free memory allocations.
335                 for(int i=0; i < pDataArray.size(); i++){
336             if (pDataArray[i]->wroteAccnos) { wroteAccnos = pDataArray[i]->wroteAccnos; nonBlankAccnosFiles.push_back(outAccnos + toString(processIDS[i]) + ".temp");  }
337                         else { m->mothurRemove((outAccnos + toString(processIDS[i]) + ".temp"));  }
338                         CloseHandle(hThreadArray[i]);
339                         delete pDataArray[i];
340                 }
341 #endif          
342                 
343                 for (int i = 0; i < processIDS.size(); i++) {
344                         m->appendFiles((outFasta + toString(processIDS[i]) + ".temp"), outFasta);
345                         m->mothurRemove((outFasta + toString(processIDS[i]) + ".temp"));
346                 }
347                 
348         if (nonBlankAccnosFiles.size() != 0) { 
349                         m->renameFile(nonBlankAccnosFiles[0], outAccnos);
350                         
351                         for (int h=1; h < nonBlankAccnosFiles.size(); h++) {
352                                 m->appendFiles(nonBlankAccnosFiles[h], outAccnos);
353                                 m->mothurRemove(nonBlankAccnosFiles[h]);
354                         }
355                 }else { //recreate the accnosfile if needed
356                         ofstream out;
357                         m->openOutputFile(outAccnos, out);
358                         out.close();
359                 }
360
361                 return wroteAccnos;
362         }
363         catch(exception& e) {
364                 m->errorOut(e, "ChopSeqsCommand", "createProcesses");
365                 exit(1);
366         }
367 }
368 /**************************************************************************************/
369 bool ChopSeqsCommand::driver(linePair filePos, string filename, string outFasta, string outAccnos) {    
370         try {
371                 
372                 ofstream out;
373                 m->openOutputFile(outFasta, out);
374         
375         ofstream outAcc;
376                 m->openOutputFile(outAccnos, outAcc);
377         
378                 ifstream in;
379                 m->openInputFile(filename, in);
380         
381                 in.seekg(filePos.start);
382         
383                 bool done = false;
384         bool wroteAccnos = false;
385                 int count = 0;
386         
387                 while (!done) {
388             
389                         if (m->control_pressed) { in.close(); out.close(); return 1; }
390             
391                         Sequence seq(in); m->gobble(in);
392                         
393                         if (m->control_pressed) {  in.close(); out.close(); outAcc.close(); m->mothurRemove(outFasta); m->mothurRemove(outAccnos); return 0;  }
394                         
395                         if (seq.getName() != "") {
396                                 string newSeqString = getChopped(seq);
397                                 
398                                 //output trimmed sequence
399                                 if (newSeqString != "") {
400                                         out << ">" << seq.getName() << endl << newSeqString << endl;
401                                 }else{
402                                         outAcc << seq.getName() << endl;
403                                         wroteAccnos = true;
404                                 }
405                 count++;
406                         }
407                         
408 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
409             unsigned long long pos = in.tellg();
410             if ((pos == -1) || (pos >= filePos.end)) { break; }
411 #else
412             if (in.eof()) { break; }
413 #endif
414             //report progress
415                         if((count) % 1000 == 0){        m->mothurOut(toString(count)); m->mothurOutEndLine();           }
416                         
417                 }
418                 //report progress
419                 if((count) % 1000 != 0){        m->mothurOut(toString(count)); m->mothurOutEndLine();           }
420
421                 
422                 in.close();
423         out.close();
424         outAcc.close();
425                 
426                 return wroteAccnos;
427         }
428         catch(exception& e) {
429                 m->errorOut(e, "ChopSeqsCommand", "driver");
430                 exit(1);
431         }
432 }
433 //**********************************************************************************************************************
434 string ChopSeqsCommand::getChopped(Sequence seq) {
435         try {
436                 string temp = seq.getAligned();
437                 string tempUnaligned = seq.getUnaligned();
438                 
439                 if (countGaps) {
440                         //if needed trim sequence
441                         if (keep == "front") {//you want to keep the beginning
442                                 int tempLength = temp.length();
443
444                                 if (tempLength > numbases) { //you have enough bases to remove some
445                                 
446                                         int stopSpot = 0;
447                                         int numBasesCounted = 0;
448                                         
449                                         for (int i = 0; i < temp.length(); i++) {
450                                                 //eliminate N's
451                                                 if (toupper(temp[i]) == 'N') { temp[i] = '.'; }
452                                                 
453                                                 numBasesCounted++; 
454                                                 
455                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
456                                         }
457                                         
458                                         if (stopSpot == 0) { temp = ""; }
459                                         else {  temp = temp.substr(0, stopSpot+1);  }
460                                                         
461                                 }else { 
462                                         if (!Short) { temp = ""; } //sequence too short
463                                 }
464                         }else { //you are keeping the back
465                                 int tempLength = temp.length();
466                                 if (tempLength > numbases) { //you have enough bases to remove some
467                                         
468                                         int stopSpot = 0;
469                                         int numBasesCounted = 0;
470                                         
471                                         for (int i = (temp.length()-1); i >= 0; i--) {
472                                                 //eliminate N's
473                                                 if (toupper(temp[i]) == 'N') { temp[i] = '.'; }
474                                                 
475                                                 numBasesCounted++; 
476
477                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
478                                         }
479                                 
480                                         if (stopSpot == 0) { temp = ""; }
481                                         else {  temp = temp.substr(stopSpot+1);  }
482                                 }else { 
483                                         if (!Short) { temp = ""; } //sequence too short
484                                 }
485                         }
486
487                 }else{
488                                 
489                         //if needed trim sequence
490                         if (keep == "front") {//you want to keep the beginning
491                                 int tempLength = tempUnaligned.length();
492
493                                 if (tempLength > numbases) { //you have enough bases to remove some
494                                         
495                                         int stopSpot = 0;
496                                         int numBasesCounted = 0;
497                                         
498                                         for (int i = 0; i < temp.length(); i++) {
499                                                 //eliminate N's
500                                                 if (toupper(temp[i]) == 'N') { 
501                                                         temp[i] = '.'; 
502                                                         tempLength--;
503                                                         if (tempLength < numbases) { stopSpot = 0; break; }
504                                                 }
505                                                 
506                                                 if(isalpha(temp[i])) { numBasesCounted++; }
507                                                 
508                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
509                                         }
510                                         
511                                         if (stopSpot == 0) { temp = ""; }
512                                         else {  temp = temp.substr(0, stopSpot+1);  }
513                                                         
514                                 }else { 
515                                         if (!Short) { temp = ""; } //sequence too short
516                                 }                               
517                         }else { //you are keeping the back
518                                 int tempLength = tempUnaligned.length();
519                                 if (tempLength > numbases) { //you have enough bases to remove some
520                                         
521                                         int stopSpot = 0;
522                                         int numBasesCounted = 0;
523                                         
524                                         for (int i = (temp.length()-1); i >= 0; i--) {
525                                                 //eliminate N's
526                                                 if (toupper(temp[i]) == 'N') { 
527                                                         temp[i] = '.'; 
528                                                         tempLength--;
529                                                         if (tempLength < numbases) { stopSpot = 0; break; }
530                                                 }
531                                                 
532                                                 if(isalpha(temp[i])) { numBasesCounted++; }
533
534                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
535                                         }
536                                 
537                                         if (stopSpot == 0) { temp = ""; }
538                                         else {  temp = temp.substr(stopSpot);  }
539                                 }else { 
540                                         if (!Short) { temp = ""; } //sequence too short
541                                 }
542                         }
543                 }
544                 
545                 return temp;
546         }
547         catch(exception& e) {
548                 m->errorOut(e, "ChopSeqsCommand", "getChopped");
549                 exit(1);
550         }
551 }
552 //**********************************************************************************************************************
553
554