]> git.donarmstrong.com Git - mothur.git/blob - shhhseqscommand.h
fixed bug in phylo.diversity rooting. added filename patterns and create filename...
[mothur.git] / shhhseqscommand.h
1 #ifndef SHHHSEQSCOMMAND_H
2 #define SHHHSEQSCOMMAND_H
3
4 /*
5  *  shhhseqscommand.h
6  *  Mothur
7  *
8  *  Created by westcott on 11/8/11.
9  *  Copyright 2011 Schloss Lab. All rights reserved.
10  *
11  */
12
13 #include "command.hpp"
14 #include "myseqdist.h"
15 #include "seqnoise.h"
16 #include "sequenceparser.h"
17 #include "deconvolutecommand.h"
18 #include "clustercommand.h"
19
20 //**********************************************************************************************************************
21
22 class ShhhSeqsCommand : public Command {
23         
24 public:
25         ShhhSeqsCommand(string);
26         ShhhSeqsCommand();
27         ~ShhhSeqsCommand() {}
28         
29         vector<string> setParameters();
30         string getCommandName()                 { return "shhh.seqs";   }
31         string getCommandCategory()             { return "Sequence Processing";         }
32         
33         string getHelpString(); 
34     string getOutputPattern(string);    
35         string getCitation() { return "Schloss PD, Gevers D, Westcott SL (2011).  Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies.  PLoS ONE.  6:e27310.\nQuince C, Lanzen A, Davenport RJ, Turnbaugh PJ (2011).  Removing noise from pyrosequenced amplicons.  BMC Bioinformatics  12:38.\nhttp://www.mothur.org/wiki/Shhh.seqs"; }
36         string getDescription()         { return "shhh.seqs"; }
37         
38         
39         int execute(); 
40         void help() { m->mothurOut(getHelpString()); }          
41         
42 private:
43         
44         struct linePair {
45                 int start;
46                 int end;
47                 linePair(int i, int j) : start(i), end(j) {}
48         };
49         
50         bool abort;
51         string outputDir, fastafile, namefile, groupfile;
52         int processors;
53         double sigma;
54         vector<string> outputNames;
55         
56         int readData(correctDist*, seqNoise&, vector<string>&, vector<string>&, vector<string>&, vector<int>&);
57         int loadData(correctDist*, seqNoise&, vector<string>&, vector<string>&, vector<string>&, vector<int>&, map<string, string>&, vector<Sequence>&);
58
59         int driver(seqNoise&, vector<string>&, vector<string>&, vector<string>&, vector<int>&, string, string, string, string);
60         vector<string> driverGroups(SequenceParser&, string, string, string, int, int, vector<string>);
61         vector<string> createProcessesGroups(SequenceParser&, string, string, string, vector<string>);
62         int deconvoluteResults(string, string);
63
64
65
66 };
67
68 /**************************************************************************************************/
69 //custom data structure for threads to use.
70 // This is passed by void pointer so it can be any data type
71 // that can be passed using a single void pointer (LPVOID).
72 struct shhhseqsData {
73         string fastafile; 
74         string namefile; 
75         string groupfile;
76         string newFName, newNName, newMName;
77         MothurOut* m;
78         int start;
79         int end;
80         int sigma, threadID;
81         vector<string> groups;
82         vector<string> mapfileNames;
83         
84         shhhseqsData(){}
85         shhhseqsData(string f, string n, string g, string nff,  string nnf, string nmf, vector<string> gr, MothurOut* mout, int st, int en, int s, int tid) {
86                 fastafile = f;
87                 namefile = n;
88                 groupfile = g;
89                 newFName = nff;
90                 newNName = nnf;
91                 newMName = nmf;
92                 m = mout;
93                 start = st;
94                 end = en;
95                 sigma = s;
96                 threadID = tid;
97                 groups = gr;
98         }
99 };
100
101 /**************************************************************************************************/
102 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
103 #else
104 static DWORD WINAPI MyShhhSeqsThreadFunction(LPVOID lpParam){ 
105         shhhseqsData* pDataArray;
106         pDataArray = (shhhseqsData*)lpParam;
107         
108         try {
109                 
110                 //parse fasta and name file by group
111                 SequenceParser parser(pDataArray->groupfile, pDataArray->fastafile, pDataArray->namefile);
112                                                 
113                 //precluster each group
114                 for (int k = pDataArray->start; k < pDataArray->end; k++) {
115                         
116                         int start = time(NULL);
117                         
118                         if (pDataArray->m->control_pressed) {  return 0; }
119                         
120                         pDataArray->m->mothurOutEndLine(); pDataArray->m->mothurOut("Processing group " + pDataArray->groups[k] + ":"); pDataArray->m->mothurOutEndLine();
121
122                         map<string, string> thisNameMap;
123                         thisNameMap = parser.getNameMap(pDataArray->groups[k]); 
124                         vector<Sequence> thisSeqs = parser.getSeqs(pDataArray->groups[k]);
125                         
126                         if (pDataArray->m->control_pressed) {  return 0; }
127                         
128                         vector<string> sequences;
129                         vector<string> uniqueNames;
130                         vector<string> redundantNames;
131                         vector<int> seqFreq;
132                         
133                         seqNoise noise;
134                         correctDist* correct = new correctDist(1); //we use one processor since we already split up the work load.
135                         
136                         //load this groups info in order
137                         //loadData(correct, noise, sequences, uniqueNames, redundantNames, seqFreq, thisNameMap, thisSeqs);
138                         //////////////////////////////////////////////////////////////////////////////////////////////////
139                         bool error = false;
140                         map<string, string>::iterator it;
141                         
142                         for (int i = 0; i < thisSeqs.size(); i++) {
143                                 
144                                 if (pDataArray->m->control_pressed) { return 0; }
145                                 
146                                 if (thisSeqs[i].getName() != "") {
147                                         correct->addSeq(thisSeqs[i].getName(), thisSeqs[i].getAligned());
148                                         
149                                         it = thisNameMap.find(thisSeqs[i].getName());
150                                         if (it != thisNameMap.end()) {
151                                                 noise.addSeq(thisSeqs[i].getAligned(), sequences);
152                                                 noise.addRedundantName(it->first, it->second, uniqueNames, redundantNames, seqFreq);
153                                         }else {
154                                                 pDataArray->m->mothurOut("[ERROR]: " + thisSeqs[i].getName() + " is in your fasta file and not in your namefile, please correct.");
155                                                 error = true;
156                                         }
157                                 }
158                         }
159                         
160                         if (error) { return 0; }
161                         //////////////////////////////////////////////////////////////////////////////////////////////////
162                         
163                         if (pDataArray->m->control_pressed) { return 0; }
164                         
165                         //calc distances for cluster
166                         string distFileName = pDataArray->m->getRootName(pDataArray->m->getSimpleName(pDataArray->fastafile)) + pDataArray->groups[k] + ".shhh.dist";
167                         correct->execute(distFileName);
168                         delete correct;
169                         
170                         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(distFileName); return 0; }
171                         
172                         //driver(noise, sequences, uniqueNames, redundantNames, seqFreq, distFileName, newFFile+groups[i], newNFile+groups[i], newMFile+groups[i]+".map"); 
173                         ///////////////////////////////////////////////////////////////////////////////////////////////////
174                         double cutOff = 0.08;
175                         int minIter = 10;
176                         int maxIter = 1000;
177                         double minDelta = 1e-6;
178                         int numIters = 0;
179                         double maxDelta = 1e6;
180                         int numSeqs = sequences.size();
181                         
182                         //run cluster command
183                         string inputString = "phylip=" + distFileName + ", method=furthest, cutoff=0.08";
184                         pDataArray->m->mothurOut("/******************************************/"); pDataArray->m->mothurOutEndLine(); 
185                         pDataArray->m->mothurOut("Running command: cluster(" + inputString + ")"); pDataArray->m->mothurOutEndLine(); 
186                         
187                         Command* clusterCommand = new ClusterCommand(inputString);
188                         clusterCommand->execute();
189                         
190                         map<string, vector<string> > filenames = clusterCommand->getOutputFiles();
191                         string listFileName = filenames["list"][0];
192                         string rabundFileName = filenames["rabund"][0]; pDataArray->m->mothurRemove(rabundFileName);
193                         string sabundFileName = filenames["sabund"][0]; pDataArray->m->mothurRemove(sabundFileName);
194                         
195                         delete clusterCommand;
196                         pDataArray->m->mothurOut("/******************************************/"); pDataArray->m->mothurOutEndLine(); 
197                         
198                         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(distFileName); pDataArray->m->mothurRemove(listFileName); return 0; }
199                         
200                         vector<double> distances(numSeqs * numSeqs);
201                         noise.getDistanceData(distFileName, distances);
202                         pDataArray->m->mothurRemove(distFileName); 
203                         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(listFileName); return 0; }
204                         
205                         vector<int> otuData(numSeqs);
206                         vector<int> otuFreq;
207                         vector<vector<int> > otuBySeqLookUp;
208                         noise.getListData(listFileName, cutOff, otuData, otuFreq, otuBySeqLookUp);
209                         pDataArray->m->mothurRemove(listFileName);
210                         if (pDataArray->m->control_pressed) { return 0; }
211                         
212                         int numOTUs = otuFreq.size();
213                         
214                         vector<double> weights(numOTUs, 0);
215                         vector<int> change(numOTUs, 1);
216                         vector<int> centroids(numOTUs, -1);
217                         vector<int> cumCount(numOTUs, 0);
218                         
219                         vector<double> tau(numSeqs, 1);
220                         vector<int> anP(numSeqs, 0);
221                         vector<int> anI(numSeqs, 0);
222                         vector<int> anN(numSeqs, 0);
223                         vector<vector<int> > aanI = otuBySeqLookUp;
224                         
225                         while(numIters < minIter || ((maxDelta > minDelta) && (numIters < maxIter))){
226                                 
227                                 if (pDataArray->m->control_pressed) { return 0; }
228                                 
229                                 noise.updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount); if (pDataArray->m->control_pressed) { return 0; }
230                                 maxDelta = noise.calcNewWeights(weights, seqFreq, anI, cumCount, anP, otuFreq, tau);  if (pDataArray->m->control_pressed) { return 0; }
231                                 
232                                 noise.calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau); if (pDataArray->m->control_pressed) { return 0; }
233                                 noise.checkCentroids(weights, centroids); if (pDataArray->m->control_pressed) { return 0; }
234                                 
235                                 otuFreq.assign(numOTUs, 0);
236                                 
237                                 int total = 0;
238                                 
239                                 for(int i=0;i<numSeqs;i++){
240                                         if (pDataArray->m->control_pressed) { return 0; }
241                                         
242                                         double offset = 1e6;
243                                         double norm = 0.0000;
244                                         double minWeight = 0.1;
245                                         vector<double> currentTau(numOTUs);
246                                         
247                                         for(int j=0;j<numOTUs;j++){
248                                                 if (pDataArray->m->control_pressed) { return 0; }
249                                                 if(weights[j] > minWeight && distances[i * numSeqs+centroids[j]] < offset){
250                                                         offset = distances[i * numSeqs+centroids[j]];
251                                                 }
252                                         }
253                                         
254                                         for(int j=0;j<numOTUs;j++){
255                                                 if (pDataArray->m->control_pressed) { return 0; }
256                                                 if(weights[j] > minWeight){
257                                                         currentTau[j] = exp(pDataArray->sigma * (-distances[(i * numSeqs + centroids[j])] + offset)) * weights[j];
258                                                         norm += currentTau[j];
259                                                 }
260                                                 else{
261                                                         currentTau[j] = 0.0000;
262                                                 }
263                                         }                       
264                                         
265                                         for(int j=0;j<numOTUs;j++){
266                                                 if (pDataArray->m->control_pressed) { return 0; }
267                                                 currentTau[j] /= norm;
268                                         }
269                                         
270                                         for(int j=0;j<numOTUs;j++){
271                                                 if (pDataArray->m->control_pressed) { return 0; }
272                                                 
273                                                 if(currentTau[j] > 1.0e-4){
274                                                         int oldTotal = total;
275                                                         total++;
276                                                         
277                                                         tau.resize(oldTotal+1);
278                                                         tau[oldTotal] = currentTau[j];
279                                                         otuBySeqLookUp[j][otuFreq[j]] = oldTotal;
280                                                         aanI[j][otuFreq[j]] = i;
281                                                         otuFreq[j]++;
282                                                         
283                                                 }
284                                         }
285                                         
286                                         anP.resize(total);
287                                         anI.resize(total);
288                                 }
289                                 
290                                 numIters++;
291                         }
292                         
293                         noise.updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);  if (pDataArray->m->control_pressed) { return 0; }
294                         
295                         vector<double> percentage(numSeqs);
296                         noise.setUpOTUData(otuData, percentage, cumCount, tau, otuFreq, anP, anI);  if (pDataArray->m->control_pressed) { return 0; }
297                         noise.finishOTUData(otuData, otuFreq, anP, anI, cumCount, otuBySeqLookUp, aanI, tau);  if (pDataArray->m->control_pressed) { return 0; }
298                         
299                         change.assign(numOTUs, 1);
300                         noise.calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau); if (pDataArray->m->control_pressed) { return 0; }
301                         
302                         
303                         vector<int> finalTau(numOTUs, 0);
304                         for(int i=0;i<numSeqs;i++){
305                                 if (pDataArray->m->control_pressed) { return 0; }
306                                 finalTau[otuData[i]] += int(seqFreq[i]);
307                         }
308                         
309                         noise.writeOutput(pDataArray->newFName+pDataArray->groups[k], pDataArray->newNName+pDataArray->groups[k], pDataArray->newMName+pDataArray->groups[k]+".map", finalTau, centroids, otuData, sequences, uniqueNames, redundantNames, seqFreq, distances);
310                         
311                         ///////////////////////////////////////////////////////////////////////////////////////////////////
312                         
313                         if (pDataArray->m->control_pressed) { return 0; }
314                         
315                         pDataArray->m->appendFiles(pDataArray->newFName+pDataArray->groups[k], pDataArray->newFName); pDataArray->m->mothurRemove(pDataArray->newFName+pDataArray->groups[k]);
316                         pDataArray->m->appendFiles(pDataArray->newNName+pDataArray->groups[k], pDataArray->newNName); pDataArray->m->mothurRemove(pDataArray->newNName+pDataArray->groups[k]);
317                         pDataArray->mapfileNames.push_back(pDataArray->newMName+pDataArray->groups[k]+".map");
318                         
319                         pDataArray->m->mothurOut("It took " + toString(time(NULL) - start) + " secs to process group " + pDataArray->groups[k] + "."); pDataArray->m->mothurOutEndLine(); 
320                 }
321                 
322                 return 0;
323                 
324                 
325         }
326         catch(exception& e) {
327                 pDataArray->m->errorOut(e, "ShhhSeqsCommand", "MyShhhSeqsThreadFunction");
328                 exit(1);
329         }
330
331 #endif
332 /**************************************************************************************************/
333
334 #endif