5 * Created by Pat Schloss on 12/27/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "shhhercommand.h"
12 #include "readcolumn.h"
13 #include "readmatrix.hpp"
14 #include "rabundvector.hpp"
15 #include "sabundvector.hpp"
16 #include "listvector.hpp"
17 #include "cluster.hpp"
18 #include "sparsematrix.hpp"
21 //**********************************************************************************************************************
26 #define MIN_WEIGHT 0.1
27 #define MIN_TAU 0.0001
29 //**********************************************************************************************************************
30 vector<string> ShhherCommand::setParameters(){
32 CommandParameter pflow("flow", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pflow);
33 CommandParameter pfile("file", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pfile);
34 CommandParameter plookup("lookup", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(plookup);
35 CommandParameter pcutoff("cutoff", "Number", "", "0.01", "", "", "",false,false); parameters.push_back(pcutoff);
36 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
37 CommandParameter pmaxiter("maxiter", "Number", "", "1000", "", "", "",false,false); parameters.push_back(pmaxiter);
38 CommandParameter psigma("sigma", "Number", "", "60", "", "", "",false,false); parameters.push_back(psigma);
39 CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "",false,false); parameters.push_back(pmindelta);
40 CommandParameter porder("order", "String", "", "", "", "", "",false,false); parameters.push_back(porder);
41 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
42 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
44 vector<string> myArray;
45 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
49 m->errorOut(e, "ShhherCommand", "setParameters");
53 //**********************************************************************************************************************
54 string ShhherCommand::getHelpString(){
56 string helpString = "";
57 helpString += "The shhh.seqs command reads a file containing flowgrams and creates a file of corrected sequences.\n";
61 m->errorOut(e, "ShhherCommand", "getHelpString");
65 //**********************************************************************************************************************
67 ShhherCommand::ShhherCommand(){
69 abort = true; calledHelp = true;
72 //initialize outputTypes
73 // vector<string> tempOutNames;
74 // outputTypes["pn.dist"] = tempOutNames;
78 m->errorOut(e, "ShhherCommand", "ShhherCommand");
83 //**********************************************************************************************************************
85 ShhherCommand::ShhherCommand(string option) {
89 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
90 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
96 abort = false; calledHelp = false;
99 //allow user to run help
100 if(option == "help") { help(); abort = true; calledHelp = true; }
101 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
104 vector<string> myArray = setParameters();
106 OptionParser parser(option);
107 map<string,string> parameters = parser.getParameters();
109 ValidParameters validParameter;
110 map<string,string>::iterator it;
112 //check to make sure all parameters are valid for command
113 for (it = parameters.begin(); it != parameters.end(); it++) {
114 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
117 //initialize outputTypes
118 vector<string> tempOutNames;
119 // outputTypes["pn.dist"] = tempOutNames;
120 // outputTypes["fasta"] = tempOutNames;
122 //if the user changes the input directory command factory will send this info to us in the output parameter
123 string inputDir = validParameter.validFile(parameters, "inputdir", false);
124 if (inputDir == "not found"){ inputDir = ""; }
127 it = parameters.find("flow");
128 //user has given a template file
129 if(it != parameters.end()){
130 path = m->hasPath(it->second);
131 //if the user has not given a path then, add inputdir. else leave path alone.
132 if (path == "") { parameters["flow"] = inputDir + it->second; }
135 it = parameters.find("lookup");
136 //user has given a template file
137 if(it != parameters.end()){
138 path = m->hasPath(it->second);
139 //if the user has not given a path then, add inputdir. else leave path alone.
140 if (path == "") { parameters["lookup"] = inputDir + it->second; }
143 it = parameters.find("file");
144 //user has given a template file
145 if(it != parameters.end()){
146 path = m->hasPath(it->second);
147 //if the user has not given a path then, add inputdir. else leave path alone.
148 if (path == "") { parameters["file"] = inputDir + it->second; }
153 //check for required parameters
154 flowFileName = validParameter.validFile(parameters, "flow", true);
155 flowFilesFileName = validParameter.validFile(parameters, "file", true);
156 if (flowFileName == "not found" && flowFilesFileName == "not found") {
157 m->mothurOut("values for either flow or file must be provided for the shhh.seqs command.");
158 m->mothurOutEndLine();
161 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
163 if(flowFileName != "not found"){
164 compositeFASTAFileName = "";
165 compositeNamesFileName = "";
170 //flow.files = 9 character offset
171 compositeFASTAFileName = flowFilesFileName.substr(0, flowFilesFileName.length()-10) + "shhh.fasta";
172 m->openOutputFile(compositeFASTAFileName, temp);
175 compositeNamesFileName = flowFilesFileName.substr(0, flowFilesFileName.length()-10) + "shhh.names";
176 m->openOutputFile(compositeNamesFileName, temp);
180 //if the user changes the output directory command factory will send this info to us in the output parameter
181 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
183 outputDir += m->hasPath(flowFileName); //if user entered a file with a path then preserve it
187 //check for optional parameter and set defaults
188 // ...at some point should added some additional type checking...
190 temp = validParameter.validFile(parameters, "lookup", true);
191 if (temp == "not found") {
192 lookupFileName = "LookUp_Titanium.pat";
196 ableToOpen = m->openInputFile(lookupFileName, in, "noerror");
199 //if you can't open it, try input location
200 if (ableToOpen == 1) {
201 if (inputDir != "") { //default path is set
202 string tryPath = inputDir + lookupFileName;
203 m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
205 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
207 lookupFileName = tryPath;
211 //if you can't open it, try default location
212 if (ableToOpen == 1) {
213 if (m->getDefaultPath() != "") { //default path is set
214 string tryPath = m->getDefaultPath() + m->getSimpleName(lookupFileName);
215 m->mothurOut("Unable to open " + lookupFileName + ". Trying default " + tryPath); m->mothurOutEndLine();
217 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
219 lookupFileName = tryPath;
223 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
224 if (ableToOpen == 1) {
225 string exepath = m->argv;
226 string tempPath = exepath;
227 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
228 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
230 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
231 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
233 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
235 lookupFileName = tryPath;
238 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
240 else if(temp == "not open") {
242 lookupFileName = validParameter.validFile(parameters, "lookup", false);
244 //if you can't open it its not inputDir, try mothur excutable location
245 string exepath = m->argv;
246 string tempPath = exepath;
247 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
248 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
250 string tryPath = m->getFullPathName(exepath) + lookupFileName;
251 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
253 int ableToOpen = m->openInputFile(tryPath, in2, "noerror");
255 lookupFileName = tryPath;
257 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
258 }else { lookupFileName = temp; }
260 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
261 m->setProcessors(temp);
262 convert(temp, processors);
264 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
265 convert(temp, cutoff);
267 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
268 convert(temp, minDelta);
270 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
271 convert(temp, maxIters);
273 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
274 convert(temp, sigma);
276 flowOrder = validParameter.validFile(parameters, "order", false);
277 if (flowOrder == "not found"){ flowOrder = "TACG"; }
278 else if(flowOrder.length() != 4){
279 m->mothurOut("The value of the order option must be four bases long\n");
289 catch(exception& e) {
290 m->errorOut(e, "ShhherCommand", "ShhherCommand");
294 //**********************************************************************************************************************
296 int ShhherCommand::execute(){
298 if (abort == true) { if (calledHelp) { return 0; } return 2; }
305 for(int i=1;i<ncpus;i++){
306 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
308 if(abort == 1){ return 0; }
312 m->mothurOut("\nGetting preliminary data...\n");
316 vector<string> flowFileVector;
317 if(flowFilesFileName != "not found"){
320 ifstream flowFilesFile;
321 m->openInputFile(flowFilesFileName, flowFilesFile);
322 while(flowFilesFile){
323 flowFilesFile >> fName;
324 flowFileVector.push_back(fName);
325 m->gobble(flowFilesFile);
329 flowFileVector.push_back(flowFileName);
331 int numFiles = flowFileVector.size();
333 for(int i=1;i<ncpus;i++){
334 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
337 for(int i=0;i<numFiles;i++){
338 double begClock = clock();
339 unsigned long int begTime = time(NULL);
341 flowFileName = flowFileVector[i];
343 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
344 m->mothurOut("Reading flowgrams...\n");
347 m->mothurOut("Identifying unique flowgrams...\n");
350 m->mothurOut("Calculating distances between flowgrams...\n");
352 strcpy(fileName, flowFileName.c_str());
354 for(int i=1;i<ncpus;i++){
355 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
357 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
358 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
359 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
360 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
361 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
362 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
363 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
364 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
365 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
366 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
369 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
372 for(int i=1;i<ncpus;i++){
373 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
375 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
376 remove((distFileName + ".temp." + toString(i)).c_str());
379 string namesFileName = createNamesFile();
381 m->mothurOut("\nClustering flowgrams...\n");
382 string listFileName = cluster(distFileName, namesFileName);
384 getOTUData(listFileName);
386 remove(distFileName.c_str());
387 remove(namesFileName.c_str());
388 remove(listFileName.c_str());
392 for(int i=1;i<ncpus;i++){
393 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
394 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
395 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
396 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
403 int numOTUsOnCPU = numOTUs / ncpus;
404 int numSeqsOnCPU = numSeqs / ncpus;
405 m->mothurOut("\nDenoising flowgrams...\n");
406 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
408 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
410 double cycClock = clock();
411 unsigned long int cycTime = time(NULL);
414 int total = singleTau.size();
415 for(int i=1;i<ncpus;i++){
416 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
417 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
418 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
420 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
421 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
422 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
423 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
424 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
427 calcCentroidsDriver(0, numOTUsOnCPU);
429 for(int i=1;i<ncpus;i++){
430 int otuStart = i * numOTUs / ncpus;
431 int otuStop = (i + 1) * numOTUs / ncpus;
433 vector<int> tempCentroids(numOTUs, 0);
434 vector<short> tempChange(numOTUs, 0);
436 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
437 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
439 for(int j=otuStart;j<otuStop;j++){
440 centroids[j] = tempCentroids[j];
441 change[j] = tempChange[j];
445 maxDelta = getNewWeights();
446 double nLL = getLikelihood();
449 for(int i=1;i<ncpus;i++){
450 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
451 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
452 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
455 calcNewDistancesParent(0, numSeqsOnCPU);
457 total = singleTau.size();
459 for(int i=1;i<ncpus;i++){
461 int seqStart = i * numSeqs / ncpus;
462 int seqStop = (i + 1) * numSeqs / ncpus;
464 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
466 vector<int> childSeqIndex(childTotal, 0);
467 vector<double> childSingleTau(childTotal, 0);
468 vector<double> childDist(numSeqs * numOTUs, 0);
469 vector<int> otuIndex(childTotal, 0);
471 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
472 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
473 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
474 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
476 int oldTotal = total;
478 singleTau.resize(total, 0);
479 seqIndex.resize(total, 0);
480 seqNumber.resize(total, 0);
484 for(int j=oldTotal;j<total;j++){
485 int otuI = otuIndex[childIndex];
486 int seqI = childSeqIndex[childIndex];
488 singleTau[j] = childSingleTau[childIndex];
490 aaP[otuI][nSeqsPerOTU[otuI]] = j;
491 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
496 int index = seqStart * numOTUs;
497 for(int j=seqStart;j<seqStop;j++){
498 for(int k=0;k<numOTUs;k++){
499 dist[index] = childDist[index];
507 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
509 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
511 for(int i=1;i<ncpus;i++){
512 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
517 for(int i=1;i<ncpus;i++){
518 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
524 m->mothurOut("\nFinalizing...\n");
528 vector<int> otuCounts(numOTUs, 0);
529 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
530 calcCentroidsDriver(0, numOTUs);
532 writeQualities(otuCounts);
533 writeSequences(otuCounts);
534 writeNames(otuCounts);
535 writeClusters(otuCounts);
539 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
545 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
546 if(abort){ return 0; }
549 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
551 for(int i=0;i<numFiles;i++){
552 //Now into the pyrodist part
556 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
557 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
558 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
559 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
561 flowDataIntI.resize(numSeqs * numFlowCells);
562 flowDataPrI.resize(numSeqs * numFlowCells);
563 mapUniqueToSeq.resize(numSeqs);
564 mapSeqToUnique.resize(numSeqs);
565 lengths.resize(numSeqs);
566 jointLookUp.resize(NUMBINS * NUMBINS);
568 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
569 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
570 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
571 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
572 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
573 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
574 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
576 flowFileName = string(fileName);
577 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
578 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
580 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
583 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
585 //Now into the pyronoise part
586 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
588 singleLookUp.resize(HOMOPS * NUMBINS);
589 uniqueFlowgrams.resize(numUniques * numFlowCells);
590 weight.resize(numOTUs);
591 centroids.resize(numOTUs);
592 change.resize(numOTUs);
593 dist.assign(numOTUs * numSeqs, 0);
594 nSeqsPerOTU.resize(numOTUs);
595 cumNumSeqs.resize(numOTUs);
597 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
598 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
599 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
601 int startOTU = pid * numOTUs / ncpus;
602 int endOTU = (pid + 1) * numOTUs / ncpus;
604 int startSeq = pid * numSeqs / ncpus;
605 int endSeq = (pid + 1) * numSeqs /ncpus;
611 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
612 singleTau.assign(total, 0.0000);
613 seqNumber.assign(total, 0);
614 seqIndex.assign(total, 0);
616 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
617 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
618 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
619 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
620 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
621 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
622 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
624 calcCentroidsDriver(startOTU, endOTU);
626 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
627 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
629 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
630 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
631 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
633 vector<int> otuIndex(total, 0);
634 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
635 total = otuIndex.size();
637 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
638 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
639 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
640 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
641 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
643 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
647 MPI_Barrier(MPI_COMM_WORLD);
650 if(compositeFASTAFileName != ""){
651 outputNames.push_back(compositeFASTAFileName);
652 outputNames.push_back(compositeNamesFileName);
655 m->mothurOutEndLine();
656 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
657 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
658 m->mothurOutEndLine();
663 catch(exception& e) {
664 m->errorOut(e, "ShhherCommand", "execute");
669 /**************************************************************************************************/
671 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
673 ostringstream outStream;
674 outStream.setf(ios::fixed, ios::floatfield);
675 outStream.setf(ios::dec, ios::basefield);
676 outStream.setf(ios::showpoint);
677 outStream.precision(6);
679 int begTime = time(NULL);
680 double begClock = clock();
682 for(int i=startSeq;i<stopSeq;i++){
683 for(int j=0;j<i;j++){
684 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
686 if(flowDistance < 1e-6){
687 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
689 else if(flowDistance <= cutoff){
690 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
694 m->mothurOut(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
698 m->mothurOut(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
700 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
701 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
703 ofstream distFile(fDistFileName.c_str());
704 distFile << outStream.str();
707 return fDistFileName;
709 catch(exception& e) {
710 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
716 //**********************************************************************************************************************
718 int ShhherCommand::execute(){
720 if (abort == true) { return 0; }
725 vector<string> flowFileVector;
726 if(flowFilesFileName != "not found"){
729 ifstream flowFilesFile;
730 m->openInputFile(flowFilesFileName, flowFilesFile);
731 while(flowFilesFile){
732 flowFilesFile >> fName;
733 flowFileVector.push_back(fName);
734 m->gobble(flowFilesFile);
738 flowFileVector.push_back(flowFileName);
740 int numFiles = flowFileVector.size();
743 for(int i=0;i<numFiles;i++){
744 flowFileName = flowFileVector[i];
746 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
747 m->mothurOut("Reading flowgrams...\n");
750 m->mothurOut("Identifying unique flowgrams...\n");
754 m->mothurOut("Calculating distances between flowgrams...\n");
755 string distFileName = createDistFile(processors);
756 string namesFileName = createNamesFile();
758 m->mothurOut("\nClustering flowgrams...\n");
759 string listFileName = cluster(distFileName, namesFileName);
760 getOTUData(listFileName);
762 remove(distFileName.c_str());
763 remove(namesFileName.c_str());
764 remove(listFileName.c_str());
771 double begClock = clock();
772 unsigned long int begTime = time(NULL);
774 m->mothurOut("\nDenoising flowgrams...\n");
775 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
777 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
779 double cycClock = clock();
780 unsigned long int cycTime = time(NULL);
785 maxDelta = getNewWeights();
786 double nLL = getLikelihood();
793 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
797 m->mothurOut("\nFinalizing...\n");
801 vector<int> otuCounts(numOTUs, 0);
802 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
804 calcCentroidsDriver(0, numOTUs);
805 writeQualities(otuCounts);
806 writeSequences(otuCounts);
807 writeNames(otuCounts);
808 writeClusters(otuCounts);
811 m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
814 if(compositeFASTAFileName != ""){
815 outputNames.push_back(compositeFASTAFileName);
816 outputNames.push_back(compositeNamesFileName);
819 m->mothurOutEndLine();
820 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
821 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
822 m->mothurOutEndLine();
826 catch(exception& e) {
827 m->errorOut(e, "ShhherCommand", "execute");
832 /**************************************************************************************************/
834 void ShhherCommand::getFlowData(){
837 m->openInputFile(flowFileName, flowFile);
840 seqNameVector.clear();
842 flowDataIntI.clear();
846 int currentNumFlowCells;
850 flowFile >> numFlowCells;
851 int index = 0;//pcluster
852 while(!flowFile.eof()){
853 flowFile >> seqName >> currentNumFlowCells;
854 lengths.push_back(currentNumFlowCells);
856 seqNameVector.push_back(seqName);
857 nameMap[seqName] = index++;//pcluster
859 for(int i=0;i<numFlowCells;i++){
860 flowFile >> intensity;
861 if(intensity > 9.99) { intensity = 9.99; }
862 int intI = int(100 * intensity + 0.0001);
863 flowDataIntI.push_back(intI);
869 numSeqs = seqNameVector.size();
871 for(int i=0;i<numSeqs;i++){
872 int iNumFlowCells = i * numFlowCells;
873 for(int j=lengths[i];j<numFlowCells;j++){
874 flowDataIntI[iNumFlowCells + j] = 0;
879 catch(exception& e) {
880 m->errorOut(e, "ShhherCommand", "getFlowData");
885 /**************************************************************************************************/
887 void ShhherCommand::getSingleLookUp(){
889 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
890 singleLookUp.assign(HOMOPS * NUMBINS, 0);
894 m->openInputFile(lookupFileName, lookUpFile);
896 for(int i=0;i<HOMOPS;i++){
898 lookUpFile >> logFracFreq;
900 for(int j=0;j<NUMBINS;j++) {
901 lookUpFile >> singleLookUp[index];
907 catch(exception& e) {
908 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
913 /**************************************************************************************************/
915 void ShhherCommand::getJointLookUp(){
918 // the most likely joint probability (-log) that two intenities have the same polymer length
919 jointLookUp.resize(NUMBINS * NUMBINS, 0);
921 for(int i=0;i<NUMBINS;i++){
922 for(int j=0;j<NUMBINS;j++){
924 double minSum = 100000000;
926 for(int k=0;k<HOMOPS;k++){
927 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
929 if(sum < minSum) { minSum = sum; }
931 jointLookUp[i * NUMBINS + j] = minSum;
935 catch(exception& e) {
936 m->errorOut(e, "ShhherCommand", "getJointLookUp");
941 /**************************************************************************************************/
943 double ShhherCommand::getProbIntensity(int intIntensity){
946 double minNegLogProb = 100000000;
949 for(int i=0;i<HOMOPS;i++){//loop signal strength
950 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
951 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
954 return minNegLogProb;
956 catch(exception& e) {
957 m->errorOut(e, "ShhherCommand", "getProbIntensity");
962 /**************************************************************************************************/
964 void ShhherCommand::getUniques(){
969 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
970 uniqueCount.assign(numSeqs, 0); // anWeights
971 uniqueLengths.assign(numSeqs, 0);
972 mapSeqToUnique.assign(numSeqs, -1);
973 mapUniqueToSeq.assign(numSeqs, -1);
975 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
977 for(int i=0;i<numSeqs;i++){
980 vector<short> current(numFlowCells);
981 for(int j=0;j<numFlowCells;j++){
982 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
985 for(int j=0;j<numUniques;j++){
986 int offset = j * numFlowCells;
990 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
991 else { shorterLength = uniqueLengths[j]; }
993 for(int k=0;k<shorterLength;k++){
994 if(current[k] != uniqueFlowgrams[offset + k]){
1001 mapSeqToUnique[i] = j;
1004 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
1010 if(index == numUniques){
1011 uniqueLengths[numUniques] = lengths[i];
1012 uniqueCount[numUniques] = 1;
1013 mapSeqToUnique[i] = numUniques;//anMap
1014 mapUniqueToSeq[numUniques] = i;//anF
1016 for(int k=0;k<numFlowCells;k++){
1017 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
1018 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
1024 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
1025 uniqueLengths.resize(numUniques);
1027 flowDataPrI.resize(numSeqs * numFlowCells, 0);
1028 for(int i=0;i<flowDataPrI.size();i++) { flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
1030 catch(exception& e) {
1031 m->errorOut(e, "ShhherCommand", "getUniques");
1036 /**************************************************************************************************/
1038 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
1040 int minLength = lengths[mapSeqToUnique[seqA]];
1041 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
1043 int ANumFlowCells = seqA * numFlowCells;
1044 int BNumFlowCells = seqB * numFlowCells;
1048 for(int i=0;i<minLength;i++){
1049 int flowAIntI = flowDataIntI[ANumFlowCells + i];
1050 float flowAPrI = flowDataPrI[ANumFlowCells + i];
1052 int flowBIntI = flowDataIntI[BNumFlowCells + i];
1053 float flowBPrI = flowDataPrI[BNumFlowCells + i];
1054 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
1057 dist /= (float) minLength;
1060 catch(exception& e) {
1061 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
1066 /**************************************************************************************************/
1068 void ShhherCommand::flowDistParentFork(string distFileName, int startSeq, int stopSeq){
1071 ostringstream outStream;
1072 outStream.setf(ios::fixed, ios::floatfield);
1073 outStream.setf(ios::dec, ios::basefield);
1074 outStream.setf(ios::showpoint);
1075 outStream.precision(6);
1077 int begTime = time(NULL);
1078 double begClock = clock();
1080 for(int i=startSeq;i<stopSeq;i++){
1081 for(int j=0;j<i;j++){
1082 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
1084 if(flowDistance < 1e-6){
1085 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
1087 else if(flowDistance <= cutoff){
1088 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
1092 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
1093 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1094 m->mothurOutEndLine();
1097 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
1098 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1099 m->mothurOutEndLine();
1101 ofstream distFile(distFileName.c_str());
1102 distFile << outStream.str();
1105 catch(exception& e) {
1106 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
1111 /**************************************************************************************************/
1113 string ShhherCommand::createDistFile(int processors){
1115 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
1117 unsigned long int begTime = time(NULL);
1118 double begClock = clock();
1123 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1124 if(processors == 1) { flowDistParentFork(fDistFileName, 0, numUniques); }
1125 else{ //you have multiple processors
1127 if (numSeqs < processors){ processors = 1; }
1129 vector<int> start(processors, 0);
1130 vector<int> end(processors, 0);
1132 for (int i = 0; i < processors; i++) {
1133 start[i] = int(sqrt(float(i)/float(processors)) * numUniques);
1134 end[i] = int(sqrt(float(i+1)/float(processors)) * numUniques);
1138 vector<int> processIDs;
1140 //loop through and create all the processes you want
1141 while (process != processors) {
1145 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1147 }else if (pid == 0){
1148 flowDistParentFork(fDistFileName + toString(getpid()) + ".temp", start[process], end[process]);
1151 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1153 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1158 //parent does its part
1159 flowDistParentFork(fDistFileName, start[0], end[0]);
1161 //force parent to wait until all the processes are done
1162 for (int i=0;i<processIDs.size();i++) {
1163 int temp = processIDs[i];
1167 //append and remove temp files
1168 for (int i=0;i<processIDs.size();i++) {
1169 m->appendFiles((fDistFileName + toString(processIDs[i]) + ".temp"), fDistFileName);
1170 remove((fDistFileName + toString(processIDs[i]) + ".temp").c_str());
1176 flowDistParentFork(fDistFileName, 0, numUniques);
1179 m->mothurOutEndLine();
1181 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
1184 return fDistFileName;
1186 catch(exception& e) {
1187 m->errorOut(e, "ShhherCommand", "createDistFile");
1193 /**************************************************************************************************/
1195 string ShhherCommand::createNamesFile(){
1198 vector<string> duplicateNames(numUniques, "");
1199 for(int i=0;i<numSeqs;i++){
1200 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
1203 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
1206 m->openOutputFile(nameFileName, nameFile);
1208 for(int i=0;i<numUniques;i++){
1209 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1210 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1214 return nameFileName;
1216 catch(exception& e) {
1217 m->errorOut(e, "ShhherCommand", "createNamesFile");
1222 //**********************************************************************************************************************
1224 string ShhherCommand::cluster(string distFileName, string namesFileName){
1227 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1228 read->setCutoff(cutoff);
1230 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1231 clusterNameMap->readMap();
1232 read->read(clusterNameMap);
1234 ListVector* list = read->getListVector();
1235 SparseMatrix* matrix = read->getMatrix();
1238 delete clusterNameMap;
1240 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1242 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1243 string tag = cluster->getTag();
1245 double clusterCutoff = cutoff;
1246 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1247 cluster->update(clusterCutoff);
1250 list->setLabel(toString(cutoff));
1252 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
1254 m->openOutputFile(listFileName, listFile);
1255 list->print(listFile);
1258 delete matrix; delete cluster; delete rabund; delete list;
1260 return listFileName;
1262 catch(exception& e) {
1263 m->errorOut(e, "ShhherCommand", "cluster");
1268 /**************************************************************************************************/
1270 void ShhherCommand::getOTUData(string listFileName){
1274 m->openInputFile(listFileName, listFile);
1277 listFile >> label >> numOTUs;
1279 otuData.assign(numSeqs, 0);
1280 cumNumSeqs.assign(numOTUs, 0);
1281 nSeqsPerOTU.assign(numOTUs, 0);
1282 aaP.clear();aaP.resize(numOTUs);
1288 string singleOTU = "";
1290 for(int i=0;i<numOTUs;i++){
1292 listFile >> singleOTU;
1294 istringstream otuString(singleOTU);
1298 string seqName = "";
1300 for(int j=0;j<singleOTU.length();j++){
1301 char letter = otuString.get();
1307 map<string,int>::iterator nmIt = nameMap.find(seqName);
1308 int index = nmIt->second;
1310 nameMap.erase(nmIt);
1314 aaP[i].push_back(index);
1319 map<string,int>::iterator nmIt = nameMap.find(seqName);
1321 int index = nmIt->second;
1322 nameMap.erase(nmIt);
1326 aaP[i].push_back(index);
1331 sort(aaP[i].begin(), aaP[i].end());
1332 for(int j=0;j<nSeqsPerOTU[i];j++){
1333 seqNumber.push_back(aaP[i][j]);
1335 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
1336 aaP[i].push_back(0);
1342 for(int i=1;i<numOTUs;i++){
1343 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
1346 seqIndex = seqNumber;
1351 catch(exception& e) {
1352 m->errorOut(e, "ShhherCommand", "getOTUData");
1357 /**************************************************************************************************/
1359 void ShhherCommand::initPyroCluster(){
1361 dist.assign(numSeqs * numOTUs, 0);
1362 change.assign(numOTUs, 1);
1363 centroids.assign(numOTUs, -1);
1364 weight.assign(numOTUs, 0);
1365 singleTau.assign(numSeqs, 1.0);
1367 nSeqsBreaks.assign(processors+1, 0);
1368 nOTUsBreaks.assign(processors+1, 0);
1371 for(int i=0;i<processors;i++){
1372 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
1373 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
1375 nSeqsBreaks[processors] = numSeqs;
1376 nOTUsBreaks[processors] = numOTUs;
1378 catch(exception& e) {
1379 m->errorOut(e, "ShhherCommand", "initPyroCluster");
1384 /**************************************************************************************************/
1386 void ShhherCommand::fill(){
1389 for(int i=0;i<numOTUs;i++){
1390 cumNumSeqs[i] = index;
1391 for(int j=0;j<nSeqsPerOTU[i];j++){
1392 seqNumber[index] = aaP[i][j];
1393 seqIndex[index] = aaI[i][j];
1399 catch(exception& e) {
1400 m->errorOut(e, "ShhherCommand", "fill");
1405 /**************************************************************************************************/
1407 void ShhherCommand::calcCentroids(){
1410 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1412 if(processors == 1) {
1413 calcCentroidsDriver(0, numOTUs);
1415 else{ //you have multiple processors
1416 if (numOTUs < processors){ processors = 1; }
1419 vector<int> processIDs;
1421 //loop through and create all the processes you want
1422 while (process != processors) {
1426 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1428 }else if (pid == 0){
1429 calcCentroidsDriver(nOTUsBreaks[process], nOTUsBreaks[process+1]);
1432 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1434 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1439 //parent does its part
1440 calcCentroidsDriver(nOTUsBreaks[0], nOTUsBreaks[1]);
1442 //force parent to wait until all the processes are done
1443 for (int i=0;i<processIDs.size();i++) {
1444 int temp = processIDs[i];
1450 calcCentroidsDriver(0, numOTUs);
1453 catch(exception& e) {
1454 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1459 /**************************************************************************************************/
1461 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1463 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1469 for(int i=start;i<finish;i++){
1473 int minFlowGram = 100000000;
1474 double minFlowValue = 1e8;
1475 change[i] = 0; //FALSE
1477 for(int j=0;j<nSeqsPerOTU[i];j++){
1478 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1481 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1482 vector<double> adF(nSeqsPerOTU[i]);
1483 vector<int> anL(nSeqsPerOTU[i]);
1485 for(int j=0;j<nSeqsPerOTU[i];j++){
1486 int index = cumNumSeqs[i] + j;
1487 int nI = seqIndex[index];
1488 int nIU = mapSeqToUnique[nI];
1491 for(k=0;k<position;k++){
1497 anL[position] = nIU;
1498 adF[position] = 0.0000;
1503 for(int j=0;j<nSeqsPerOTU[i];j++){
1504 int index = cumNumSeqs[i] + j;
1505 int nI = seqIndex[index];
1507 double tauValue = singleTau[seqNumber[index]];
1509 for(int k=0;k<position;k++){
1510 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1511 adF[k] += dist * tauValue;
1515 for(int j=0;j<position;j++){
1516 if(adF[j] < minFlowValue){
1518 minFlowValue = adF[j];
1522 if(centroids[i] != anL[minFlowGram]){
1524 centroids[i] = anL[minFlowGram];
1527 else if(centroids[i] != -1){
1533 catch(exception& e) {
1534 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1539 /**************************************************************************************************/
1541 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1544 int flowAValue = cent * numFlowCells;
1545 int flowBValue = flow * numFlowCells;
1549 for(int i=0;i<length;i++){
1550 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1555 return dist / (double)length;
1557 catch(exception& e) {
1558 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1563 /**************************************************************************************************/
1565 double ShhherCommand::getNewWeights(){
1568 double maxChange = 0;
1570 for(int i=0;i<numOTUs;i++){
1572 double difference = weight[i];
1575 for(int j=0;j<nSeqsPerOTU[i];j++){
1576 int index = cumNumSeqs[i] + j;
1577 double tauValue = singleTau[seqNumber[index]];
1578 weight[i] += tauValue;
1581 difference = fabs(weight[i] - difference);
1582 if(difference > maxChange){ maxChange = difference; }
1586 catch(exception& e) {
1587 m->errorOut(e, "ShhherCommand", "getNewWeights");
1592 /**************************************************************************************************/
1594 double ShhherCommand::getLikelihood(){
1598 vector<long double> P(numSeqs, 0);
1601 for(int i=0;i<numOTUs;i++){
1602 if(weight[i] > MIN_WEIGHT){
1608 for(int i=0;i<numOTUs;i++){
1609 for(int j=0;j<nSeqsPerOTU[i];j++){
1610 int index = cumNumSeqs[i] + j;
1611 int nI = seqIndex[index];
1612 double singleDist = dist[seqNumber[index]];
1614 P[nI] += weight[i] * exp(-singleDist * sigma);
1618 for(int i=0;i<numSeqs;i++){
1619 if(P[i] == 0){ P[i] = DBL_EPSILON; }
1624 nLL = nLL -(double)numSeqs * log(sigma);
1628 catch(exception& e) {
1629 m->errorOut(e, "ShhherCommand", "getNewWeights");
1634 /**************************************************************************************************/
1636 void ShhherCommand::checkCentroids(){
1638 vector<int> unique(numOTUs, 1);
1640 for(int i=0;i<numOTUs;i++){
1641 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1646 for(int i=0;i<numOTUs;i++){
1648 for(int j=i+1;j<numOTUs;j++){
1651 if(centroids[j] == centroids[i]){
1655 weight[i] += weight[j];
1663 catch(exception& e) {
1664 m->errorOut(e, "ShhherCommand", "checkCentroids");
1669 /**************************************************************************************************/
1671 void ShhherCommand::calcNewDistances(){
1674 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1676 if(processors == 1) {
1677 calcNewDistancesParent(0, numSeqs);
1679 else{ //you have multiple processors
1680 if (numSeqs < processors){ processors = 1; }
1682 vector<vector<int> > child_otuIndex(processors);
1683 vector<vector<int> > child_seqIndex(processors);
1684 vector<vector<double> > child_singleTau(processors);
1685 vector<int> totals(processors);
1688 vector<int> processIDs;
1690 //loop through and create all the processes you want
1691 while (process != processors) {
1695 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1697 }else if (pid == 0){
1698 calcNewDistancesChild(nSeqsBreaks[process], nSeqsBreaks[process+1], child_otuIndex[process], child_seqIndex[process], child_singleTau[process]);
1699 totals[process] = child_otuIndex[process].size();
1703 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1705 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1710 //parent does its part
1711 calcNewDistancesParent(nSeqsBreaks[0], nSeqsBreaks[1]);
1712 int total = seqIndex.size();
1714 //force parent to wait until all the processes are done
1715 for (int i=0;i<processIDs.size();i++) {
1716 int temp = processIDs[i];
1720 for(int i=1;i<processors;i++){
1721 int oldTotal = total;
1724 singleTau.resize(total, 0);
1725 seqIndex.resize(total, 0);
1726 seqNumber.resize(total, 0);
1730 for(int j=oldTotal;j<total;j++){
1731 int otuI = child_otuIndex[i][childIndex];
1732 int seqI = child_seqIndex[i][childIndex];
1734 singleTau[j] = child_singleTau[i][childIndex];
1735 aaP[otuI][nSeqsPerOTU[otuI]] = j;
1736 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
1737 nSeqsPerOTU[otuI]++;
1744 calcNewDistancesParent(0, numSeqs);
1747 catch(exception& e) {
1748 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1753 /**************************************************************************************************/
1755 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1758 vector<double> newTau(numOTUs,0);
1759 vector<double> norms(numSeqs, 0);
1766 for(int i=startSeq;i<stopSeq;i++){
1767 double offset = 1e8;
1768 int indexOffset = i * numOTUs;
1770 for(int j=0;j<numOTUs;j++){
1772 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1773 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1775 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1776 offset = dist[indexOffset + j];
1780 for(int j=0;j<numOTUs;j++){
1781 if(weight[j] > MIN_WEIGHT){
1782 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1783 norms[i] += newTau[j];
1790 for(int j=0;j<numOTUs;j++){
1792 newTau[j] /= norms[i];
1794 if(newTau[j] > MIN_TAU){
1795 otuIndex.push_back(j);
1796 seqIndex.push_back(i);
1797 singleTau.push_back(newTau[j]);
1803 catch(exception& e) {
1804 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1809 /**************************************************************************************************/
1811 void ShhherCommand::calcNewDistancesChild(int startSeq, int stopSeq, vector<int>& child_otuIndex, vector<int>& child_seqIndex, vector<double>& child_singleTau){
1814 vector<double> newTau(numOTUs,0);
1815 vector<double> norms(numSeqs, 0);
1816 child_otuIndex.resize(0);
1817 child_seqIndex.resize(0);
1818 child_singleTau.resize(0);
1820 for(int i=startSeq;i<stopSeq;i++){
1821 double offset = 1e8;
1822 int indexOffset = i * numOTUs;
1825 for(int j=0;j<numOTUs;j++){
1826 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1827 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1830 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1831 offset = dist[indexOffset + j];
1835 for(int j=0;j<numOTUs;j++){
1836 if(weight[j] > MIN_WEIGHT){
1837 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1838 norms[i] += newTau[j];
1845 for(int j=0;j<numOTUs;j++){
1846 newTau[j] /= norms[i];
1848 if(newTau[j] > MIN_TAU){
1849 child_otuIndex.push_back(j);
1850 child_seqIndex.push_back(i);
1851 child_singleTau.push_back(newTau[j]);
1856 catch(exception& e) {
1857 m->errorOut(e, "ShhherCommand", "calcNewDistancesChild");
1862 /**************************************************************************************************/
1864 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1869 vector<double> newTau(numOTUs,0);
1870 vector<double> norms(numSeqs, 0);
1871 nSeqsPerOTU.assign(numOTUs, 0);
1873 for(int i=startSeq;i<stopSeq;i++){
1874 int indexOffset = i * numOTUs;
1876 double offset = 1e8;
1878 for(int j=0;j<numOTUs;j++){
1879 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1880 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1883 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1884 offset = dist[indexOffset + j];
1888 for(int j=0;j<numOTUs;j++){
1889 if(weight[j] > MIN_WEIGHT){
1890 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1891 norms[i] += newTau[j];
1898 for(int j=0;j<numOTUs;j++){
1899 newTau[j] /= norms[i];
1902 for(int j=0;j<numOTUs;j++){
1903 if(newTau[j] > MIN_TAU){
1905 int oldTotal = total;
1909 singleTau.resize(total, 0);
1910 seqNumber.resize(total, 0);
1911 seqIndex.resize(total, 0);
1913 singleTau[oldTotal] = newTau[j];
1915 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1916 aaI[j][nSeqsPerOTU[j]] = i;
1922 catch(exception& e) {
1923 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1928 /**************************************************************************************************/
1930 void ShhherCommand::setOTUs(){
1933 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1935 for(int i=0;i<numOTUs;i++){
1936 for(int j=0;j<nSeqsPerOTU[i];j++){
1937 int index = cumNumSeqs[i] + j;
1938 double tauValue = singleTau[seqNumber[index]];
1939 int sIndex = seqIndex[index];
1940 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1944 for(int i=0;i<numSeqs;i++){
1945 double maxTau = -1.0000;
1947 for(int j=0;j<numOTUs;j++){
1948 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1949 maxTau = bigTauMatrix[i * numOTUs + j];
1954 otuData[i] = maxOTU;
1957 nSeqsPerOTU.assign(numOTUs, 0);
1959 for(int i=0;i<numSeqs;i++){
1960 int index = otuData[i];
1962 singleTau[i] = 1.0000;
1965 aaP[index][nSeqsPerOTU[index]] = i;
1966 aaI[index][nSeqsPerOTU[index]] = i;
1968 nSeqsPerOTU[index]++;
1972 catch(exception& e) {
1973 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1978 /**************************************************************************************************/
1980 void ShhherCommand::writeQualities(vector<int> otuCounts){
1983 string qualityFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.qual";
1985 ofstream qualityFile;
1986 m->openOutputFile(qualityFileName, qualityFile);
1988 qualityFile.setf(ios::fixed, ios::floatfield);
1989 qualityFile.setf(ios::showpoint);
1990 qualityFile << setprecision(6);
1992 vector<vector<int> > qualities(numOTUs);
1993 vector<double> pr(HOMOPS, 0);
1996 for(int i=0;i<numOTUs;i++){
2000 if(nSeqsPerOTU[i] > 0){
2001 qualities[i].assign(1024, -1);
2003 while(index < numFlowCells){
2004 double maxPrValue = 1e8;
2005 short maxPrIndex = -1;
2006 double count = 0.0000;
2008 pr.assign(HOMOPS, 0);
2010 for(int j=0;j<nSeqsPerOTU[i];j++){
2011 int lIndex = cumNumSeqs[i] + j;
2012 double tauValue = singleTau[seqNumber[lIndex]];
2013 int sequenceIndex = aaI[i][j];
2014 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
2018 for(int s=0;s<HOMOPS;s++){
2019 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
2023 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
2024 maxPrValue = pr[maxPrIndex];
2026 if(count > MIN_COUNT){
2028 double norm = 0.0000;
2030 for(int s=0;s<HOMOPS;s++){
2031 norm += exp(-(pr[s] - maxPrValue));
2034 for(int s=1;s<=maxPrIndex;s++){
2036 double temp = 0.0000;
2038 U += exp(-(pr[s-1]-maxPrValue))/norm;
2046 temp = floor(-10 * temp);
2047 value = (int)floor(temp);
2048 if(value > 100){ value = 100; }
2050 qualities[i][base] = (int)value;
2060 if(otuCounts[i] > 0){
2061 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
2063 int j=4; //need to get past the first four bases
2064 while(qualities[i][j] != -1){
2065 qualityFile << qualities[i][j] << ' ';
2068 qualityFile << endl;
2071 qualityFile.close();
2072 outputNames.push_back(qualityFileName);
2075 catch(exception& e) {
2076 m->errorOut(e, "ShhherCommand", "writeQualities");
2081 /**************************************************************************************************/
2083 void ShhherCommand::writeSequences(vector<int> otuCounts){
2086 string fastaFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.fasta";
2088 m->openOutputFile(fastaFileName, fastaFile);
2090 vector<string> names(numOTUs, "");
2092 for(int i=0;i<numOTUs;i++){
2093 int index = centroids[i];
2095 if(otuCounts[i] > 0){
2096 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
2100 for(int j=0;j<numFlowCells;j++){
2102 char base = flowOrder[j % 4];
2103 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
2108 fastaFile << newSeq.substr(4) << endl;
2113 outputNames.push_back(fastaFileName);
2115 if(compositeFASTAFileName != ""){
2116 m->appendFiles(fastaFileName, compositeFASTAFileName);
2119 catch(exception& e) {
2120 m->errorOut(e, "ShhherCommand", "writeSequences");
2125 /**************************************************************************************************/
2127 void ShhherCommand::writeNames(vector<int> otuCounts){
2129 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
2131 m->openOutputFile(nameFileName, nameFile);
2133 for(int i=0;i<numOTUs;i++){
2134 if(otuCounts[i] > 0){
2135 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
2137 for(int j=1;j<nSeqsPerOTU[i];j++){
2138 nameFile << ',' << seqNameVector[aaI[i][j]];
2145 outputNames.push_back(nameFileName);
2148 if(compositeNamesFileName != ""){
2149 m->appendFiles(nameFileName, compositeNamesFileName);
2152 catch(exception& e) {
2153 m->errorOut(e, "ShhherCommand", "writeNames");
2158 /**************************************************************************************************/
2160 void ShhherCommand::writeGroups(){
2162 string fileRoot = flowFileName.substr(0,flowFileName.find_last_of('.'));
2163 string groupFileName = fileRoot + ".shhh.groups";
2165 m->openOutputFile(groupFileName, groupFile);
2167 for(int i=0;i<numSeqs;i++){
2168 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
2171 outputNames.push_back(groupFileName);
2174 catch(exception& e) {
2175 m->errorOut(e, "ShhherCommand", "writeGroups");
2180 /**************************************************************************************************/
2182 void ShhherCommand::writeClusters(vector<int> otuCounts){
2184 string otuCountsFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.counts";
2185 ofstream otuCountsFile;
2186 m->openOutputFile(otuCountsFileName, otuCountsFile);
2188 string bases = flowOrder;
2190 for(int i=0;i<numOTUs;i++){
2191 //output the translated version of the centroid sequence for the otu
2192 if(otuCounts[i] > 0){
2193 int index = centroids[i];
2195 otuCountsFile << "ideal\t";
2196 for(int j=8;j<numFlowCells;j++){
2197 char base = bases[j % 4];
2198 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
2199 otuCountsFile << base;
2202 otuCountsFile << endl;
2204 for(int j=0;j<nSeqsPerOTU[i];j++){
2205 int sequence = aaI[i][j];
2206 otuCountsFile << seqNameVector[sequence] << '\t';
2210 for(int k=0;k<lengths[sequence];k++){
2211 char base = bases[k % 4];
2212 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
2214 for(int s=0;s<freq;s++){
2216 //otuCountsFile << base;
2219 otuCountsFile << newSeq.substr(4) << endl;
2221 otuCountsFile << endl;
2224 otuCountsFile.close();
2225 outputNames.push_back(otuCountsFileName);
2228 catch(exception& e) {
2229 m->errorOut(e, "ShhherCommand", "writeClusters");
2234 //**********************************************************************************************************************