]> git.donarmstrong.com Git - mothur.git/blob - aligncommand.h
testing public repo change
[mothur.git] / aligncommand.h
1 #ifndef ALIGNCOMMAND_H
2 #define ALIGNCOMMAND_H
3
4 /*
5  *  aligncommand.h
6  *  Mothur
7  *
8  *  Created by Sarah Westcott on 5/15/09.
9  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10  *
11  */
12
13 #include "mothur.h"
14 #include "command.hpp"
15 #include "database.hpp"
16 #include "alignment.hpp"
17 #include "alignmentdb.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 "nast.hpp"
26 #include "nastreport.hpp"
27
28
29 class AlignCommand : public Command {
30         
31 public:
32         AlignCommand(string);   
33         AlignCommand();
34         ~AlignCommand();
35         
36         vector<string> setParameters();
37         string getCommandName()                 { return "align.seqs";                  }
38         string getCommandCategory()             { return "Sequence Processing"; }
39         string getHelpString(); 
40         string getCitation() { return "DeSantis TZ, Jr., Hugenholtz P, Keller K, Brodie EL, Larsen N, Piceno YM, Phan R, Andersen GL (2006). NAST: a multiple sequence alignment server for comparative analysis of 16S rRNA genes. Nucleic Acids Res 34: W394-9.\nSchloss PD (2009). A high-throughput DNA sequence aligner for microbial ecology studies. PLoS ONE 4: e8230.\nSchloss PD (2010). The effects of alignment quality, distance calculation method, sequence filtering, and region on the analysis of 16S rRNA gene-based studies. PLoS Comput Biol 6: e1000844.\nhttp://www.mothur.org/wiki/Align.seqs http://www.mothur.org/wiki/Align.seqs"; }
41         string getDescription()         { return "align sequences"; }
42         
43         int execute(); 
44         void help() { m->mothurOut(getHelpString()); }  
45         
46 private:
47         struct linePair {
48                 unsigned long long start;
49                 unsigned long long end;
50                 linePair(unsigned long long i, unsigned long long j) : start(i), end(j) {}
51         };
52         vector<int> processIDS;   //processid
53         vector<linePair*> lines;
54         bool MPIWroteAccnos;
55         
56         AlignmentDB* templateDB;
57         
58         int driver(linePair*, string, string, string, string);
59         int createProcesses(string, string, string, string);
60         void appendAlignFiles(string, string); 
61         void appendReportFiles(string, string);
62         
63         #ifdef USE_MPI
64         int driverMPI(int, int, MPI_File&, MPI_File&, MPI_File&, MPI_File&, vector<unsigned long long>&);
65         #endif
66         
67         string candidateFileName, templateFileName, distanceFileName, search, align, outputDir;
68         float match, misMatch, gapOpen, gapExtend, threshold;
69         int processors, kmerSize;
70         vector<string> candidateFileNames;
71         vector<string> outputNames;
72         
73         bool abort, flip, calledHelp, save;
74
75 };
76
77 /**************************************************************************************************/
78 //custom data structure for threads to use.
79 // This is passed by void pointer so it can be any data type
80 // that can be passed using a single void pointer (LPVOID).
81 struct alignData {
82         string templateFileName;
83         string alignFName; 
84         string reportFName; 
85         string accnosFName;
86         string filename;
87         string align;
88         string search;
89         bool flip;
90         unsigned long long start;
91         unsigned long long end;
92         MothurOut* m;
93         //AlignmentDB* templateDB;
94         float match, misMatch, gapOpen, gapExtend, threshold;
95         int count, kmerSize, threadID;
96         
97         alignData(){}
98         alignData(string te, string a, string r, string ac, string f, string al, string se, int ks, MothurOut* mout, unsigned long long st, unsigned long long en, bool fl, float ma, float misMa, float gapO, float gapE, float thr, int tid) {
99                 templateFileName = te;
100                 alignFName = a;
101                 reportFName = r;
102                 accnosFName = ac;
103                 filename = f;
104                 flip = fl;
105                 m = mout;
106                 start = st;
107                 end = en;
108                 //templateDB = tdb;
109                 match = ma; 
110                 misMatch = misMa;
111                 gapOpen = gapO; 
112                 gapExtend = gapE; 
113                 threshold = thr;
114                 align = al;
115                 search = se;
116                 count = 0;
117                 kmerSize = ks;
118                 threadID = tid;
119         }
120 };
121
122 /**************************************************************************************************/
123 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
124 #else
125 static DWORD WINAPI MyAlignThreadFunction(LPVOID lpParam){ 
126         alignData* pDataArray;
127         pDataArray = (alignData*)lpParam;
128         
129         try {
130                 ofstream alignmentFile;
131                 pDataArray->m->openOutputFile(pDataArray->alignFName, alignmentFile);
132                 
133                 ofstream accnosFile;
134                 pDataArray->m->openOutputFile(pDataArray->accnosFName, accnosFile);
135                 
136                 NastReport report(pDataArray->reportFName);
137                 
138                 ifstream inFASTA;
139                 pDataArray->m->openInputFile(pDataArray->filename, inFASTA);
140                 
141                 //print header if you are process 0
142                 if ((pDataArray->start == 0) || (pDataArray->start == 1)) {
143                         inFASTA.seekg(0);
144                 }else { //this accounts for the difference in line endings. 
145                         inFASTA.seekg(pDataArray->start-1); pDataArray->m->gobble(inFASTA); 
146                 }
147                 
148                 pDataArray->count = pDataArray->end;
149                 
150                 AlignmentDB* templateDB = new AlignmentDB(pDataArray->templateFileName, pDataArray->search, pDataArray->kmerSize, pDataArray->gapOpen, pDataArray->gapExtend, pDataArray->match, pDataArray->misMatch, pDataArray->threadID);
151                 
152                 //moved this into driver to avoid deep copies in windows paralellized version
153                 Alignment* alignment;
154                 int longestBase = templateDB->getLongestBase();
155                 if(pDataArray->align == "gotoh")                        {       alignment = new GotohOverlap(pDataArray->gapOpen, pDataArray->gapExtend, pDataArray->match, pDataArray->misMatch, longestBase);                 }
156                 else if(pDataArray->align == "needleman")       {       alignment = new NeedlemanOverlap(pDataArray->gapOpen, pDataArray->match, pDataArray->misMatch, longestBase);                            }
157                 else if(pDataArray->align == "blast")           {       alignment = new BlastAlignment(pDataArray->gapOpen, pDataArray->gapExtend, pDataArray->match, pDataArray->misMatch);            }
158                 else if(pDataArray->align == "noalign")         {       alignment = new NoAlign();                                                                                                      }
159                 else {
160                         pDataArray->m->mothurOut(pDataArray->align + " is not a valid alignment option. I will run the command using needleman.");
161                         pDataArray->m->mothurOutEndLine();
162                         alignment = new NeedlemanOverlap(pDataArray->gapOpen, pDataArray->match, pDataArray->misMatch, longestBase);
163                 }
164                 
165                 int count = 0;
166                 for(int i = 0; i < pDataArray->end; i++){ //end is the number of sequences to process
167                         
168                         if (pDataArray->m->control_pressed) {  break; }
169                         
170                         Sequence* candidateSeq = new Sequence(inFASTA);  pDataArray->m->gobble(inFASTA);
171                         report.setCandidate(candidateSeq);
172                         
173                         int origNumBases = candidateSeq->getNumBases();
174                         string originalUnaligned = candidateSeq->getUnaligned();
175                         int numBasesNeeded = origNumBases * pDataArray->threshold;
176                         
177                         if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
178                                 if (candidateSeq->getUnaligned().length() > alignment->getnRows()) {
179                                         alignment->resize(candidateSeq->getUnaligned().length()+1);
180                                 }
181                                 
182                                 Sequence temp = templateDB->findClosestSequence(candidateSeq);
183                                 Sequence* templateSeq = &temp;
184                                 
185                                 float searchScore = templateDB->getSearchScore();
186                                 
187                                 Nast* nast = new Nast(alignment, candidateSeq, templateSeq);
188                                 
189                                 Sequence* copy;
190                                 
191                                 Nast* nast2;
192                                 bool needToDeleteCopy = false;  //this is needed in case you have you enter the ifs below
193                                 //since nast does not make a copy of hte sequence passed, and it is used by the reporter below
194                                 //you can't delete the copy sequence til after you report, but you may choose not to create it in the first place
195                                 //so this bool tells you if you need to delete it
196                                 
197                                 //if there is a possibility that this sequence should be reversed
198                                 if (candidateSeq->getNumBases() < numBasesNeeded) {
199                                         
200                                         string wasBetter =  "";
201                                         //if the user wants you to try the reverse
202                                         if (pDataArray->flip) {
203                                                 
204                                                 //get reverse compliment
205                                                 copy = new Sequence(candidateSeq->getName(), originalUnaligned);
206                                                 copy->reverseComplement();
207                                                 
208                                                 //rerun alignment
209                                                 Sequence temp2 = templateDB->findClosestSequence(copy);
210                                                 Sequence* templateSeq2 = &temp2;
211                                                 
212                                                 searchScore = templateDB->getSearchScore();
213                                                 
214                                                 nast2 = new Nast(alignment, copy, templateSeq2);
215                                                 
216                                                 //check if any better
217                                                 if (copy->getNumBases() > candidateSeq->getNumBases()) {
218                                                         candidateSeq->setAligned(copy->getAligned());  //use reverse compliments alignment since its better
219                                                         templateSeq = templateSeq2; 
220                                                         delete nast;
221                                                         nast = nast2;
222                                                         needToDeleteCopy = true;
223                                                         wasBetter = "\treverse complement produced a better alignment, so mothur used the reverse complement.";
224                                                 }else{  
225                                                         wasBetter = "\treverse complement did NOT produce a better alignment so it was not used, please check sequence.";
226                                                         delete nast2;
227                                                         delete copy;    
228                                                 }
229                                         }
230                                         
231                                         //create accnos file with names
232                                         accnosFile << candidateSeq->getName() << wasBetter << endl;
233                                 }
234                                 
235                                 report.setTemplate(templateSeq);
236                                 report.setSearchParameters(pDataArray->search, searchScore);
237                                 report.setAlignmentParameters(pDataArray->align, alignment);
238                                 report.setNastParameters(*nast);
239                                 
240                                 alignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl;
241                                 
242                                 report.print();
243                                 delete nast;
244                                 if (needToDeleteCopy) {   delete copy;   }
245                                 
246                                 count++;
247                         }
248                         delete candidateSeq;
249                         
250                         //report progress
251                         if((count) % 100 == 0){ pDataArray->m->mothurOut(toString(count)); pDataArray->m->mothurOutEndLine();           }
252                         
253                 }
254                 //report progress
255                 if((count) % 100 != 0){ pDataArray->m->mothurOut(toString(count)); pDataArray->m->mothurOutEndLine();           }
256                 
257                 delete alignment;
258                 delete templateDB;
259                 alignmentFile.close();
260                 inFASTA.close();
261                 accnosFile.close();
262         }
263         catch(exception& e) {
264                 pDataArray->m->errorOut(e, "AlignCommand", "MyAlignThreadFunction");
265                 exit(1);
266         }
267
268 #endif
269
270
271
272 #endif