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"
20 //**********************************************************************************************************************
25 #define MIN_WEIGHT 0.1
26 #define MIN_TAU 0.0001
29 //**********************************************************************************************************************
31 vector<string> ShhherCommand::getValidParameters(){
34 "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors"
37 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
41 m->errorOut(e, "ShhherCommand", "getValidParameters");
46 //**********************************************************************************************************************
48 ShhherCommand::ShhherCommand(){
52 //initialize outputTypes
53 vector<string> tempOutNames;
54 outputTypes["pn.dist"] = tempOutNames;
58 m->errorOut(e, "ShhherCommand", "ShhherCommand");
63 //**********************************************************************************************************************
65 vector<string> ShhherCommand::getRequiredParameters(){
67 string Array[] = {"flow"};
68 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
72 m->errorOut(e, "ShhherCommand", "getRequiredParameters");
77 //**********************************************************************************************************************
79 vector<string> ShhherCommand::getRequiredFiles(){
81 vector<string> myArray;
85 m->errorOut(e, "ShhherCommand", "getRequiredFiles");
90 //**********************************************************************************************************************
92 ShhherCommand::ShhherCommand(string option) {
96 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
97 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
106 //allow user to run help
107 if(option == "help") { help(); abort = true; }
111 //valid paramters for this command
112 string AlignArray[] = {
113 "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors"
116 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
118 OptionParser parser(option);
119 map<string,string> parameters = parser.getParameters();
121 ValidParameters validParameter;
122 map<string,string>::iterator it;
124 //check to make sure all parameters are valid for command
125 for (it = parameters.begin(); it != parameters.end(); it++) {
126 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
129 //initialize outputTypes
130 vector<string> tempOutNames;
131 outputTypes["pn.dist"] = tempOutNames;
132 // outputTypes["fasta"] = tempOutNames;
134 //if the user changes the input directory command factory will send this info to us in the output parameter
135 string inputDir = validParameter.validFile(parameters, "inputdir", false);
136 if (inputDir == "not found"){ inputDir = ""; }
139 it = parameters.find("flow");
140 //user has given a template file
141 if(it != parameters.end()){
142 path = m->hasPath(it->second);
143 //if the user has not given a path then, add inputdir. else leave path alone.
144 if (path == "") { parameters["flow"] = inputDir + it->second; }
147 it = parameters.find("lookup");
148 //user has given a template file
149 if(it != parameters.end()){
150 path = m->hasPath(it->second);
151 //if the user has not given a path then, add inputdir. else leave path alone.
152 if (path == "") { parameters["lookup"] = inputDir + it->second; }
155 it = parameters.find("file");
156 //user has given a template file
157 if(it != parameters.end()){
158 path = m->hasPath(it->second);
159 //if the user has not given a path then, add inputdir. else leave path alone.
160 if (path == "") { parameters["file"] = inputDir + it->second; }
165 //check for required parameters
166 flowFileName = validParameter.validFile(parameters, "flow", true);
167 flowFilesFileName = validParameter.validFile(parameters, "file", true);
168 if (flowFileName == "not found" && flowFilesFileName == "not found") {
169 m->mothurOut("values for either flow or file must be provided for the shhh.seqs command.");
170 m->mothurOutEndLine();
173 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
175 //if the user changes the output directory command factory will send this info to us in the output parameter
176 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
178 outputDir += m->hasPath(flowFileName); //if user entered a file with a path then preserve it
182 //check for optional parameter and set defaults
183 // ...at some point should added some additional type checking...
185 temp = validParameter.validFile(parameters, "lookup", true);
186 if (temp == "not found") { lookupFileName = "LookUp_E123.pat"; }
187 else if(temp == "not open") { abort = true; }
188 else { lookupFileName = temp; }
190 temp = validParameter.validFile(parameters, "processors", false);if (temp == "not found"){ temp = "1"; }
191 convert(temp, processors);
193 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
194 convert(temp, cutoff);
196 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
197 convert(temp, minDelta);
199 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
200 convert(temp, maxIters);
202 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
203 convert(temp, sigma);
205 globaldata = GlobalData::getInstance();
213 catch(exception& e) {
214 m->errorOut(e, "ShhherCommand", "ShhherCommand");
219 //**********************************************************************************************************************
221 ShhherCommand::~ShhherCommand(){}
223 //**********************************************************************************************************************
225 void ShhherCommand::help(){
227 m->mothurOut("The shhher command reads a file containing flowgrams and creates a file of corrected sequences.\n");
229 catch(exception& e) {
230 m->errorOut(e, "ShhherCommand", "help");
235 //**********************************************************************************************************************
237 int ShhherCommand::execute(){
242 double begClock = clock();
243 unsigned long int begTime = time(NULL);
245 cout.setf(ios::fixed, ios::floatfield);
246 cout.setf(ios::showpoint);
247 cout << setprecision(2);
252 for(int i=1;i<ncpus;i++){
253 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
255 if(abort == 1){ return 0; }
259 cout << "\nGetting preliminary data..." << endl;
263 vector<string> flowFileVector;
264 if(flowFilesFileName != "not found"){
267 ifstream flowFilesFile;
268 m->openInputFile(flowFilesFileName, flowFilesFile);
269 while(flowFilesFile){
270 flowFilesFile >> fName;
271 flowFileVector.push_back(fName);
272 m->gobble(flowFilesFile);
276 flowFileVector.push_back(flowFileName);
278 int numFiles = flowFileVector.size();
280 for(int i=1;i<ncpus;i++){
281 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
284 for(int i=0;i<numFiles;i++){
285 flowFileName = flowFileVector[i];
287 cout << "\n>>>>>\tProcessing " << flowFileName << " (file " << i+1 << " of " << numFiles << ")\t<<<<<" << endl;
288 cout << "Reading flowgrams..." << endl;
290 cout << "Identifying unique flowgrams..." << endl;
293 cout << "Calculating distances between flowgrams..." << endl;
295 strcpy(fileName, flowFileName.c_str());
297 for(int i=1;i<ncpus;i++){
298 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
300 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
301 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
302 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
303 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
304 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
305 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
306 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
307 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
308 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
309 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
312 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
315 for(int i=1;i<ncpus;i++){
316 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
318 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
319 remove((distFileName + ".temp." + toString(i)).c_str());
322 string namesFileName = createNamesFile();
324 cout << "\nClustering flowgrams..." << endl;
325 string listFileName = cluster(distFileName, namesFileName);
326 // string listFileName = "PriestPot_C7.pn.list";
327 // string listFileName = "test.mock_rep3.v69.pn.list";
329 getOTUData(listFileName);
332 for(int i=1;i<ncpus;i++){
333 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
334 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
335 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
336 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
343 int numOTUsOnCPU = numOTUs / ncpus;
344 int numSeqsOnCPU = numSeqs / ncpus;
346 cout << "\nDenoising flowgrams..." << endl;
347 cout << "iter\tmaxDelta\tnLL\t\tcycletime" << endl;
349 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
351 double cycClock = clock();
352 unsigned long int cycTime = time(NULL);
355 int total = singleTau.size();
356 for(int i=1;i<ncpus;i++){
357 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
358 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
359 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
361 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
362 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
363 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
364 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
365 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
368 calcCentroidsDriver(0, numOTUsOnCPU);
370 for(int i=1;i<ncpus;i++){
371 int otuStart = i * numOTUs / ncpus;
372 int otuStop = (i + 1) * numOTUs / ncpus;
374 vector<int> tempCentroids(numOTUs, 0);
375 vector<short> tempChange(numOTUs, 0);
377 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
378 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
380 for(int j=otuStart;j<otuStop;j++){
381 centroids[j] = tempCentroids[j];
382 change[j] = tempChange[j];
386 maxDelta = getNewWeights();
387 double nLL = getLikelihood();
390 for(int i=1;i<ncpus;i++){
391 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
392 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
393 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
396 calcNewDistancesParent(0, numSeqsOnCPU);
398 total = singleTau.size();
400 for(int i=1;i<ncpus;i++){
402 int seqStart = i * numSeqs / ncpus;
403 int seqStop = (i + 1) * numSeqs / ncpus;
405 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
407 vector<int> childSeqIndex(childTotal, 0);
408 vector<double> childSingleTau(childTotal, 0);
409 vector<double> childDist(numSeqs * numOTUs, 0);
410 vector<int> otuIndex(childTotal, 0);
412 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
413 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
414 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
415 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
417 int oldTotal = total;
419 singleTau.resize(total, 0);
420 seqIndex.resize(total, 0);
421 seqNumber.resize(total, 0);
425 for(int j=oldTotal;j<total;j++){
426 int otuI = otuIndex[childIndex];
427 int seqI = childSeqIndex[childIndex];
429 singleTau[j] = childSingleTau[childIndex];
431 aaP[otuI][nSeqsPerOTU[otuI]] = j;
432 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
437 int index = seqStart * numOTUs;
438 for(int j=seqStart;j<seqStop;j++){
439 for(int k=0;k<numOTUs;k++){
440 dist[index] = childDist[index];
448 cout << iter << '\t' << maxDelta << '\t' << setprecision(2) << nLL << '\t' << time(NULL) - cycTime << '\t' << setprecision(6) << (clock() - cycClock)/(double)CLOCKS_PER_SEC << endl;
450 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
452 for(int i=1;i<ncpus;i++){
453 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
458 for(int i=1;i<ncpus;i++){
459 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
465 cout << "\nFinalizing..." << endl;
468 vector<int> otuCounts(numOTUs, 0);
469 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
470 calcCentroidsDriver(0, numOTUs);
471 writeQualities(otuCounts);
472 writeSequences(otuCounts);
473 writeNames(otuCounts);
474 writeClusters(otuCounts);
477 remove(distFileName.c_str());
478 remove(namesFileName.c_str());
479 remove(listFileName.c_str());
481 cout << "Total time to process " << flowFileName << ":\t" << time(NULL) - begTime << '\t' << setprecision(6) << (clock() - begClock)/(double)CLOCKS_PER_SEC << endl;
489 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
490 if(abort){ return 0; }
493 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
495 for(int i=0;i<numFiles;i++){
496 //Now into the pyrodist part
498 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
499 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
500 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
501 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
503 flowDataIntI.resize(numSeqs * numFlowCells);
504 flowDataPrI.resize(numSeqs * numFlowCells);
505 mapUniqueToSeq.resize(numSeqs);
506 mapSeqToUnique.resize(numSeqs);
507 lengths.resize(numSeqs);
508 jointLookUp.resize(NUMBINS * NUMBINS);
510 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
511 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
512 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
513 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
514 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
515 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
516 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
518 flowFileName = string(fileName);
519 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
520 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
522 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
525 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
527 //Now into the pyronoise part
528 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
530 singleLookUp.resize(HOMOPS * NUMBINS);
531 uniqueFlowgrams.resize(numUniques * numFlowCells);
532 weight.resize(numOTUs);
533 centroids.resize(numOTUs);
534 change.resize(numOTUs);
535 dist.assign(numOTUs * numSeqs, 0);
536 nSeqsPerOTU.resize(numOTUs);
537 cumNumSeqs.resize(numOTUs);
539 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
540 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
541 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
543 int startOTU = pid * numOTUs / ncpus;
544 int endOTU = (pid + 1) * numOTUs / ncpus;
546 int startSeq = pid * numSeqs / ncpus;
547 int endSeq = (pid + 1) * numSeqs /ncpus;
553 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
554 singleTau.assign(total, 0.0000);
555 seqNumber.assign(total, 0);
556 seqIndex.assign(total, 0);
558 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
559 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
560 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
561 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
562 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
563 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
564 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
566 calcCentroidsDriver(startOTU, endOTU);
568 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
569 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
572 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
573 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
574 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
576 vector<int> otuIndex(total, 0);
577 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
578 total = otuIndex.size();
580 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
581 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
582 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
583 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
584 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
586 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
591 MPI_Barrier(MPI_COMM_WORLD);
595 catch(exception& e) {
596 m->errorOut(e, "ShhherCommand", "execute");
601 /**************************************************************************************************/
603 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
605 ostringstream outStream;
606 outStream.setf(ios::fixed, ios::floatfield);
607 outStream.setf(ios::dec, ios::basefield);
608 outStream.setf(ios::showpoint);
609 outStream.precision(6);
611 int begTime = time(NULL);
612 double begClock = clock();
614 for(int i=startSeq;i<stopSeq;i++){
615 for(int j=0;j<i;j++){
616 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
618 if(flowDistance < 1e-6){
619 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
621 else if(flowDistance <= cutoff){
622 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
626 cout << i << "\t" << (time(NULL) - begTime) << "\t" << (clock()-begClock)/CLOCKS_PER_SEC << endl;
629 cout << stopSeq << "\t" << (time(NULL) - begTime) << "\t" << (clock()-begClock)/CLOCKS_PER_SEC << endl;
631 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.dist";
632 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
634 ofstream distFile(fDistFileName.c_str());
635 distFile << outStream.str();
638 return fDistFileName;
640 catch(exception& e) {
641 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
647 //**********************************************************************************************************************
649 int ShhherCommand::execute(){
651 if (abort == true) { return 0; }
653 cout.setf(ios::fixed, ios::floatfield);
654 cout.setf(ios::showpoint);
660 vector<string> flowFileVector;
661 if(flowFilesFileName != "not found"){
664 ifstream flowFilesFile;
665 m->openInputFile(flowFilesFileName, flowFilesFile);
666 while(flowFilesFile){
667 flowFilesFile >> fName;
668 flowFileVector.push_back(fName);
669 m->gobble(flowFilesFile);
673 flowFileVector.push_back(flowFileName);
675 int numFiles = flowFileVector.size();
678 for(int i=0;i<numFiles;i++){
679 flowFileName = flowFileVector[i];
681 cout << "\n>>>>>\tProcessing " << flowFileName << " (file " << i+1 << " of " << numFiles << ")\t<<<<<" << endl;
682 cout << "Reading flowgrams..." << endl;
684 cout << "Identifying unique flowgrams..." << endl;
688 cout << "Calculating distances between flowgrams..." << endl;
689 string distFileName = createDistFile(processors);
690 string namesFileName = createNamesFile();
692 cout << "\nClustering flowgrams..." << endl;
693 string listFileName = cluster(distFileName, namesFileName);
694 getOTUData(listFileName);
701 double begClock = clock();
702 unsigned long int begTime = time(NULL);
704 cout << "\nDenoising flowgrams..." << endl;
705 cout << "iter\tmaxDelta\tnLL\t\tcycletime" << endl;
707 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
709 double cycClock = clock();
710 unsigned long int cycTime = time(NULL);
715 maxDelta = getNewWeights();
716 double nLL = getLikelihood();
723 cout << iter << '\t' << maxDelta << '\t' << setprecision(2) << nLL << '\t' << time(NULL) - cycTime << '\t' << setprecision(6) << (clock() - cycClock)/(double)CLOCKS_PER_SEC << endl;
726 cout << "\nFinalizing..." << endl;
730 vector<int> otuCounts(numOTUs, 0);
731 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
733 calcCentroidsDriver(0, numOTUs);
734 writeQualities(otuCounts);
735 writeSequences(otuCounts);
736 writeNames(otuCounts);
737 writeClusters(otuCounts);
740 remove(distFileName.c_str());
741 remove(namesFileName.c_str());
742 remove(listFileName.c_str());
744 cout << "Total time to process " << flowFileName << ":\t" << time(NULL) - begTime << '\t' << setprecision(6) << (clock() - begClock)/(double)CLOCKS_PER_SEC << endl;
748 catch(exception& e) {
749 m->errorOut(e, "ShhherCommand", "execute");
754 /**************************************************************************************************/
756 void ShhherCommand::getFlowData(){
759 m->openInputFile(flowFileName, flowFile);
763 int currentNumFlowCells;
767 flowFile >> numFlowCells;
768 int index = 0;//pcluster
769 while(!flowFile.eof()){
770 flowFile >> seqName >> currentNumFlowCells;
771 lengths.push_back(currentNumFlowCells);
773 seqNameVector.push_back(seqName);
774 nameMap[seqName] = index++;//pcluster
776 for(int i=0;i<numFlowCells;i++){
777 flowFile >> intensity;
778 if(intensity > 9.99) { intensity = 9.99; }
779 int intI = int(100 * intensity + 0.0001);
780 flowDataIntI.push_back(intI);
786 numSeqs = seqNameVector.size();
788 for(int i=0;i<numSeqs;i++){
789 int iNumFlowCells = i * numFlowCells;
790 for(int j=lengths[i];j<numFlowCells;j++){
791 flowDataIntI[iNumFlowCells + j] = 0;
796 catch(exception& e) {
797 m->errorOut(e, "ShhherCommand", "getFlowData");
802 /**************************************************************************************************/
804 void ShhherCommand::getSingleLookUp(){
806 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
807 singleLookUp.assign(HOMOPS * NUMBINS, 0);
811 m->openInputFile(lookupFileName, lookUpFile);
813 for(int i=0;i<HOMOPS;i++){
815 lookUpFile >> logFracFreq;
817 for(int j=0;j<NUMBINS;j++) {
818 lookUpFile >> singleLookUp[index];
824 catch(exception& e) {
825 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
830 /**************************************************************************************************/
832 void ShhherCommand::getJointLookUp(){
835 // the most likely joint probability (-log) that two intenities have the same polymer length
836 jointLookUp.resize(NUMBINS * NUMBINS, 0);
838 for(int i=0;i<NUMBINS;i++){
839 for(int j=0;j<NUMBINS;j++){
841 float minSum = 10000000000;
843 for(int k=0;k<HOMOPS;k++){
844 float sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
846 if(sum < minSum) { minSum = sum; }
848 jointLookUp[i * NUMBINS + j] = minSum;
852 catch(exception& e) {
853 m->errorOut(e, "ShhherCommand", "getJointLookUp");
858 /**************************************************************************************************/
860 float ShhherCommand::getProbIntensity(int intIntensity){
862 float minNegLogProb = 10000000000;
864 for(int i=0;i<HOMOPS;i++){//loop signal strength
865 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
866 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
869 return minNegLogProb;
871 catch(exception& e) {
872 m->errorOut(e, "ShhherCommand", "getProbIntensity");
877 /**************************************************************************************************/
879 void ShhherCommand::getUniques(){
884 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
885 uniqueCount.assign(numSeqs, 0); // anWeights
886 uniqueLengths.assign(numSeqs, 0);
887 mapSeqToUnique.assign(numSeqs, -1);
888 mapUniqueToSeq.assign(numSeqs, -1);
890 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
892 for(int i=0;i<numSeqs;i++){
895 vector<short> current(numFlowCells);
896 for(int j=0;j<numFlowCells;j++){ current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0)); }
898 for(int j=0;j<numUniques;j++){
899 int offset = j * numFlowCells;
902 for(int k=0;k<numFlowCells;k++){
903 if(current[k] != uniqueFlowgrams[offset + k]){
910 mapSeqToUnique[i] = j;
918 if(index == numUniques){
919 uniqueLengths[numUniques] = lengths[i];
920 uniqueCount[numUniques] = 1;
921 mapSeqToUnique[i] = numUniques;//anMap
922 mapUniqueToSeq[numUniques] = i;//anF
924 for(int k=0;k<numFlowCells;k++){
925 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
926 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
932 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
933 uniqueLengths.resize(numUniques);
935 flowDataPrI.assign(numSeqs * numFlowCells, 0);
936 for(int i=0;i<flowDataPrI.size();i++) { flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
938 catch(exception& e) {
939 m->errorOut(e, "ShhherCommand", "getUniques");
944 /**************************************************************************************************/
946 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
948 int minLength = lengths[mapSeqToUnique[seqA]];
949 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
951 int ANumFlowCells = seqA * numFlowCells;
952 int BNumFlowCells = seqB * numFlowCells;
956 for(int i=0;i<minLength;i++){
957 int flowAIntI = flowDataIntI[ANumFlowCells + i];
958 float flowAPrI = flowDataPrI[ANumFlowCells + i];
960 int flowBIntI = flowDataIntI[BNumFlowCells + i];
961 float flowBPrI = flowDataPrI[BNumFlowCells + i];
962 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
965 dist /= (float) minLength;
968 catch(exception& e) {
969 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
974 /**************************************************************************************************/
976 void ShhherCommand::flowDistParentFork(string distFileName, int startSeq, int stopSeq){
979 ostringstream outStream;
980 outStream.setf(ios::fixed, ios::floatfield);
981 outStream.setf(ios::dec, ios::basefield);
982 outStream.setf(ios::showpoint);
983 outStream.precision(6);
985 int begTime = time(NULL);
986 double begClock = clock();
988 for(int i=startSeq;i<stopSeq;i++){
989 for(int j=0;j<i;j++){
990 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
992 if(flowDistance < 1e-6){
993 outStream << seqNameVector[mapUniqueToSeq[i]] << '\t' << seqNameVector[mapUniqueToSeq[j]] << '\t' << 0.000000 << endl;
995 else if(flowDistance <= cutoff){
996 outStream << seqNameVector[mapUniqueToSeq[i]] << '\t' << seqNameVector[mapUniqueToSeq[j]] << '\t' << flowDistance << endl;
1000 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
1001 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1002 m->mothurOutEndLine();
1005 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
1006 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1007 m->mothurOutEndLine();
1009 ofstream distFile(distFileName.c_str());
1010 distFile << outStream.str();
1013 catch(exception& e) {
1014 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
1019 /**************************************************************************************************/
1021 string ShhherCommand::createDistFile(int processors){
1023 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.dist";
1025 unsigned long int begTime = time(NULL);
1026 double begClock = clock();
1031 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1032 if(processors == 1) { flowDistParentFork(fDistFileName, 0, numUniques); }
1033 else{ //you have multiple processors
1035 if (numSeqs < processors){ processors = 1; }
1037 vector<int> start(processors, 0);
1038 vector<int> end(processors, 0);
1040 for (int i = 0; i < processors; i++) {
1041 start[i] = int(sqrt(float(i)/float(processors)) * numUniques);
1042 end[i] = int(sqrt(float(i+1)/float(processors)) * numUniques);
1046 vector<int> processIDs;
1048 //loop through and create all the processes you want
1049 while (process != processors) {
1053 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1055 }else if (pid == 0){
1056 flowDistParentFork(fDistFileName + toString(getpid()) + ".temp", start[process], end[process]);
1059 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1061 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1066 //parent does its part
1067 flowDistParentFork(fDistFileName, start[0], end[0]);
1069 //force parent to wait until all the processes are done
1070 for (int i=0;i<processIDs.size();i++) {
1071 int temp = processIDs[i];
1075 //append and remove temp files
1076 for (int i=0;i<processIDs.size();i++) {
1077 m->appendFiles((fDistFileName + toString(processIDs[i]) + ".temp"), fDistFileName);
1078 remove((fDistFileName + toString(processIDs[i]) + ".temp").c_str());
1084 flowDistParentFork(fDistFileName, 0, numUniques);
1087 m->mothurOutEndLine();
1089 cout << "Total time: " << (time(NULL) - begTime) << "\t" << (clock() - begClock)/CLOCKS_PER_SEC << endl;;
1091 return fDistFileName;
1093 catch(exception& e) {
1094 m->errorOut(e, "ShhherCommand", "createDistFile");
1100 /**************************************************************************************************/
1102 string ShhherCommand::createNamesFile(){
1105 vector<string> duplicateNames(numUniques, "");
1106 for(int i=0;i<numSeqs;i++){
1107 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
1110 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.names";
1113 m->openOutputFile(nameFileName, nameFile);
1115 for(int i=0;i<numUniques;i++){
1116 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1117 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1121 return nameFileName;
1123 catch(exception& e) {
1124 m->errorOut(e, "ShhherCommand", "createNamesFile");
1129 //**********************************************************************************************************************
1131 string ShhherCommand::cluster(string distFileName, string namesFileName){
1134 SparseMatrix* matrix;
1136 RAbundVector* rabund;
1138 globaldata->setNameFile(namesFileName);
1139 globaldata->setColumnFile(distFileName);
1140 globaldata->setFormat("column");
1142 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1143 read->setCutoff(cutoff);
1145 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1146 clusterNameMap->readMap();
1147 read->read(clusterNameMap);
1149 list = read->getListVector();
1150 matrix = read->getMatrix();
1153 delete clusterNameMap;
1155 rabund = new RAbundVector(list->getRAbundVector());
1157 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1158 string tag = cluster->getTag();
1160 double clusterCutoff = cutoff;
1161 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1162 cluster->update(clusterCutoff);
1163 float dist = matrix->getSmallDist();
1166 list->setLabel(toString(cutoff));
1168 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.list";
1170 m->openOutputFile(listFileName, listFile);
1171 list->print(listFile);
1174 delete matrix; delete cluster; delete rabund; delete list;
1176 return listFileName;
1178 catch(exception& e) {
1179 m->errorOut(e, "ShhherCommand", "cluster");
1184 /**************************************************************************************************/
1186 void ShhherCommand::getOTUData(string listFileName){
1190 m->openInputFile(listFileName, listFile);
1193 listFile >> label >> numOTUs;
1195 otuData.assign(numSeqs, 0);
1196 cumNumSeqs.assign(numOTUs, 0);
1197 nSeqsPerOTU.assign(numOTUs, 0);
1198 aaP.resize(numOTUs);
1200 string singleOTU = "";
1202 for(int i=0;i<numOTUs;i++){
1204 listFile >> singleOTU;
1206 istringstream otuString(singleOTU);
1210 string seqName = "";
1212 for(int j=0;j<singleOTU.length();j++){
1213 char letter = otuString.get();
1219 map<string,int>::iterator nmIt = nameMap.find(seqName);
1220 int index = nmIt->second;
1222 nameMap.erase(nmIt);
1226 aaP[i].push_back(index);
1231 map<string,int>::iterator nmIt = nameMap.find(seqName);
1233 int index = nmIt->second;
1234 nameMap.erase(nmIt);
1238 aaP[i].push_back(index);
1243 sort(aaP[i].begin(), aaP[i].end());
1244 for(int j=0;j<nSeqsPerOTU[i];j++){
1245 seqNumber.push_back(aaP[i][j]);
1247 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
1248 aaP[i].push_back(0);
1252 for(int i=1;i<numOTUs;i++){
1253 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
1256 seqIndex = seqNumber;
1260 catch(exception& e) {
1261 m->errorOut(e, "ShhherCommand", "getOTUData");
1266 /**************************************************************************************************/
1268 void ShhherCommand::initPyroCluster(){
1270 dist.assign(numSeqs * numOTUs, 0);
1271 change.assign(numOTUs, 1);
1272 centroids.assign(numOTUs, -1);
1273 weight.assign(numOTUs, 0);
1274 singleTau.assign(numSeqs, 1.0);
1276 nSeqsBreaks.assign(processors+1, 0);
1277 nOTUsBreaks.assign(processors+1, 0);
1280 for(int i=0;i<processors;i++){
1281 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
1282 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
1284 nSeqsBreaks[processors] = numSeqs;
1285 nOTUsBreaks[processors] = numOTUs;
1287 catch(exception& e) {
1288 m->errorOut(e, "ShhherCommand", "initPyroCluster");
1293 /**************************************************************************************************/
1295 void ShhherCommand::fill(){
1298 for(int i=0;i<numOTUs;i++){
1299 cumNumSeqs[i] = index;
1300 for(int j=0;j<nSeqsPerOTU[i];j++){
1301 seqNumber[index] = aaP[i][j];
1302 seqIndex[index] = aaI[i][j];
1308 catch(exception& e) {
1309 m->errorOut(e, "ShhherCommand", "fill");
1314 /**************************************************************************************************/
1316 void ShhherCommand::calcCentroids(){
1319 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1321 if(processors == 1) {
1322 calcCentroidsDriver(0, numOTUs);
1324 else{ //you have multiple processors
1325 if (numOTUs < processors){ processors = 1; }
1328 vector<int> processIDs;
1330 //loop through and create all the processes you want
1331 while (process != processors) {
1335 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1337 }else if (pid == 0){
1338 calcCentroidsDriver(nOTUsBreaks[process], nOTUsBreaks[process+1]);
1341 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1343 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1348 //parent does its part
1349 calcCentroidsDriver(nOTUsBreaks[0], nOTUsBreaks[1]);
1351 //force parent to wait until all the processes are done
1352 for (int i=0;i<processIDs.size();i++) {
1353 int temp = processIDs[i];
1359 calcCentroidsDriver(0, numOTUs);
1362 catch(exception& e) {
1363 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1368 /**************************************************************************************************/
1370 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1372 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1377 for(int i=start;i<finish;i++){
1381 int minFlowGram = 100000000;
1382 double minFlowValue = 1e8;
1383 change[i] = 0; //FALSE
1385 for(int j=0;j<nSeqsPerOTU[i];j++){
1386 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1389 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1390 vector<double> adF(nSeqsPerOTU[i]);
1391 vector<int> anL(nSeqsPerOTU[i]);
1393 for(int j=0;j<nSeqsPerOTU[i];j++){
1394 int index = cumNumSeqs[i] + j;
1395 int nI = seqIndex[index];
1396 int nIU = mapSeqToUnique[nI];
1399 for(k=0;k<position;k++){
1405 anL[position] = nIU;
1406 adF[position] = 0.0000;
1411 for(int j=0;j<nSeqsPerOTU[i];j++){
1412 int index = cumNumSeqs[i] + j;
1413 int nI = seqIndex[index];
1414 int nIU = mapSeqToUnique[nI];
1416 double tauValue = singleTau[seqNumber[index]];
1418 for(int k=0;k<position;k++){
1419 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1420 adF[k] += dist * tauValue;
1424 for(int j=0;j<position;j++){
1425 if(adF[j] < minFlowValue){
1427 minFlowValue = adF[j];
1431 if(centroids[i] != anL[minFlowGram]){
1433 centroids[i] = anL[minFlowGram];
1436 else if(centroids[i] != -1){
1442 catch(exception& e) {
1443 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1448 /**************************************************************************************************/
1450 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1453 int flowAValue = cent * numFlowCells;
1454 int flowBValue = flow * numFlowCells;
1458 for(int i=0;i<length;i++){
1459 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1464 return dist / (double)length;
1466 catch(exception& e) {
1467 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1472 /**************************************************************************************************/
1474 double ShhherCommand::getNewWeights(){
1477 double maxChange = 0;
1479 for(int i=0;i<numOTUs;i++){
1481 double difference = weight[i];
1484 for(int j=0;j<nSeqsPerOTU[i];j++){
1485 int index = cumNumSeqs[i] + j;
1486 int nI = seqIndex[index];
1487 double tauValue = singleTau[seqNumber[index]];
1488 weight[i] += tauValue;
1491 difference = fabs(weight[i] - difference);
1492 if(difference > maxChange){ maxChange = difference; }
1496 catch(exception& e) {
1497 m->errorOut(e, "ShhherCommand", "getNewWeights");
1502 /**************************************************************************************************/
1504 double ShhherCommand::getLikelihood(){
1508 vector<double> P(numSeqs, 0);
1511 for(int i=0;i<numOTUs;i++){
1512 if(weight[i] > MIN_WEIGHT){
1517 for(int i=0;i<numOTUs;i++){
1518 for(int j=0;j<nSeqsPerOTU[i];j++){
1519 int index = cumNumSeqs[i] + j;
1520 int nI = seqIndex[index];
1521 double singleDist = dist[seqNumber[index]];
1523 P[nI] += weight[i] * exp(-singleDist * sigma);
1528 for(int i=0;i<numSeqs;i++){
1532 nLL = nLL -(double)numSeqs * log(sigma);
1536 catch(exception& e) {
1537 m->errorOut(e, "ShhherCommand", "getNewWeights");
1542 /**************************************************************************************************/
1544 void ShhherCommand::checkCentroids(){
1546 vector<int> unique(numOTUs, 1);
1548 for(int i=0;i<numOTUs;i++){
1549 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1554 for(int i=0;i<numOTUs;i++){
1556 for(int j=i+1;j<numOTUs;j++){
1559 if(centroids[j] == centroids[i]){
1563 weight[i] += weight[j];
1571 catch(exception& e) {
1572 m->errorOut(e, "ShhherCommand", "checkCentroids");
1577 /**************************************************************************************************/
1579 void ShhherCommand::calcNewDistances(){
1582 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1584 if(processors == 1) {
1585 calcNewDistancesParent(0, numSeqs);
1587 else{ //you have multiple processors
1588 if (numSeqs < processors){ processors = 1; }
1590 vector<vector<int> > child_otuIndex(processors);
1591 vector<vector<int> > child_seqIndex(processors);
1592 vector<vector<double> > child_singleTau(processors);
1593 vector<int> totals(processors);
1596 vector<int> processIDs;
1598 //loop through and create all the processes you want
1599 while (process != processors) {
1603 processIDs.push_back(pid); //create map from line number to pid so you can append files in correct order later
1605 }else if (pid == 0){
1606 calcNewDistancesChild(nSeqsBreaks[process], nSeqsBreaks[process+1], child_otuIndex[process], child_seqIndex[process], child_singleTau[process]);
1607 totals[process] = child_otuIndex[process].size();
1611 m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1613 for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill (temp, SIGINT); }
1618 //parent does its part
1619 calcNewDistancesParent(nSeqsBreaks[0], nSeqsBreaks[1]);
1620 int total = seqIndex.size();
1622 //force parent to wait until all the processes are done
1623 for (int i=0;i<processIDs.size();i++) {
1624 int temp = processIDs[i];
1628 for(int i=1;i<processors;i++){
1629 int oldTotal = total;
1632 singleTau.resize(total, 0);
1633 seqIndex.resize(total, 0);
1634 seqNumber.resize(total, 0);
1638 for(int j=oldTotal;j<total;j++){
1639 int otuI = child_otuIndex[i][childIndex];
1640 int seqI = child_seqIndex[i][childIndex];
1642 singleTau[j] = child_singleTau[i][childIndex];
1643 aaP[otuI][nSeqsPerOTU[otuI]] = j;
1644 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
1645 nSeqsPerOTU[otuI]++;
1652 calcNewDistancesParent(0, numSeqs);
1655 catch(exception& e) {
1656 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1661 /**************************************************************************************************/
1663 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1666 vector<double> newTau(numOTUs,0);
1667 vector<double> norms(numSeqs, 0);
1670 singleTau.resize(0);
1674 for(int i=startSeq;i<stopSeq;i++){
1675 double offset = 1e8;
1676 int indexOffset = i * numOTUs;
1678 for(int j=0;j<numOTUs;j++){
1680 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1681 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1683 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1684 offset = dist[indexOffset + j];
1688 for(int j=0;j<numOTUs;j++){
1689 if(weight[j] > MIN_WEIGHT){
1690 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1691 norms[i] += newTau[j];
1698 for(int j=0;j<numOTUs;j++){
1700 newTau[j] /= norms[i];
1702 if(newTau[j] > MIN_TAU){
1703 otuIndex.push_back(j);
1704 seqIndex.push_back(i);
1705 singleTau.push_back(newTau[j]);
1711 catch(exception& e) {
1712 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1717 /**************************************************************************************************/
1719 void ShhherCommand::calcNewDistancesChild(int startSeq, int stopSeq, vector<int>& child_otuIndex, vector<int>& child_seqIndex, vector<double>& child_singleTau){
1722 vector<double> newTau(numOTUs,0);
1723 vector<double> norms(numSeqs, 0);
1724 child_otuIndex.resize(0);
1725 child_seqIndex.resize(0);
1726 child_singleTau.resize(0);
1728 for(int i=startSeq;i<stopSeq;i++){
1729 double offset = 1e8;
1730 int indexOffset = i * numOTUs;
1733 for(int j=0;j<numOTUs;j++){
1734 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1735 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1738 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1739 offset = dist[indexOffset + j];
1743 for(int j=0;j<numOTUs;j++){
1744 if(weight[j] > MIN_WEIGHT){
1745 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1746 norms[i] += newTau[j];
1753 for(int j=0;j<numOTUs;j++){
1754 newTau[j] /= norms[i];
1756 if(newTau[j] > MIN_TAU){
1757 child_otuIndex.push_back(j);
1758 child_seqIndex.push_back(i);
1759 child_singleTau.push_back(newTau[j]);
1764 catch(exception& e) {
1765 m->errorOut(e, "ShhherCommand", "calcNewDistancesChild");
1770 /**************************************************************************************************/
1772 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1777 vector<double> newTau(numOTUs,0);
1778 vector<double> norms(numSeqs, 0);
1779 nSeqsPerOTU.assign(numOTUs, 0);
1781 for(int i=startSeq;i<stopSeq;i++){
1782 int indexOffset = i * numOTUs;
1784 double offset = 1e8;
1786 for(int j=0;j<numOTUs;j++){
1787 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1788 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1791 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1792 offset = dist[indexOffset + j];
1796 for(int j=0;j<numOTUs;j++){
1797 if(weight[j] > MIN_WEIGHT){
1798 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1799 norms[i] += newTau[j];
1806 for(int j=0;j<numOTUs;j++){
1807 newTau[j] /= norms[i];
1810 for(int j=0;j<numOTUs;j++){
1811 if(newTau[j] > MIN_TAU){
1813 int oldTotal = total;
1817 singleTau.resize(total, 0);
1818 seqNumber.resize(total, 0);
1819 seqIndex.resize(total, 0);
1821 singleTau[oldTotal] = newTau[j];
1823 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1824 aaI[j][nSeqsPerOTU[j]] = i;
1830 catch(exception& e) {
1831 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1836 /**************************************************************************************************/
1838 void ShhherCommand::setOTUs(){
1841 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1843 for(int i=0;i<numOTUs;i++){
1844 for(int j=0;j<nSeqsPerOTU[i];j++){
1845 int index = cumNumSeqs[i] + j;
1846 double tauValue = singleTau[seqNumber[index]];
1847 int sIndex = seqIndex[index];
1848 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1852 for(int i=0;i<numSeqs;i++){
1853 double maxTau = -1.0000;
1855 for(int j=0;j<numOTUs;j++){
1856 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1857 maxTau = bigTauMatrix[i * numOTUs + j];
1862 otuData[i] = maxOTU;
1865 nSeqsPerOTU.assign(numOTUs, 0);
1867 for(int i=0;i<numSeqs;i++){
1868 int index = otuData[i];
1870 singleTau[i] = 1.0000;
1873 aaP[index][nSeqsPerOTU[index]] = i;
1874 aaI[index][nSeqsPerOTU[index]] = i;
1876 nSeqsPerOTU[index]++;
1880 catch(exception& e) {
1881 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1886 /**************************************************************************************************/
1888 void ShhherCommand::writeQualities(vector<int> otuCounts){
1891 string qualityFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.qual";
1893 ofstream qualityFile;
1894 m->openOutputFile(qualityFileName, qualityFile);
1896 qualityFile.setf(ios::fixed, ios::floatfield);
1897 qualityFile.setf(ios::showpoint);
1898 qualityFile << setprecision(6);
1900 vector<vector<int> > qualities(numOTUs);
1901 vector<double> pr(HOMOPS, 0);
1905 for(int i=0;i<numOTUs;i++){
1909 if(nSeqsPerOTU[i] > 0){
1910 qualities[i].assign(1024, -1);
1912 while(index < numFlowCells){
1913 double maxPrValue = 1e8;
1914 short maxPrIndex = -1;
1915 double count = 0.0000;
1917 pr.assign(HOMOPS, 0);
1919 for(int j=0;j<nSeqsPerOTU[i];j++){
1920 int lIndex = cumNumSeqs[i] + j;
1921 double tauValue = singleTau[seqNumber[lIndex]];
1922 int sequenceIndex = aaI[i][j];
1923 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1927 for(int s=0;s<HOMOPS;s++){
1928 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1932 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1933 maxPrValue = pr[maxPrIndex];
1935 if(count > MIN_COUNT){
1937 double norm = 0.0000;
1939 for(int s=0;s<HOMOPS;s++){
1940 norm += exp(-(pr[s] - maxPrValue));
1943 for(int s=1;s<=maxPrIndex;s++){
1945 double temp = 0.0000;
1947 U += exp(-(pr[s-1]-maxPrValue))/norm;
1955 temp = floor(-10 * temp);
1956 value = (int)floor(temp);
1957 if(value > 100){ value = 100; }
1959 qualities[i][base] = (int)value;
1969 if(otuCounts[i] > 0){
1970 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1972 int j=4; //need to get past the first four bases
1973 while(qualities[i][j] != -1){
1974 qualityFile << qualities[i][j] << ' ';
1977 qualityFile << endl;
1980 qualityFile.close();
1983 catch(exception& e) {
1984 m->errorOut(e, "ShhherCommand", "writeQualities");
1989 /**************************************************************************************************/
1991 void ShhherCommand::writeSequences(vector<int> otuCounts){
1994 string bases = "TACG";
1996 string fastaFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.fasta";
1998 m->openOutputFile(fastaFileName, fastaFile);
2000 vector<string> names(numOTUs, "");
2002 for(int i=0;i<numOTUs;i++){
2003 int index = centroids[i];
2005 if(otuCounts[i] > 0){
2006 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
2008 for(int j=8;j<numFlowCells;j++){
2010 char base = bases[j % 4];
2011 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
2020 catch(exception& e) {
2021 m->errorOut(e, "ShhherCommand", "writeSequences");
2026 /**************************************************************************************************/
2028 void ShhherCommand::writeNames(vector<int> otuCounts){
2030 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.final.names";
2032 m->openOutputFile(nameFileName, nameFile);
2034 for(int i=0;i<numOTUs;i++){
2035 if(otuCounts[i] > 0){
2036 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
2038 for(int j=1;j<nSeqsPerOTU[i];j++){
2039 nameFile << ',' << seqNameVector[aaI[i][j]];
2047 catch(exception& e) {
2048 m->errorOut(e, "ShhherCommand", "writeNames");
2053 /**************************************************************************************************/
2055 void ShhherCommand::writeGroups(){
2057 string fileRoot = flowFileName.substr(0,flowFileName.find_last_of('.'));
2058 string groupFileName = fileRoot + ".pn.groups";
2060 m->openOutputFile(groupFileName, groupFile);
2062 for(int i=0;i<numSeqs;i++){
2063 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
2067 catch(exception& e) {
2068 m->errorOut(e, "ShhherCommand", "writeGroups");
2073 /**************************************************************************************************/
2075 void ShhherCommand::writeClusters(vector<int> otuCounts){
2077 string otuCountsFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.counts";
2078 ofstream otuCountsFile;
2079 m->openOutputFile(otuCountsFileName, otuCountsFile);
2081 string bases = "TACG";
2083 for(int i=0;i<numOTUs;i++){
2084 //output the translated version of the centroid sequence for the otu
2085 if(otuCounts[i] > 0){
2086 int index = centroids[i];
2088 otuCountsFile << "ideal\t";
2089 for(int j=8;j<numFlowCells;j++){
2090 char base = bases[j % 4];
2091 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
2092 otuCountsFile << base;
2095 otuCountsFile << endl;
2097 for(int j=0;j<nSeqsPerOTU[i];j++){
2098 int sequence = aaI[i][j];
2099 otuCountsFile << seqNameVector[sequence] << '\t';
2101 for(int k=8;k<lengths[sequence];k++){
2102 char base = bases[k % 4];
2103 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
2105 for(int s=0;s<freq;s++){
2106 otuCountsFile << base;
2109 otuCountsFile << endl;
2111 otuCountsFile << endl;
2114 otuCountsFile.close();
2116 catch(exception& e) {
2117 m->errorOut(e, "ShhherCommand", "writeClusters");
2122 //**********************************************************************************************************************