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
30 //**********************************************************************************************************************
32 vector<string> ShhherCommand::getValidParameters(){
35 "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors"
38 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
42 m->errorOut(e, "ShhherCommand", "getValidParameters");
47 //**********************************************************************************************************************
49 ShhherCommand::ShhherCommand(){
51 abort = true; calledHelp = true;
53 //initialize outputTypes
54 vector<string> tempOutNames;
55 outputTypes["pn.dist"] = tempOutNames;
59 m->errorOut(e, "ShhherCommand", "ShhherCommand");
64 //**********************************************************************************************************************
66 vector<string> ShhherCommand::getRequiredParameters(){
68 string Array[] = {"flow"};
69 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
73 m->errorOut(e, "ShhherCommand", "getRequiredParameters");
78 //**********************************************************************************************************************
80 vector<string> ShhherCommand::getRequiredFiles(){
82 vector<string> myArray;
86 m->errorOut(e, "ShhherCommand", "getRequiredFiles");
91 //**********************************************************************************************************************
93 ShhherCommand::ShhherCommand(string option) {
97 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
98 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
104 abort = false; calledHelp = false;
107 //allow user to run help
108 if(option == "help") { help(); abort = true; calledHelp = true; }
112 //valid paramters for this command
113 string AlignArray[] = {
114 "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors"
117 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
119 OptionParser parser(option);
120 map<string,string> parameters = parser.getParameters();
122 ValidParameters validParameter;
123 map<string,string>::iterator it;
125 //check to make sure all parameters are valid for command
126 for (it = parameters.begin(); it != parameters.end(); it++) {
127 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
130 //initialize outputTypes
131 vector<string> tempOutNames;
132 outputTypes["pn.dist"] = tempOutNames;
133 // outputTypes["fasta"] = tempOutNames;
135 //if the user changes the input directory command factory will send this info to us in the output parameter
136 string inputDir = validParameter.validFile(parameters, "inputdir", false);
137 if (inputDir == "not found"){ inputDir = ""; }
140 it = parameters.find("flow");
141 //user has given a template file
142 if(it != parameters.end()){
143 path = m->hasPath(it->second);
144 //if the user has not given a path then, add inputdir. else leave path alone.
145 if (path == "") { parameters["flow"] = inputDir + it->second; }
148 it = parameters.find("lookup");
149 //user has given a template file
150 if(it != parameters.end()){
151 path = m->hasPath(it->second);
152 //if the user has not given a path then, add inputdir. else leave path alone.
153 if (path == "") { parameters["lookup"] = inputDir + it->second; }
156 it = parameters.find("file");
157 //user has given a template file
158 if(it != parameters.end()){
159 path = m->hasPath(it->second);
160 //if the user has not given a path then, add inputdir. else leave path alone.
161 if (path == "") { parameters["file"] = inputDir + it->second; }
166 //check for required parameters
167 flowFileName = validParameter.validFile(parameters, "flow", true);
168 flowFilesFileName = validParameter.validFile(parameters, "file", true);
169 if (flowFileName == "not found" && flowFilesFileName == "not found") {
170 m->mothurOut("values for either flow or file must be provided for the shhh.seqs command.");
171 m->mothurOutEndLine();
174 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
176 //if the user changes the output directory command factory will send this info to us in the output parameter
177 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
179 outputDir += m->hasPath(flowFileName); //if user entered a file with a path then preserve it
183 //check for optional parameter and set defaults
184 // ...at some point should added some additional type checking...
186 temp = validParameter.validFile(parameters, "lookup", true);
187 if (temp == "not found") { lookupFileName = "LookUp_Titanium.pat"; }
188 else if(temp == "not open") { abort = true; }
189 else { lookupFileName = temp; }
191 temp = validParameter.validFile(parameters, "processors", false);if (temp == "not found"){ temp = "1"; }
192 convert(temp, processors);
194 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
195 convert(temp, cutoff);
197 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
198 convert(temp, minDelta);
200 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
201 convert(temp, maxIters);
203 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
204 convert(temp, sigma);
206 globaldata = GlobalData::getInstance();
214 catch(exception& e) {
215 m->errorOut(e, "ShhherCommand", "ShhherCommand");
220 //**********************************************************************************************************************
222 ShhherCommand::~ShhherCommand(){}
224 //**********************************************************************************************************************
226 void ShhherCommand::help(){
228 m->mothurOut("The shhher command reads a file containing flowgrams and creates a file of corrected sequences.\n");
230 catch(exception& e) {
231 m->errorOut(e, "ShhherCommand", "help");
236 //**********************************************************************************************************************
238 int ShhherCommand::execute(){
240 if (abort == true) { if (calledHelp) { return 0; } return 2; }
245 double begClock = clock();
246 unsigned long int begTime = time(NULL);
248 cout.setf(ios::fixed, ios::floatfield);
249 cout.setf(ios::showpoint);
250 cout << setprecision(2);
255 for(int i=1;i<ncpus;i++){
256 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
258 if(abort == 1){ return 0; }
262 cout << "\nGetting preliminary data..." << endl;
266 vector<string> flowFileVector;
267 if(flowFilesFileName != "not found"){
270 ifstream flowFilesFile;
271 m->openInputFile(flowFilesFileName, flowFilesFile);
272 while(flowFilesFile){
273 flowFilesFile >> fName;
274 flowFileVector.push_back(fName);
275 m->gobble(flowFilesFile);
279 flowFileVector.push_back(flowFileName);
281 int numFiles = flowFileVector.size();
283 for(int i=1;i<ncpus;i++){
284 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
287 for(int i=0;i<numFiles;i++){
288 flowFileName = flowFileVector[i];
290 cout << "\n>>>>>\tProcessing " << flowFileName << " (file " << i+1 << " of " << numFiles << ")\t<<<<<" << endl;
291 cout << "Reading flowgrams..." << endl;
293 cout << "Identifying unique flowgrams..." << endl;
296 cout << "Calculating distances between flowgrams..." << endl;
298 strcpy(fileName, flowFileName.c_str());
300 for(int i=1;i<ncpus;i++){
301 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
303 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
304 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
305 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
306 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
307 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
308 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
309 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
310 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
311 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
312 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
315 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
318 for(int i=1;i<ncpus;i++){
319 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
321 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
322 remove((distFileName + ".temp." + toString(i)).c_str());
325 string namesFileName = createNamesFile();
327 cout << "\nClustering flowgrams..." << endl;
328 string listFileName = cluster(distFileName, namesFileName);
329 // string listFileName = "PriestPot_C7.pn.list";
330 // string listFileName = "test.mock_rep3.v69.pn.list";
332 getOTUData(listFileName);
335 for(int i=1;i<ncpus;i++){
336 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
337 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
338 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
339 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
346 int numOTUsOnCPU = numOTUs / ncpus;
347 int numSeqsOnCPU = numSeqs / ncpus;
349 cout << "\nDenoising flowgrams..." << endl;
350 cout << "iter\tmaxDelta\tnLL\t\tcycletime" << endl;
352 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
354 double cycClock = clock();
355 unsigned long int cycTime = time(NULL);
358 int total = singleTau.size();
359 for(int i=1;i<ncpus;i++){
360 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
361 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
362 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
364 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
365 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
366 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
367 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
368 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
371 calcCentroidsDriver(0, numOTUsOnCPU);
373 for(int i=1;i<ncpus;i++){
374 int otuStart = i * numOTUs / ncpus;
375 int otuStop = (i + 1) * numOTUs / ncpus;
377 vector<int> tempCentroids(numOTUs, 0);
378 vector<short> tempChange(numOTUs, 0);
380 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
381 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
383 for(int j=otuStart;j<otuStop;j++){
384 centroids[j] = tempCentroids[j];
385 change[j] = tempChange[j];
389 maxDelta = getNewWeights();
390 double nLL = getLikelihood();
393 for(int i=1;i<ncpus;i++){
394 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
395 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
396 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
399 calcNewDistancesParent(0, numSeqsOnCPU);
401 total = singleTau.size();
403 for(int i=1;i<ncpus;i++){
405 int seqStart = i * numSeqs / ncpus;
406 int seqStop = (i + 1) * numSeqs / ncpus;
408 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
410 vector<int> childSeqIndex(childTotal, 0);
411 vector<double> childSingleTau(childTotal, 0);
412 vector<double> childDist(numSeqs * numOTUs, 0);
413 vector<int> otuIndex(childTotal, 0);
415 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
416 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
417 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
418 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
420 int oldTotal = total;
422 singleTau.resize(total, 0);
423 seqIndex.resize(total, 0);
424 seqNumber.resize(total, 0);
428 for(int j=oldTotal;j<total;j++){
429 int otuI = otuIndex[childIndex];
430 int seqI = childSeqIndex[childIndex];
432 singleTau[j] = childSingleTau[childIndex];
434 aaP[otuI][nSeqsPerOTU[otuI]] = j;
435 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
440 int index = seqStart * numOTUs;
441 for(int j=seqStart;j<seqStop;j++){
442 for(int k=0;k<numOTUs;k++){
443 dist[index] = childDist[index];
451 cout << iter << '\t' << maxDelta << '\t' << setprecision(2) << nLL << '\t' << time(NULL) - cycTime << '\t' << setprecision(6) << (clock() - cycClock)/(double)CLOCKS_PER_SEC << endl;
453 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
455 for(int i=1;i<ncpus;i++){
456 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
461 for(int i=1;i<ncpus;i++){
462 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
468 cout << "\nFinalizing..." << endl;
471 vector<int> otuCounts(numOTUs, 0);
472 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
473 calcCentroidsDriver(0, numOTUs);
474 writeQualities(otuCounts);
475 writeSequences(otuCounts);
476 writeNames(otuCounts);
477 writeClusters(otuCounts);
480 remove(distFileName.c_str());
481 remove(namesFileName.c_str());
482 remove(listFileName.c_str());
484 cout << "Total time to process " << flowFileName << ":\t" << time(NULL) - begTime << '\t' << setprecision(6) << (clock() - begClock)/(double)CLOCKS_PER_SEC << endl;
492 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
493 if(abort){ return 0; }
496 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
498 for(int i=0;i<numFiles;i++){
499 //Now into the pyrodist part
501 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
502 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
503 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
504 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
506 flowDataIntI.resize(numSeqs * numFlowCells);
507 flowDataPrI.resize(numSeqs * numFlowCells);
508 mapUniqueToSeq.resize(numSeqs);
509 mapSeqToUnique.resize(numSeqs);
510 lengths.resize(numSeqs);
511 jointLookUp.resize(NUMBINS * NUMBINS);
513 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
514 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
515 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
516 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
517 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
518 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
519 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
521 flowFileName = string(fileName);
522 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
523 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
525 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
528 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
530 //Now into the pyronoise part
531 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
533 singleLookUp.resize(HOMOPS * NUMBINS);
534 uniqueFlowgrams.resize(numUniques * numFlowCells);
535 weight.resize(numOTUs);
536 centroids.resize(numOTUs);
537 change.resize(numOTUs);
538 dist.assign(numOTUs * numSeqs, 0);
539 nSeqsPerOTU.resize(numOTUs);
540 cumNumSeqs.resize(numOTUs);
542 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
543 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
544 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
546 int startOTU = pid * numOTUs / ncpus;
547 int endOTU = (pid + 1) * numOTUs / ncpus;
549 int startSeq = pid * numSeqs / ncpus;
550 int endSeq = (pid + 1) * numSeqs /ncpus;
556 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
557 singleTau.assign(total, 0.0000);
558 seqNumber.assign(total, 0);
559 seqIndex.assign(total, 0);
561 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
562 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
563 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
564 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
565 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
566 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
567 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
569 calcCentroidsDriver(startOTU, endOTU);
571 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
572 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
575 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
576 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
577 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
579 vector<int> otuIndex(total, 0);
580 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
581 total = otuIndex.size();
583 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
584 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
585 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
586 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
587 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
589 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
594 MPI_Barrier(MPI_COMM_WORLD);
598 catch(exception& e) {
599 m->errorOut(e, "ShhherCommand", "execute");
604 /**************************************************************************************************/
606 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
608 ostringstream outStream;
609 outStream.setf(ios::fixed, ios::floatfield);
610 outStream.setf(ios::dec, ios::basefield);
611 outStream.setf(ios::showpoint);
612 outStream.precision(6);
614 int begTime = time(NULL);
615 double begClock = clock();
617 for(int i=startSeq;i<stopSeq;i++){
618 for(int j=0;j<i;j++){
619 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
621 if(flowDistance < 1e-6){
622 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
624 else if(flowDistance <= cutoff){
625 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
629 cout << i << "\t" << (time(NULL) - begTime) << "\t" << (clock()-begClock)/CLOCKS_PER_SEC << endl;
632 cout << stopSeq << "\t" << (time(NULL) - begTime) << "\t" << (clock()-begClock)/CLOCKS_PER_SEC << endl;
634 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.dist";
635 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
637 ofstream distFile(fDistFileName.c_str());
638 distFile << outStream.str();
641 return fDistFileName;
643 catch(exception& e) {
644 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
650 //**********************************************************************************************************************
652 int ShhherCommand::execute(){
654 if (abort == true) { return 0; }
656 cout.setf(ios::fixed, ios::floatfield);
657 cout.setf(ios::showpoint);
662 vector<string> flowFileVector;
663 if(flowFilesFileName != "not found"){
666 ifstream flowFilesFile;
667 m->openInputFile(flowFilesFileName, flowFilesFile);
668 while(flowFilesFile){
669 flowFilesFile >> fName;
670 flowFileVector.push_back(fName);
671 m->gobble(flowFilesFile);
675 flowFileVector.push_back(flowFileName);
677 int numFiles = flowFileVector.size();
680 for(int i=0;i<numFiles;i++){
681 flowFileName = flowFileVector[i];
683 cout << "\n>>>>>\tProcessing " << flowFileName << " (file " << i+1 << " of " << numFiles << ")\t<<<<<" << endl;
684 cout << "Reading flowgrams..." << endl;
686 cout << "Identifying unique flowgrams..." << endl;
690 cout << "Calculating distances between flowgrams..." << endl;
691 string distFileName = createDistFile(processors);
692 string namesFileName = createNamesFile();
694 cout << "\nClustering flowgrams..." << endl;
695 string listFileName = cluster(distFileName, namesFileName);
696 getOTUData(listFileName);
703 double begClock = clock();
704 unsigned long int begTime = time(NULL);
706 cout << "\nDenoising flowgrams..." << endl;
707 cout << "iter\tmaxDelta\tnLL\t\tcycletime" << endl;
709 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
711 double cycClock = clock();
712 unsigned long int cycTime = time(NULL);
717 maxDelta = getNewWeights();
718 double nLL = getLikelihood();
725 cout << iter << '\t' << maxDelta << '\t' << setprecision(2) << nLL << '\t' << time(NULL) - cycTime << '\t' << setprecision(6) << (clock() - cycClock)/(double)CLOCKS_PER_SEC << endl;
728 cout << "\nFinalizing..." << endl;
732 vector<int> otuCounts(numOTUs, 0);
733 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
735 calcCentroidsDriver(0, numOTUs);
736 writeQualities(otuCounts);
737 writeSequences(otuCounts);
738 writeNames(otuCounts);
739 writeClusters(otuCounts);
742 remove(distFileName.c_str());
743 remove(namesFileName.c_str());
744 remove(listFileName.c_str());
746 cout << "Total time to process " << flowFileName << ":\t" << time(NULL) - begTime << '\t' << setprecision(6) << (clock() - begClock)/(double)CLOCKS_PER_SEC << endl;
750 catch(exception& e) {
751 m->errorOut(e, "ShhherCommand", "execute");
756 /**************************************************************************************************/
758 void ShhherCommand::getFlowData(){
761 m->openInputFile(flowFileName, flowFile);
765 int currentNumFlowCells;
769 flowFile >> numFlowCells;
770 int index = 0;//pcluster
771 while(!flowFile.eof()){
772 flowFile >> seqName >> currentNumFlowCells;
773 lengths.push_back(currentNumFlowCells);
775 seqNameVector.push_back(seqName);
776 nameMap[seqName] = index++;//pcluster
778 for(int i=0;i<numFlowCells;i++){
779 flowFile >> intensity;
780 if(intensity > 9.99) { intensity = 9.99; }
781 int intI = int(100 * intensity + 0.0001);
782 flowDataIntI.push_back(intI);
788 numSeqs = seqNameVector.size();
790 for(int i=0;i<numSeqs;i++){
791 int iNumFlowCells = i * numFlowCells;
792 for(int j=lengths[i];j<numFlowCells;j++){
793 flowDataIntI[iNumFlowCells + j] = 0;
798 catch(exception& e) {
799 m->errorOut(e, "ShhherCommand", "getFlowData");
804 /**************************************************************************************************/
806 void ShhherCommand::getSingleLookUp(){
808 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
809 singleLookUp.assign(HOMOPS * NUMBINS, 0);
813 m->openInputFile(lookupFileName, lookUpFile);
815 for(int i=0;i<HOMOPS;i++){
817 lookUpFile >> logFracFreq;
819 for(int j=0;j<NUMBINS;j++) {
820 lookUpFile >> singleLookUp[index];
826 catch(exception& e) {
827 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
832 /**************************************************************************************************/
834 void ShhherCommand::getJointLookUp(){
837 // the most likely joint probability (-log) that two intenities have the same polymer length
838 jointLookUp.resize(NUMBINS * NUMBINS, 0);
840 for(int i=0;i<NUMBINS;i++){
841 for(int j=0;j<NUMBINS;j++){
843 double minSum = 100000000;
845 for(int k=0;k<HOMOPS;k++){
846 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
848 if(sum < minSum) { minSum = sum; }
850 jointLookUp[i * NUMBINS + j] = minSum;
854 catch(exception& e) {
855 m->errorOut(e, "ShhherCommand", "getJointLookUp");
860 /**************************************************************************************************/
862 double ShhherCommand::getProbIntensity(int intIntensity){
865 double minNegLogProb = 100000000;
868 for(int i=0;i<HOMOPS;i++){//loop signal strength
869 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
870 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
873 return minNegLogProb;
875 catch(exception& e) {
876 m->errorOut(e, "ShhherCommand", "getProbIntensity");
881 /**************************************************************************************************/
883 void ShhherCommand::getUniques(){
888 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
889 uniqueCount.assign(numSeqs, 0); // anWeights
890 uniqueLengths.assign(numSeqs, 0);
891 mapSeqToUnique.assign(numSeqs, -1);
892 mapUniqueToSeq.assign(numSeqs, -1);
894 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
896 for(int i=0;i<numSeqs;i++){
899 vector<short> current(numFlowCells);
900 for(int j=0;j<numFlowCells;j++){ current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0)); }
902 for(int j=0;j<numUniques;j++){
903 int offset = j * numFlowCells;
906 for(int k=0;k<numFlowCells;k++){
907 if(current[k] != uniqueFlowgrams[offset + k]){
914 mapSeqToUnique[i] = j;
922 if(index == numUniques){
923 uniqueLengths[numUniques] = lengths[i];
924 uniqueCount[numUniques] = 1;
925 mapSeqToUnique[i] = numUniques;//anMap
926 mapUniqueToSeq[numUniques] = i;//anF
928 for(int k=0;k<numFlowCells;k++){
929 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
930 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
936 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
937 uniqueLengths.resize(numUniques);
939 flowDataPrI.assign(numSeqs * numFlowCells, 0);
940 for(int i=0;i<flowDataPrI.size();i++) { flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
942 catch(exception& e) {
943 m->errorOut(e, "ShhherCommand", "getUniques");
948 /**************************************************************************************************/
950 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
952 int minLength = lengths[mapSeqToUnique[seqA]];
953 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
955 int ANumFlowCells = seqA * numFlowCells;
956 int BNumFlowCells = seqB * numFlowCells;
960 for(int i=0;i<minLength;i++){
961 int flowAIntI = flowDataIntI[ANumFlowCells + i];
962 float flowAPrI = flowDataPrI[ANumFlowCells + i];
964 int flowBIntI = flowDataIntI[BNumFlowCells + i];
965 float flowBPrI = flowDataPrI[BNumFlowCells + i];
966 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
969 dist /= (float) minLength;
972 catch(exception& e) {
973 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
978 /**************************************************************************************************/
980 void ShhherCommand::flowDistParentFork(string distFileName, int startSeq, int stopSeq){
983 ostringstream outStream;
984 outStream.setf(ios::fixed, ios::floatfield);
985 outStream.setf(ios::dec, ios::basefield);
986 outStream.setf(ios::showpoint);
987 outStream.precision(6);
989 int begTime = time(NULL);
990 double begClock = clock();
992 for(int i=startSeq;i<stopSeq;i++){
993 for(int j=0;j<i;j++){
994 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
996 if(flowDistance < 1e-6){
997 outStream << seqNameVector[mapUniqueToSeq[i]] << '\t' << seqNameVector[mapUniqueToSeq[j]] << '\t' << 0.000000 << endl;
999 else if(flowDistance <= cutoff){
1000 outStream << seqNameVector[mapUniqueToSeq[i]] << '\t' << seqNameVector[mapUniqueToSeq[j]] << '\t' << flowDistance << endl;
1004 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
1005 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1006 m->mothurOutEndLine();
1009 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
1010 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1011 m->mothurOutEndLine();
1013 ofstream distFile(distFileName.c_str());
1014 distFile << outStream.str();
1017 catch(exception& e) {
1018 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
1023 /**************************************************************************************************/
1025 string ShhherCommand::createDistFile(int processors){
1027 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.dist";
1029 unsigned long int begTime = time(NULL);
1030 double begClock = clock();
1035 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1036 if(processors == 1) { flowDistParentFork(fDistFileName, 0, numUniques); }
1037 else{ //you have multiple processors
1039 if (numSeqs < processors){ processors = 1; }
1041 vector<int> start(processors, 0);
1042 vector<int> end(processors, 0);
1044 for (int i = 0; i < processors; i++) {
1045 start[i] = int(sqrt(float(i)/float(processors)) * numUniques);
1046 end[i] = int(sqrt(float(i+1)/float(processors)) * numUniques);
1050 vector<int> processIDs;
1052 //loop through and create all the processes you want
1053 while (process != processors) {
1057 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1059 }else if (pid == 0){
1060 flowDistParentFork(fDistFileName + toString(getpid()) + ".temp", start[process], end[process]);
1063 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1065 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1070 //parent does its part
1071 flowDistParentFork(fDistFileName, start[0], end[0]);
1073 //force parent to wait until all the processes are done
1074 for (int i=0;i<processIDs.size();i++) {
1075 int temp = processIDs[i];
1079 //append and remove temp files
1080 for (int i=0;i<processIDs.size();i++) {
1081 m->appendFiles((fDistFileName + toString(processIDs[i]) + ".temp"), fDistFileName);
1082 remove((fDistFileName + toString(processIDs[i]) + ".temp").c_str());
1088 flowDistParentFork(fDistFileName, 0, numUniques);
1091 m->mothurOutEndLine();
1093 cout << "Total time: " << (time(NULL) - begTime) << "\t" << (clock() - begClock)/CLOCKS_PER_SEC << endl;;
1095 return fDistFileName;
1097 catch(exception& e) {
1098 m->errorOut(e, "ShhherCommand", "createDistFile");
1104 /**************************************************************************************************/
1106 string ShhherCommand::createNamesFile(){
1109 vector<string> duplicateNames(numUniques, "");
1110 for(int i=0;i<numSeqs;i++){
1111 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
1114 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.names";
1117 m->openOutputFile(nameFileName, nameFile);
1119 for(int i=0;i<numUniques;i++){
1120 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1121 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1125 return nameFileName;
1127 catch(exception& e) {
1128 m->errorOut(e, "ShhherCommand", "createNamesFile");
1133 //**********************************************************************************************************************
1135 string ShhherCommand::cluster(string distFileName, string namesFileName){
1138 SparseMatrix* matrix;
1140 RAbundVector* rabund;
1142 globaldata->setNameFile(namesFileName);
1143 globaldata->setColumnFile(distFileName);
1144 globaldata->setFormat("column");
1146 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1147 read->setCutoff(cutoff);
1149 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1150 clusterNameMap->readMap();
1151 read->read(clusterNameMap);
1153 list = read->getListVector();
1154 matrix = read->getMatrix();
1157 delete clusterNameMap;
1159 rabund = new RAbundVector(list->getRAbundVector());
1161 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1162 string tag = cluster->getTag();
1164 double clusterCutoff = cutoff;
1165 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1166 cluster->update(clusterCutoff);
1169 list->setLabel(toString(cutoff));
1171 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.list";
1173 m->openOutputFile(listFileName, listFile);
1174 list->print(listFile);
1177 delete matrix; delete cluster; delete rabund; delete list;
1179 return listFileName;
1181 catch(exception& e) {
1182 m->errorOut(e, "ShhherCommand", "cluster");
1187 /**************************************************************************************************/
1189 void ShhherCommand::getOTUData(string listFileName){
1193 m->openInputFile(listFileName, listFile);
1196 listFile >> label >> numOTUs;
1198 otuData.assign(numSeqs, 0);
1199 cumNumSeqs.assign(numOTUs, 0);
1200 nSeqsPerOTU.assign(numOTUs, 0);
1201 aaP.resize(numOTUs);
1203 string singleOTU = "";
1205 for(int i=0;i<numOTUs;i++){
1207 listFile >> singleOTU;
1209 istringstream otuString(singleOTU);
1213 string seqName = "";
1215 for(int j=0;j<singleOTU.length();j++){
1216 char letter = otuString.get();
1222 map<string,int>::iterator nmIt = nameMap.find(seqName);
1223 int index = nmIt->second;
1225 nameMap.erase(nmIt);
1229 aaP[i].push_back(index);
1234 map<string,int>::iterator nmIt = nameMap.find(seqName);
1236 int index = nmIt->second;
1237 nameMap.erase(nmIt);
1241 aaP[i].push_back(index);
1246 sort(aaP[i].begin(), aaP[i].end());
1247 for(int j=0;j<nSeqsPerOTU[i];j++){
1248 seqNumber.push_back(aaP[i][j]);
1250 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
1251 aaP[i].push_back(0);
1255 for(int i=1;i<numOTUs;i++){
1256 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
1259 seqIndex = seqNumber;
1263 catch(exception& e) {
1264 m->errorOut(e, "ShhherCommand", "getOTUData");
1269 /**************************************************************************************************/
1271 void ShhherCommand::initPyroCluster(){
1273 dist.assign(numSeqs * numOTUs, 0);
1274 change.assign(numOTUs, 1);
1275 centroids.assign(numOTUs, -1);
1276 weight.assign(numOTUs, 0);
1277 singleTau.assign(numSeqs, 1.0);
1279 nSeqsBreaks.assign(processors+1, 0);
1280 nOTUsBreaks.assign(processors+1, 0);
1283 for(int i=0;i<processors;i++){
1284 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
1285 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
1287 nSeqsBreaks[processors] = numSeqs;
1288 nOTUsBreaks[processors] = numOTUs;
1290 catch(exception& e) {
1291 m->errorOut(e, "ShhherCommand", "initPyroCluster");
1296 /**************************************************************************************************/
1298 void ShhherCommand::fill(){
1301 for(int i=0;i<numOTUs;i++){
1302 cumNumSeqs[i] = index;
1303 for(int j=0;j<nSeqsPerOTU[i];j++){
1304 seqNumber[index] = aaP[i][j];
1305 seqIndex[index] = aaI[i][j];
1311 catch(exception& e) {
1312 m->errorOut(e, "ShhherCommand", "fill");
1317 /**************************************************************************************************/
1319 void ShhherCommand::calcCentroids(){
1322 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1324 if(processors == 1) {
1325 calcCentroidsDriver(0, numOTUs);
1327 else{ //you have multiple processors
1328 if (numOTUs < processors){ processors = 1; }
1331 vector<int> processIDs;
1333 //loop through and create all the processes you want
1334 while (process != processors) {
1338 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1340 }else if (pid == 0){
1341 calcCentroidsDriver(nOTUsBreaks[process], nOTUsBreaks[process+1]);
1344 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1346 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1351 //parent does its part
1352 calcCentroidsDriver(nOTUsBreaks[0], nOTUsBreaks[1]);
1354 //force parent to wait until all the processes are done
1355 for (int i=0;i<processIDs.size();i++) {
1356 int temp = processIDs[i];
1362 calcCentroidsDriver(0, numOTUs);
1365 catch(exception& e) {
1366 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1371 /**************************************************************************************************/
1373 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1375 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1380 for(int i=start;i<finish;i++){
1384 int minFlowGram = 100000000;
1385 double minFlowValue = 1e8;
1386 change[i] = 0; //FALSE
1388 for(int j=0;j<nSeqsPerOTU[i];j++){
1389 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1392 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1393 vector<double> adF(nSeqsPerOTU[i]);
1394 vector<int> anL(nSeqsPerOTU[i]);
1396 for(int j=0;j<nSeqsPerOTU[i];j++){
1397 int index = cumNumSeqs[i] + j;
1398 int nI = seqIndex[index];
1399 int nIU = mapSeqToUnique[nI];
1402 for(k=0;k<position;k++){
1408 anL[position] = nIU;
1409 adF[position] = 0.0000;
1414 for(int j=0;j<nSeqsPerOTU[i];j++){
1415 int index = cumNumSeqs[i] + j;
1416 int nI = seqIndex[index];
1418 double tauValue = singleTau[seqNumber[index]];
1420 for(int k=0;k<position;k++){
1421 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1422 adF[k] += dist * tauValue;
1426 for(int j=0;j<position;j++){
1427 if(adF[j] < minFlowValue){
1429 minFlowValue = adF[j];
1433 if(centroids[i] != anL[minFlowGram]){
1435 centroids[i] = anL[minFlowGram];
1438 else if(centroids[i] != -1){
1444 catch(exception& e) {
1445 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1450 /**************************************************************************************************/
1452 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1455 int flowAValue = cent * numFlowCells;
1456 int flowBValue = flow * numFlowCells;
1460 for(int i=0;i<length;i++){
1461 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1466 return dist / (double)length;
1468 catch(exception& e) {
1469 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1474 /**************************************************************************************************/
1476 double ShhherCommand::getNewWeights(){
1479 double maxChange = 0;
1481 for(int i=0;i<numOTUs;i++){
1483 double difference = weight[i];
1486 for(int j=0;j<nSeqsPerOTU[i];j++){
1487 int index = cumNumSeqs[i] + j;
1488 double tauValue = singleTau[seqNumber[index]];
1489 weight[i] += tauValue;
1492 difference = fabs(weight[i] - difference);
1493 if(difference > maxChange){ maxChange = difference; }
1497 catch(exception& e) {
1498 m->errorOut(e, "ShhherCommand", "getNewWeights");
1503 /**************************************************************************************************/
1505 double ShhherCommand::getLikelihood(){
1509 vector<long double> P(numSeqs, 0);
1512 for(int i=0;i<numOTUs;i++){
1513 if(weight[i] > MIN_WEIGHT){
1519 for(int i=0;i<numOTUs;i++){
1520 for(int j=0;j<nSeqsPerOTU[i];j++){
1521 int index = cumNumSeqs[i] + j;
1522 int nI = seqIndex[index];
1523 double singleDist = dist[seqNumber[index]];
1525 P[nI] += weight[i] * exp(-singleDist * sigma);
1529 for(int i=0;i<numSeqs;i++){
1530 if(P[i] == 0){ P[i] = DBL_EPSILON; }
1535 nLL = nLL -(double)numSeqs * log(sigma);
1539 catch(exception& e) {
1540 m->errorOut(e, "ShhherCommand", "getNewWeights");
1545 /**************************************************************************************************/
1547 void ShhherCommand::checkCentroids(){
1549 vector<int> unique(numOTUs, 1);
1551 for(int i=0;i<numOTUs;i++){
1552 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1557 for(int i=0;i<numOTUs;i++){
1559 for(int j=i+1;j<numOTUs;j++){
1562 if(centroids[j] == centroids[i]){
1566 weight[i] += weight[j];
1574 catch(exception& e) {
1575 m->errorOut(e, "ShhherCommand", "checkCentroids");
1580 /**************************************************************************************************/
1582 void ShhherCommand::calcNewDistances(){
1585 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1587 if(processors == 1) {
1588 calcNewDistancesParent(0, numSeqs);
1590 else{ //you have multiple processors
1591 if (numSeqs < processors){ processors = 1; }
1593 vector<vector<int> > child_otuIndex(processors);
1594 vector<vector<int> > child_seqIndex(processors);
1595 vector<vector<double> > child_singleTau(processors);
1596 vector<int> totals(processors);
1599 vector<int> processIDs;
1601 //loop through and create all the processes you want
1602 while (process != processors) {
1606 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1608 }else if (pid == 0){
1609 calcNewDistancesChild(nSeqsBreaks[process], nSeqsBreaks[process+1], child_otuIndex[process], child_seqIndex[process], child_singleTau[process]);
1610 totals[process] = child_otuIndex[process].size();
1614 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1616 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1621 //parent does its part
1622 calcNewDistancesParent(nSeqsBreaks[0], nSeqsBreaks[1]);
1623 int total = seqIndex.size();
1625 //force parent to wait until all the processes are done
1626 for (int i=0;i<processIDs.size();i++) {
1627 int temp = processIDs[i];
1631 for(int i=1;i<processors;i++){
1632 int oldTotal = total;
1635 singleTau.resize(total, 0);
1636 seqIndex.resize(total, 0);
1637 seqNumber.resize(total, 0);
1641 for(int j=oldTotal;j<total;j++){
1642 int otuI = child_otuIndex[i][childIndex];
1643 int seqI = child_seqIndex[i][childIndex];
1645 singleTau[j] = child_singleTau[i][childIndex];
1646 aaP[otuI][nSeqsPerOTU[otuI]] = j;
1647 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
1648 nSeqsPerOTU[otuI]++;
1655 calcNewDistancesParent(0, numSeqs);
1658 catch(exception& e) {
1659 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1664 /**************************************************************************************************/
1666 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1669 vector<double> newTau(numOTUs,0);
1670 vector<double> norms(numSeqs, 0);
1673 singleTau.resize(0);
1677 for(int i=startSeq;i<stopSeq;i++){
1678 double offset = 1e8;
1679 int indexOffset = i * numOTUs;
1681 for(int j=0;j<numOTUs;j++){
1683 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1684 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1686 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1687 offset = dist[indexOffset + j];
1691 for(int j=0;j<numOTUs;j++){
1692 if(weight[j] > MIN_WEIGHT){
1693 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1694 norms[i] += newTau[j];
1701 for(int j=0;j<numOTUs;j++){
1703 newTau[j] /= norms[i];
1705 if(newTau[j] > MIN_TAU){
1706 otuIndex.push_back(j);
1707 seqIndex.push_back(i);
1708 singleTau.push_back(newTau[j]);
1714 catch(exception& e) {
1715 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1720 /**************************************************************************************************/
1722 void ShhherCommand::calcNewDistancesChild(int startSeq, int stopSeq, vector<int>& child_otuIndex, vector<int>& child_seqIndex, vector<double>& child_singleTau){
1725 vector<double> newTau(numOTUs,0);
1726 vector<double> norms(numSeqs, 0);
1727 child_otuIndex.resize(0);
1728 child_seqIndex.resize(0);
1729 child_singleTau.resize(0);
1731 for(int i=startSeq;i<stopSeq;i++){
1732 double offset = 1e8;
1733 int indexOffset = i * numOTUs;
1736 for(int j=0;j<numOTUs;j++){
1737 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1738 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1741 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1742 offset = dist[indexOffset + j];
1746 for(int j=0;j<numOTUs;j++){
1747 if(weight[j] > MIN_WEIGHT){
1748 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1749 norms[i] += newTau[j];
1756 for(int j=0;j<numOTUs;j++){
1757 newTau[j] /= norms[i];
1759 if(newTau[j] > MIN_TAU){
1760 child_otuIndex.push_back(j);
1761 child_seqIndex.push_back(i);
1762 child_singleTau.push_back(newTau[j]);
1767 catch(exception& e) {
1768 m->errorOut(e, "ShhherCommand", "calcNewDistancesChild");
1773 /**************************************************************************************************/
1775 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1780 vector<double> newTau(numOTUs,0);
1781 vector<double> norms(numSeqs, 0);
1782 nSeqsPerOTU.assign(numOTUs, 0);
1784 for(int i=startSeq;i<stopSeq;i++){
1785 int indexOffset = i * numOTUs;
1787 double offset = 1e8;
1789 for(int j=0;j<numOTUs;j++){
1790 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1791 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1794 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1795 offset = dist[indexOffset + j];
1799 for(int j=0;j<numOTUs;j++){
1800 if(weight[j] > MIN_WEIGHT){
1801 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1802 norms[i] += newTau[j];
1809 for(int j=0;j<numOTUs;j++){
1810 newTau[j] /= norms[i];
1813 for(int j=0;j<numOTUs;j++){
1814 if(newTau[j] > MIN_TAU){
1816 int oldTotal = total;
1820 singleTau.resize(total, 0);
1821 seqNumber.resize(total, 0);
1822 seqIndex.resize(total, 0);
1824 singleTau[oldTotal] = newTau[j];
1826 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1827 aaI[j][nSeqsPerOTU[j]] = i;
1833 catch(exception& e) {
1834 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1839 /**************************************************************************************************/
1841 void ShhherCommand::setOTUs(){
1844 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1846 for(int i=0;i<numOTUs;i++){
1847 for(int j=0;j<nSeqsPerOTU[i];j++){
1848 int index = cumNumSeqs[i] + j;
1849 double tauValue = singleTau[seqNumber[index]];
1850 int sIndex = seqIndex[index];
1851 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1855 for(int i=0;i<numSeqs;i++){
1856 double maxTau = -1.0000;
1858 for(int j=0;j<numOTUs;j++){
1859 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1860 maxTau = bigTauMatrix[i * numOTUs + j];
1865 otuData[i] = maxOTU;
1868 nSeqsPerOTU.assign(numOTUs, 0);
1870 for(int i=0;i<numSeqs;i++){
1871 int index = otuData[i];
1873 singleTau[i] = 1.0000;
1876 aaP[index][nSeqsPerOTU[index]] = i;
1877 aaI[index][nSeqsPerOTU[index]] = i;
1879 nSeqsPerOTU[index]++;
1883 catch(exception& e) {
1884 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1889 /**************************************************************************************************/
1891 void ShhherCommand::writeQualities(vector<int> otuCounts){
1894 string qualityFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.qual";
1896 ofstream qualityFile;
1897 m->openOutputFile(qualityFileName, qualityFile);
1899 qualityFile.setf(ios::fixed, ios::floatfield);
1900 qualityFile.setf(ios::showpoint);
1901 qualityFile << setprecision(6);
1903 vector<vector<int> > qualities(numOTUs);
1904 vector<double> pr(HOMOPS, 0);
1907 for(int i=0;i<numOTUs;i++){
1911 if(nSeqsPerOTU[i] > 0){
1912 qualities[i].assign(1024, -1);
1914 while(index < numFlowCells){
1915 double maxPrValue = 1e8;
1916 short maxPrIndex = -1;
1917 double count = 0.0000;
1919 pr.assign(HOMOPS, 0);
1921 for(int j=0;j<nSeqsPerOTU[i];j++){
1922 int lIndex = cumNumSeqs[i] + j;
1923 double tauValue = singleTau[seqNumber[lIndex]];
1924 int sequenceIndex = aaI[i][j];
1925 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1929 for(int s=0;s<HOMOPS;s++){
1930 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1934 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1935 maxPrValue = pr[maxPrIndex];
1937 if(count > MIN_COUNT){
1939 double norm = 0.0000;
1941 for(int s=0;s<HOMOPS;s++){
1942 norm += exp(-(pr[s] - maxPrValue));
1945 for(int s=1;s<=maxPrIndex;s++){
1947 double temp = 0.0000;
1949 U += exp(-(pr[s-1]-maxPrValue))/norm;
1957 temp = floor(-10 * temp);
1958 value = (int)floor(temp);
1959 if(value > 100){ value = 100; }
1961 qualities[i][base] = (int)value;
1971 if(otuCounts[i] > 0){
1972 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1974 int j=4; //need to get past the first four bases
1975 while(qualities[i][j] != -1){
1976 qualityFile << qualities[i][j] << ' ';
1979 qualityFile << endl;
1982 qualityFile.close();
1985 catch(exception& e) {
1986 m->errorOut(e, "ShhherCommand", "writeQualities");
1991 /**************************************************************************************************/
1993 void ShhherCommand::writeSequences(vector<int> otuCounts){
1996 string bases = "TACG";
1998 string fastaFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.fasta";
2000 m->openOutputFile(fastaFileName, fastaFile);
2002 vector<string> names(numOTUs, "");
2004 for(int i=0;i<numOTUs;i++){
2005 int index = centroids[i];
2007 if(otuCounts[i] > 0){
2008 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
2010 for(int j=8;j<numFlowCells;j++){
2012 char base = bases[j % 4];
2013 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
2022 catch(exception& e) {
2023 m->errorOut(e, "ShhherCommand", "writeSequences");
2028 /**************************************************************************************************/
2030 void ShhherCommand::writeNames(vector<int> otuCounts){
2032 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.final.names";
2034 m->openOutputFile(nameFileName, nameFile);
2036 for(int i=0;i<numOTUs;i++){
2037 if(otuCounts[i] > 0){
2038 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
2040 for(int j=1;j<nSeqsPerOTU[i];j++){
2041 nameFile << ',' << seqNameVector[aaI[i][j]];
2049 catch(exception& e) {
2050 m->errorOut(e, "ShhherCommand", "writeNames");
2055 /**************************************************************************************************/
2057 void ShhherCommand::writeGroups(){
2059 string fileRoot = flowFileName.substr(0,flowFileName.find_last_of('.'));
2060 string groupFileName = fileRoot + ".pn.groups";
2062 m->openOutputFile(groupFileName, groupFile);
2064 for(int i=0;i<numSeqs;i++){
2065 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
2069 catch(exception& e) {
2070 m->errorOut(e, "ShhherCommand", "writeGroups");
2075 /**************************************************************************************************/
2077 void ShhherCommand::writeClusters(vector<int> otuCounts){
2079 string otuCountsFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.counts";
2080 ofstream otuCountsFile;
2081 m->openOutputFile(otuCountsFileName, otuCountsFile);
2083 string bases = "TACG";
2085 for(int i=0;i<numOTUs;i++){
2086 //output the translated version of the centroid sequence for the otu
2087 if(otuCounts[i] > 0){
2088 int index = centroids[i];
2090 otuCountsFile << "ideal\t";
2091 for(int j=8;j<numFlowCells;j++){
2092 char base = bases[j % 4];
2093 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
2094 otuCountsFile << base;
2097 otuCountsFile << endl;
2099 for(int j=0;j<nSeqsPerOTU[i];j++){
2100 int sequence = aaI[i][j];
2101 otuCountsFile << seqNameVector[sequence] << '\t';
2103 for(int k=8;k<lengths[sequence];k++){
2104 char base = bases[k % 4];
2105 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
2107 for(int s=0;s<freq;s++){
2108 otuCountsFile << base;
2111 otuCountsFile << endl;
2113 otuCountsFile << endl;
2116 otuCountsFile.close();
2118 catch(exception& e) {
2119 m->errorOut(e, "ShhherCommand", "writeClusters");
2124 //**********************************************************************************************************************