5 * Created by Pat Schloss on 12/27/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "shhhercommand.h"
12 //**********************************************************************************************************************
13 vector<string> ShhherCommand::setParameters(){
15 CommandParameter pflow("flow", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pflow);
16 CommandParameter pfile("file", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pfile);
17 CommandParameter plookup("lookup", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(plookup);
18 CommandParameter pcutoff("cutoff", "Number", "", "0.01", "", "", "",false,false); parameters.push_back(pcutoff);
19 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
20 CommandParameter pmaxiter("maxiter", "Number", "", "1000", "", "", "",false,false); parameters.push_back(pmaxiter);
21 CommandParameter psigma("sigma", "Number", "", "60", "", "", "",false,false); parameters.push_back(psigma);
22 CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "",false,false); parameters.push_back(pmindelta);
23 CommandParameter porder("order", "String", "", "", "", "", "",false,false); parameters.push_back(porder);
24 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
25 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
27 vector<string> myArray;
28 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
32 m->errorOut(e, "ShhherCommand", "setParameters");
36 //**********************************************************************************************************************
37 string ShhherCommand::getHelpString(){
39 string helpString = "";
40 helpString += "The shhh.flows command reads a file containing flowgrams and creates a file of corrected sequences.\n";
44 m->errorOut(e, "ShhherCommand", "getHelpString");
48 //**********************************************************************************************************************
50 ShhherCommand::ShhherCommand(){
52 abort = true; calledHelp = true;
55 //initialize outputTypes
56 // vector<string> tempOutNames;
57 // outputTypes["pn.dist"] = tempOutNames;
61 m->errorOut(e, "ShhherCommand", "ShhherCommand");
66 //**********************************************************************************************************************
68 ShhherCommand::ShhherCommand(string option) {
72 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
73 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
77 abort = false; calledHelp = false;
79 //allow user to run help
80 if(option == "help") { help(); abort = true; calledHelp = true; }
81 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
84 vector<string> myArray = setParameters();
86 OptionParser parser(option);
87 map<string,string> parameters = parser.getParameters();
89 ValidParameters validParameter;
90 map<string,string>::iterator it;
92 //check to make sure all parameters are valid for command
93 for (it = parameters.begin(); it != parameters.end(); it++) {
94 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
97 //initialize outputTypes
98 vector<string> tempOutNames;
99 // outputTypes["pn.dist"] = tempOutNames;
100 // outputTypes["fasta"] = tempOutNames;
102 //if the user changes the input directory command factory will send this info to us in the output parameter
103 string inputDir = validParameter.validFile(parameters, "inputdir", false);
104 if (inputDir == "not found"){ inputDir = ""; }
107 it = parameters.find("flow");
108 //user has given a template file
109 if(it != parameters.end()){
110 path = m->hasPath(it->second);
111 //if the user has not given a path then, add inputdir. else leave path alone.
112 if (path == "") { parameters["flow"] = inputDir + it->second; }
115 it = parameters.find("lookup");
116 //user has given a template file
117 if(it != parameters.end()){
118 path = m->hasPath(it->second);
119 //if the user has not given a path then, add inputdir. else leave path alone.
120 if (path == "") { parameters["lookup"] = inputDir + it->second; }
123 it = parameters.find("file");
124 //user has given a template file
125 if(it != parameters.end()){
126 path = m->hasPath(it->second);
127 //if the user has not given a path then, add inputdir. else leave path alone.
128 if (path == "") { parameters["file"] = inputDir + it->second; }
132 //if the user changes the output directory command factory will send this info to us in the output parameter
133 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
135 //check for required parameters
136 flowFileName = validParameter.validFile(parameters, "flow", true);
137 flowFilesFileName = validParameter.validFile(parameters, "file", true);
138 if (flowFileName == "not found" && flowFilesFileName == "not found") {
139 m->mothurOut("values for either flow or file must be provided for the shhh.seqs command.");
140 m->mothurOutEndLine();
143 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
145 if(flowFileName != "not found"){
146 compositeFASTAFileName = "";
147 compositeNamesFileName = "";
152 string thisoutputDir = m->hasPath(flowFilesFileName); //if user entered a file with a path then preserve it
154 //flow.files = 9 character offset
155 compositeFASTAFileName = thisoutputDir + flowFilesFileName.substr(0, flowFilesFileName.length()-10) + "shhh.fasta";
156 m->openOutputFile(compositeFASTAFileName, temp);
159 compositeNamesFileName = thisoutputDir + flowFilesFileName.substr(0, flowFilesFileName.length()-10) + "shhh.names";
160 m->openOutputFile(compositeNamesFileName, temp);
164 if(flowFilesFileName != "not found"){
167 ifstream flowFilesFile;
168 m->openInputFile(flowFilesFileName, flowFilesFile);
169 while(flowFilesFile){
170 fName = m->getline(flowFilesFile);
172 //test if file is valid
174 int ableToOpen = m->openInputFile(fName, in, "noerror");
176 if (ableToOpen == 1) {
177 if (inputDir != "") { //default path is set
178 string tryPath = inputDir + fName;
179 m->mothurOut("Unable to open " + fName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
181 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
187 if (ableToOpen == 1) {
188 if (m->getDefaultPath() != "") { //default path is set
189 string tryPath = m->getDefaultPath() + m->getSimpleName(fName);
190 m->mothurOut("Unable to open " + fName + ". Trying default " + tryPath); m->mothurOutEndLine();
192 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
198 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
199 if (ableToOpen == 1) {
200 string exepath = m->argv;
201 string tempPath = exepath;
202 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
203 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
205 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(fName);
206 m->mothurOut("Unable to open " + fName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
208 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
213 if (ableToOpen == 1) { m->mothurOut("Unable to open " + fName + ". Disregarding. "); m->mothurOutEndLine(); }
214 else { flowFileVector.push_back(fName); }
215 m->gobble(flowFilesFile);
217 flowFilesFile.close();
218 if (flowFileVector.size() == 0) { m->mothurOut("[ERROR]: no valid files."); m->mothurOutEndLine(); abort = true; }
221 outputDir += m->hasPath(flowFileName);
222 flowFileVector.push_back(flowFileName);
225 //check for optional parameter and set defaults
226 // ...at some point should added some additional type checking...
228 temp = validParameter.validFile(parameters, "lookup", true);
229 if (temp == "not found") {
230 lookupFileName = "LookUp_Titanium.pat";
234 ableToOpen = m->openInputFile(lookupFileName, in, "noerror");
237 //if you can't open it, try input location
238 if (ableToOpen == 1) {
239 if (inputDir != "") { //default path is set
240 string tryPath = inputDir + lookupFileName;
241 m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
243 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
245 lookupFileName = tryPath;
249 //if you can't open it, try default location
250 if (ableToOpen == 1) {
251 if (m->getDefaultPath() != "") { //default path is set
252 string tryPath = m->getDefaultPath() + m->getSimpleName(lookupFileName);
253 m->mothurOut("Unable to open " + lookupFileName + ". Trying default " + tryPath); m->mothurOutEndLine();
255 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
257 lookupFileName = tryPath;
261 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
262 if (ableToOpen == 1) {
263 string exepath = m->argv;
264 string tempPath = exepath;
265 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
266 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
268 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
269 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
271 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
273 lookupFileName = tryPath;
276 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
278 else if(temp == "not open") {
280 lookupFileName = validParameter.validFile(parameters, "lookup", false);
282 //if you can't open it its not inputDir, try mothur excutable location
283 string exepath = m->argv;
284 string tempPath = exepath;
285 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
286 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
288 string tryPath = m->getFullPathName(exepath) + lookupFileName;
289 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
291 int ableToOpen = m->openInputFile(tryPath, in2, "noerror");
293 lookupFileName = tryPath;
295 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
296 }else { lookupFileName = temp; }
298 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
299 m->setProcessors(temp);
300 m->mothurConvert(temp, processors);
302 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
303 m->mothurConvert(temp, cutoff);
305 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
306 m->mothurConvert(temp, minDelta);
308 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
309 m->mothurConvert(temp, maxIters);
311 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
312 m->mothurConvert(temp, sigma);
314 flowOrder = validParameter.validFile(parameters, "order", false);
315 if (flowOrder == "not found"){ flowOrder = "TACG"; }
316 else if(flowOrder.length() != 4){
317 m->mothurOut("The value of the order option must be four bases long\n");
325 catch(exception& e) {
326 m->errorOut(e, "ShhherCommand", "ShhherCommand");
330 //**********************************************************************************************************************
332 int ShhherCommand::execute(){
334 if (abort == true) { if (calledHelp) { return 0; } return 2; }
341 for(int i=1;i<ncpus;i++){
342 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
344 if(abort == 1){ return 0; }
348 m->mothurOut("\nGetting preliminary data...\n");
349 getSingleLookUp(); if (m->control_pressed) { return 0; }
350 getJointLookUp(); if (m->control_pressed) { return 0; }
352 vector<string> flowFileVector;
353 if(flowFilesFileName != "not found"){
356 ifstream flowFilesFile;
357 m->openInputFile(flowFilesFileName, flowFilesFile);
358 while(flowFilesFile){
359 fName = m->getline(flowFilesFile);
360 flowFileVector.push_back(fName);
361 m->gobble(flowFilesFile);
365 flowFileVector.push_back(flowFileName);
368 int numFiles = flowFileVector.size();
370 for(int i=1;i<ncpus;i++){
371 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
374 for(int i=0;i<numFiles;i++){
376 if (m->control_pressed) { break; }
378 double begClock = clock();
379 unsigned long long begTime = time(NULL);
381 flowFileName = flowFileVector[i];
383 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
384 m->mothurOut("Reading flowgrams...\n");
387 if (m->control_pressed) { break; }
389 m->mothurOut("Identifying unique flowgrams...\n");
392 if (m->control_pressed) { break; }
394 m->mothurOut("Calculating distances between flowgrams...\n");
396 strcpy(fileName, flowFileName.c_str());
398 for(int i=1;i<ncpus;i++){
399 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
401 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
402 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
403 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
404 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
405 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
406 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
407 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
408 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
409 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
410 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
413 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
415 if (m->control_pressed) { break; }
418 for(int i=1;i<ncpus;i++){
419 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
421 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
422 m->mothurRemove((distFileName + ".temp." + toString(i)));
425 string namesFileName = createNamesFile();
427 if (m->control_pressed) { break; }
429 m->mothurOut("\nClustering flowgrams...\n");
430 string listFileName = cluster(distFileName, namesFileName);
432 if (m->control_pressed) { break; }
434 getOTUData(listFileName);
436 m->mothurRemove(distFileName);
437 m->mothurRemove(namesFileName);
438 m->mothurRemove(listFileName);
440 if (m->control_pressed) { break; }
444 if (m->control_pressed) { break; }
446 for(int i=1;i<ncpus;i++){
447 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
448 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
449 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
450 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
453 if (m->control_pressed) { break; }
458 int numOTUsOnCPU = numOTUs / ncpus;
459 int numSeqsOnCPU = numSeqs / ncpus;
460 m->mothurOut("\nDenoising flowgrams...\n");
461 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
463 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
465 double cycClock = clock();
466 unsigned long long cycTime = time(NULL);
469 if (m->control_pressed) { break; }
471 int total = singleTau.size();
472 for(int i=1;i<ncpus;i++){
473 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
474 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
475 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
477 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
478 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
479 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
480 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
481 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
484 calcCentroidsDriver(0, numOTUsOnCPU);
486 for(int i=1;i<ncpus;i++){
487 int otuStart = i * numOTUs / ncpus;
488 int otuStop = (i + 1) * numOTUs / ncpus;
490 vector<int> tempCentroids(numOTUs, 0);
491 vector<short> tempChange(numOTUs, 0);
493 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
494 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
496 for(int j=otuStart;j<otuStop;j++){
497 centroids[j] = tempCentroids[j];
498 change[j] = tempChange[j];
502 maxDelta = getNewWeights(); if (m->control_pressed) { break; }
503 double nLL = getLikelihood(); if (m->control_pressed) { break; }
504 checkCentroids(); if (m->control_pressed) { break; }
506 for(int i=1;i<ncpus;i++){
507 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
508 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
509 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
512 calcNewDistancesParent(0, numSeqsOnCPU);
514 total = singleTau.size();
516 for(int i=1;i<ncpus;i++){
518 int seqStart = i * numSeqs / ncpus;
519 int seqStop = (i + 1) * numSeqs / ncpus;
521 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
523 vector<int> childSeqIndex(childTotal, 0);
524 vector<double> childSingleTau(childTotal, 0);
525 vector<double> childDist(numSeqs * numOTUs, 0);
526 vector<int> otuIndex(childTotal, 0);
528 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
529 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
530 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
531 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
533 int oldTotal = total;
535 singleTau.resize(total, 0);
536 seqIndex.resize(total, 0);
537 seqNumber.resize(total, 0);
541 for(int j=oldTotal;j<total;j++){
542 int otuI = otuIndex[childIndex];
543 int seqI = childSeqIndex[childIndex];
545 singleTau[j] = childSingleTau[childIndex];
547 aaP[otuI][nSeqsPerOTU[otuI]] = j;
548 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
553 int index = seqStart * numOTUs;
554 for(int j=seqStart;j<seqStop;j++){
555 for(int k=0;k<numOTUs;k++){
556 dist[index] = childDist[index];
564 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
566 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
568 for(int i=1;i<ncpus;i++){
569 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
574 for(int i=1;i<ncpus;i++){
575 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
581 if (m->control_pressed) { break; }
583 m->mothurOut("\nFinalizing...\n");
586 if (m->control_pressed) { break; }
590 vector<int> otuCounts(numOTUs, 0);
591 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
592 calcCentroidsDriver(0, numOTUs);
594 if (m->control_pressed) { break; }
596 writeQualities(otuCounts); if (m->control_pressed) { break; }
597 writeSequences(otuCounts); if (m->control_pressed) { break; }
598 writeNames(otuCounts); if (m->control_pressed) { break; }
599 writeClusters(otuCounts); if (m->control_pressed) { break; }
600 writeGroups(); if (m->control_pressed) { break; }
603 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
609 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
610 if(abort){ return 0; }
613 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
615 for(int i=0;i<numFiles;i++){
617 if (m->control_pressed) { break; }
619 //Now into the pyrodist part
623 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
624 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
625 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
626 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
628 flowDataIntI.resize(numSeqs * numFlowCells);
629 flowDataPrI.resize(numSeqs * numFlowCells);
630 mapUniqueToSeq.resize(numSeqs);
631 mapSeqToUnique.resize(numSeqs);
632 lengths.resize(numSeqs);
633 jointLookUp.resize(NUMBINS * NUMBINS);
635 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
636 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
637 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
638 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
639 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
640 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
641 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
643 flowFileName = string(fileName);
644 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
645 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
647 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
649 if (m->control_pressed) { break; }
652 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
654 //Now into the pyronoise part
655 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
657 singleLookUp.resize(HOMOPS * NUMBINS);
658 uniqueFlowgrams.resize(numUniques * numFlowCells);
659 weight.resize(numOTUs);
660 centroids.resize(numOTUs);
661 change.resize(numOTUs);
662 dist.assign(numOTUs * numSeqs, 0);
663 nSeqsPerOTU.resize(numOTUs);
664 cumNumSeqs.resize(numOTUs);
666 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
667 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
668 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
670 int startOTU = pid * numOTUs / ncpus;
671 int endOTU = (pid + 1) * numOTUs / ncpus;
673 int startSeq = pid * numSeqs / ncpus;
674 int endSeq = (pid + 1) * numSeqs /ncpus;
680 if (m->control_pressed) { break; }
682 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
683 singleTau.assign(total, 0.0000);
684 seqNumber.assign(total, 0);
685 seqIndex.assign(total, 0);
687 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
688 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
689 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
690 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
691 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
692 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
693 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
695 calcCentroidsDriver(startOTU, endOTU);
697 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
698 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
700 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
701 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
702 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
704 vector<int> otuIndex(total, 0);
705 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
706 total = otuIndex.size();
708 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
709 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
710 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
711 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
712 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
714 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
719 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
721 MPI_Barrier(MPI_COMM_WORLD);
724 if(compositeFASTAFileName != ""){
725 outputNames.push_back(compositeFASTAFileName);
726 outputNames.push_back(compositeNamesFileName);
729 m->mothurOutEndLine();
730 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
731 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
732 m->mothurOutEndLine();
737 catch(exception& e) {
738 m->errorOut(e, "ShhherCommand", "execute");
742 /**************************************************************************************************/
743 string ShhherCommand::createNamesFile(){
746 vector<string> duplicateNames(numUniques, "");
747 for(int i=0;i<numSeqs;i++){
748 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
751 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
754 m->openOutputFile(nameFileName, nameFile);
756 for(int i=0;i<numUniques;i++){
758 if (m->control_pressed) { break; }
760 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
761 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
767 catch(exception& e) {
768 m->errorOut(e, "ShhherCommand", "createNamesFile");
772 /**************************************************************************************************/
774 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
776 ostringstream outStream;
777 outStream.setf(ios::fixed, ios::floatfield);
778 outStream.setf(ios::dec, ios::basefield);
779 outStream.setf(ios::showpoint);
780 outStream.precision(6);
782 int begTime = time(NULL);
783 double begClock = clock();
785 for(int i=startSeq;i<stopSeq;i++){
787 if (m->control_pressed) { break; }
789 for(int j=0;j<i;j++){
790 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
792 if(flowDistance < 1e-6){
793 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
795 else if(flowDistance <= cutoff){
796 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
800 m->mothurOut(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
804 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
805 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
807 if (m->control_pressed) { return fDistFileName; }
809 m->mothurOut(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
811 ofstream distFile(fDistFileName.c_str());
812 distFile << outStream.str();
815 return fDistFileName;
817 catch(exception& e) {
818 m->errorOut(e, "ShhherCommand", "flowDistMPI");
822 /**************************************************************************************************/
824 void ShhherCommand::getOTUData(string listFileName){
828 m->openInputFile(listFileName, listFile);
831 listFile >> label >> numOTUs;
833 otuData.assign(numSeqs, 0);
834 cumNumSeqs.assign(numOTUs, 0);
835 nSeqsPerOTU.assign(numOTUs, 0);
836 aaP.clear();aaP.resize(numOTUs);
842 string singleOTU = "";
844 for(int i=0;i<numOTUs;i++){
846 if (m->control_pressed) { break; }
848 listFile >> singleOTU;
850 istringstream otuString(singleOTU);
856 for(int j=0;j<singleOTU.length();j++){
857 char letter = otuString.get();
863 map<string,int>::iterator nmIt = nameMap.find(seqName);
864 int index = nmIt->second;
870 aaP[i].push_back(index);
875 map<string,int>::iterator nmIt = nameMap.find(seqName);
877 int index = nmIt->second;
882 aaP[i].push_back(index);
887 sort(aaP[i].begin(), aaP[i].end());
888 for(int j=0;j<nSeqsPerOTU[i];j++){
889 seqNumber.push_back(aaP[i][j]);
891 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
898 for(int i=1;i<numOTUs;i++){
899 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
902 seqIndex = seqNumber;
907 catch(exception& e) {
908 m->errorOut(e, "ShhherCommand", "getOTUData");
913 /**************************************************************************************************/
915 void ShhherCommand::initPyroCluster(){
917 if (numOTUs < processors) { processors = 1; }
919 dist.assign(numSeqs * numOTUs, 0);
920 change.assign(numOTUs, 1);
921 centroids.assign(numOTUs, -1);
922 weight.assign(numOTUs, 0);
923 singleTau.assign(numSeqs, 1.0);
925 nSeqsBreaks.assign(processors+1, 0);
926 nOTUsBreaks.assign(processors+1, 0);
929 for(int i=0;i<processors;i++){
930 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
931 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
933 nSeqsBreaks[processors] = numSeqs;
934 nOTUsBreaks[processors] = numOTUs;
936 catch(exception& e) {
937 m->errorOut(e, "ShhherCommand", "initPyroCluster");
942 /**************************************************************************************************/
944 void ShhherCommand::fill(){
947 for(int i=0;i<numOTUs;i++){
949 if (m->control_pressed) { break; }
951 cumNumSeqs[i] = index;
952 for(int j=0;j<nSeqsPerOTU[i];j++){
953 seqNumber[index] = aaP[i][j];
954 seqIndex[index] = aaI[i][j];
960 catch(exception& e) {
961 m->errorOut(e, "ShhherCommand", "fill");
966 /**************************************************************************************************/
968 void ShhherCommand::getFlowData(){
971 m->openInputFile(flowFileName, flowFile);
974 seqNameVector.clear();
976 flowDataIntI.clear();
980 int currentNumFlowCells;
984 flowFile >> numFlowCells;
985 int index = 0;//pcluster
986 while(!flowFile.eof()){
988 if (m->control_pressed) { break; }
990 flowFile >> seqName >> currentNumFlowCells;
991 lengths.push_back(currentNumFlowCells);
993 seqNameVector.push_back(seqName);
994 nameMap[seqName] = index++;//pcluster
996 for(int i=0;i<numFlowCells;i++){
997 flowFile >> intensity;
998 if(intensity > 9.99) { intensity = 9.99; }
999 int intI = int(100 * intensity + 0.0001);
1000 flowDataIntI.push_back(intI);
1002 m->gobble(flowFile);
1006 numSeqs = seqNameVector.size();
1008 for(int i=0;i<numSeqs;i++){
1010 if (m->control_pressed) { break; }
1012 int iNumFlowCells = i * numFlowCells;
1013 for(int j=lengths[i];j<numFlowCells;j++){
1014 flowDataIntI[iNumFlowCells + j] = 0;
1019 catch(exception& e) {
1020 m->errorOut(e, "ShhherCommand", "getFlowData");
1024 /**************************************************************************************************/
1025 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1028 vector<double> newTau(numOTUs,0);
1029 vector<double> norms(numSeqs, 0);
1034 for(int i=startSeq;i<stopSeq;i++){
1036 if (m->control_pressed) { break; }
1038 double offset = 1e8;
1039 int indexOffset = i * numOTUs;
1041 for(int j=0;j<numOTUs;j++){
1043 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1044 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1046 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1047 offset = dist[indexOffset + j];
1051 for(int j=0;j<numOTUs;j++){
1052 if(weight[j] > MIN_WEIGHT){
1053 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1054 norms[i] += newTau[j];
1061 for(int j=0;j<numOTUs;j++){
1063 newTau[j] /= norms[i];
1065 if(newTau[j] > MIN_TAU){
1066 otuIndex.push_back(j);
1067 seqIndex.push_back(i);
1068 singleTau.push_back(newTau[j]);
1074 catch(exception& e) {
1075 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1080 /**************************************************************************************************/
1082 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1087 vector<double> newTau(numOTUs,0);
1088 vector<double> norms(numSeqs, 0);
1089 nSeqsPerOTU.assign(numOTUs, 0);
1091 for(int i=startSeq;i<stopSeq;i++){
1093 if (m->control_pressed) { break; }
1095 int indexOffset = i * numOTUs;
1097 double offset = 1e8;
1099 for(int j=0;j<numOTUs;j++){
1101 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1102 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1105 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1106 offset = dist[indexOffset + j];
1110 for(int j=0;j<numOTUs;j++){
1111 if(weight[j] > MIN_WEIGHT){
1112 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1113 norms[i] += newTau[j];
1120 for(int j=0;j<numOTUs;j++){
1121 newTau[j] /= norms[i];
1124 for(int j=0;j<numOTUs;j++){
1125 if(newTau[j] > MIN_TAU){
1127 int oldTotal = total;
1131 singleTau.resize(total, 0);
1132 seqNumber.resize(total, 0);
1133 seqIndex.resize(total, 0);
1135 singleTau[oldTotal] = newTau[j];
1137 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1138 aaI[j][nSeqsPerOTU[j]] = i;
1146 catch(exception& e) {
1147 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1152 /**************************************************************************************************/
1154 void ShhherCommand::setOTUs(){
1157 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1159 for(int i=0;i<numOTUs;i++){
1161 if (m->control_pressed) { break; }
1163 for(int j=0;j<nSeqsPerOTU[i];j++){
1164 int index = cumNumSeqs[i] + j;
1165 double tauValue = singleTau[seqNumber[index]];
1166 int sIndex = seqIndex[index];
1167 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1171 for(int i=0;i<numSeqs;i++){
1172 double maxTau = -1.0000;
1174 for(int j=0;j<numOTUs;j++){
1175 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1176 maxTau = bigTauMatrix[i * numOTUs + j];
1181 otuData[i] = maxOTU;
1184 nSeqsPerOTU.assign(numOTUs, 0);
1186 for(int i=0;i<numSeqs;i++){
1187 int index = otuData[i];
1189 singleTau[i] = 1.0000;
1192 aaP[index][nSeqsPerOTU[index]] = i;
1193 aaI[index][nSeqsPerOTU[index]] = i;
1195 nSeqsPerOTU[index]++;
1199 catch(exception& e) {
1200 m->errorOut(e, "ShhherCommand", "setOTUs");
1205 /**************************************************************************************************/
1207 void ShhherCommand::getUniques(){
1212 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
1213 uniqueCount.assign(numSeqs, 0); // anWeights
1214 uniqueLengths.assign(numSeqs, 0);
1215 mapSeqToUnique.assign(numSeqs, -1);
1216 mapUniqueToSeq.assign(numSeqs, -1);
1218 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
1220 for(int i=0;i<numSeqs;i++){
1222 if (m->control_pressed) { break; }
1226 vector<short> current(numFlowCells);
1227 for(int j=0;j<numFlowCells;j++){
1228 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
1231 for(int j=0;j<numUniques;j++){
1232 int offset = j * numFlowCells;
1236 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
1237 else { shorterLength = uniqueLengths[j]; }
1239 for(int k=0;k<shorterLength;k++){
1240 if(current[k] != uniqueFlowgrams[offset + k]){
1247 mapSeqToUnique[i] = j;
1250 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
1256 if(index == numUniques){
1257 uniqueLengths[numUniques] = lengths[i];
1258 uniqueCount[numUniques] = 1;
1259 mapSeqToUnique[i] = numUniques;//anMap
1260 mapUniqueToSeq[numUniques] = i;//anF
1262 for(int k=0;k<numFlowCells;k++){
1263 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
1264 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
1270 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
1271 uniqueLengths.resize(numUniques);
1273 flowDataPrI.resize(numSeqs * numFlowCells, 0);
1274 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
1276 catch(exception& e) {
1277 m->errorOut(e, "ShhherCommand", "getUniques");
1282 /**************************************************************************************************/
1284 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
1286 int minLength = lengths[mapSeqToUnique[seqA]];
1287 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
1289 int ANumFlowCells = seqA * numFlowCells;
1290 int BNumFlowCells = seqB * numFlowCells;
1294 for(int i=0;i<minLength;i++){
1296 if (m->control_pressed) { break; }
1298 int flowAIntI = flowDataIntI[ANumFlowCells + i];
1299 float flowAPrI = flowDataPrI[ANumFlowCells + i];
1301 int flowBIntI = flowDataIntI[BNumFlowCells + i];
1302 float flowBPrI = flowDataPrI[BNumFlowCells + i];
1303 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
1306 dist /= (float) minLength;
1309 catch(exception& e) {
1310 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
1315 //**********************************************************************************************************************/
1317 string ShhherCommand::cluster(string distFileName, string namesFileName){
1320 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1321 read->setCutoff(cutoff);
1323 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1324 clusterNameMap->readMap();
1325 read->read(clusterNameMap);
1327 ListVector* list = read->getListVector();
1328 SparseMatrix* matrix = read->getMatrix();
1331 delete clusterNameMap;
1333 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1335 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1336 string tag = cluster->getTag();
1338 double clusterCutoff = cutoff;
1339 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1341 if (m->control_pressed) { break; }
1343 cluster->update(clusterCutoff);
1346 list->setLabel(toString(cutoff));
1348 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
1350 m->openOutputFile(listFileName, listFile);
1351 list->print(listFile);
1354 delete matrix; delete cluster; delete rabund; delete list;
1356 return listFileName;
1358 catch(exception& e) {
1359 m->errorOut(e, "ShhherCommand", "cluster");
1364 /**************************************************************************************************/
1366 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1368 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1373 for(int i=start;i<finish;i++){
1375 if (m->control_pressed) { break; }
1379 int minFlowGram = 100000000;
1380 double minFlowValue = 1e8;
1381 change[i] = 0; //FALSE
1383 for(int j=0;j<nSeqsPerOTU[i];j++){
1384 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1387 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1388 vector<double> adF(nSeqsPerOTU[i]);
1389 vector<int> anL(nSeqsPerOTU[i]);
1391 for(int j=0;j<nSeqsPerOTU[i];j++){
1392 int index = cumNumSeqs[i] + j;
1393 int nI = seqIndex[index];
1394 int nIU = mapSeqToUnique[nI];
1397 for(k=0;k<position;k++){
1403 anL[position] = nIU;
1404 adF[position] = 0.0000;
1409 for(int j=0;j<nSeqsPerOTU[i];j++){
1410 int index = cumNumSeqs[i] + j;
1411 int nI = seqIndex[index];
1413 double tauValue = singleTau[seqNumber[index]];
1415 for(int k=0;k<position;k++){
1416 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1417 adF[k] += dist * tauValue;
1421 for(int j=0;j<position;j++){
1422 if(adF[j] < minFlowValue){
1424 minFlowValue = adF[j];
1428 if(centroids[i] != anL[minFlowGram]){
1430 centroids[i] = anL[minFlowGram];
1433 else if(centroids[i] != -1){
1439 catch(exception& e) {
1440 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1445 /**************************************************************************************************/
1447 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1450 int flowAValue = cent * numFlowCells;
1451 int flowBValue = flow * numFlowCells;
1455 for(int i=0;i<length;i++){
1456 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1461 return dist / (double)length;
1463 catch(exception& e) {
1464 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1469 /**************************************************************************************************/
1471 double ShhherCommand::getNewWeights(){
1474 double maxChange = 0;
1476 for(int i=0;i<numOTUs;i++){
1478 if (m->control_pressed) { break; }
1480 double difference = weight[i];
1483 for(int j=0;j<nSeqsPerOTU[i];j++){
1484 int index = cumNumSeqs[i] + j;
1485 double tauValue = singleTau[seqNumber[index]];
1486 weight[i] += tauValue;
1489 difference = fabs(weight[i] - difference);
1490 if(difference > maxChange){ maxChange = difference; }
1494 catch(exception& e) {
1495 m->errorOut(e, "ShhherCommand", "getNewWeights");
1500 /**************************************************************************************************/
1502 double ShhherCommand::getLikelihood(){
1506 vector<long double> P(numSeqs, 0);
1509 for(int i=0;i<numOTUs;i++){
1510 if(weight[i] > MIN_WEIGHT){
1516 for(int i=0;i<numOTUs;i++){
1518 if (m->control_pressed) { break; }
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 if (m->control_pressed) { break; }
1562 for(int j=i+1;j<numOTUs;j++){
1565 if(centroids[j] == centroids[i]){
1569 weight[i] += weight[j];
1577 catch(exception& e) {
1578 m->errorOut(e, "ShhherCommand", "checkCentroids");
1582 /**************************************************************************************************/
1586 void ShhherCommand::writeQualities(vector<int> otuCounts){
1589 string thisOutputDir = outputDir;
1590 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1591 string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.qual";
1593 ofstream qualityFile;
1594 m->openOutputFile(qualityFileName, qualityFile);
1596 qualityFile.setf(ios::fixed, ios::floatfield);
1597 qualityFile.setf(ios::showpoint);
1598 qualityFile << setprecision(6);
1600 vector<vector<int> > qualities(numOTUs);
1601 vector<double> pr(HOMOPS, 0);
1604 for(int i=0;i<numOTUs;i++){
1606 if (m->control_pressed) { break; }
1611 if(nSeqsPerOTU[i] > 0){
1612 qualities[i].assign(1024, -1);
1614 while(index < numFlowCells){
1615 double maxPrValue = 1e8;
1616 short maxPrIndex = -1;
1617 double count = 0.0000;
1619 pr.assign(HOMOPS, 0);
1621 for(int j=0;j<nSeqsPerOTU[i];j++){
1622 int lIndex = cumNumSeqs[i] + j;
1623 double tauValue = singleTau[seqNumber[lIndex]];
1624 int sequenceIndex = aaI[i][j];
1625 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1629 for(int s=0;s<HOMOPS;s++){
1630 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1634 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1635 maxPrValue = pr[maxPrIndex];
1637 if(count > MIN_COUNT){
1639 double norm = 0.0000;
1641 for(int s=0;s<HOMOPS;s++){
1642 norm += exp(-(pr[s] - maxPrValue));
1645 for(int s=1;s<=maxPrIndex;s++){
1647 double temp = 0.0000;
1649 U += exp(-(pr[s-1]-maxPrValue))/norm;
1657 temp = floor(-10 * temp);
1658 value = (int)floor(temp);
1659 if(value > 100){ value = 100; }
1661 qualities[i][base] = (int)value;
1671 if(otuCounts[i] > 0){
1672 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1674 int j=4; //need to get past the first four bases
1675 while(qualities[i][j] != -1){
1676 qualityFile << qualities[i][j] << ' ';
1679 qualityFile << endl;
1682 qualityFile.close();
1683 outputNames.push_back(qualityFileName);
1686 catch(exception& e) {
1687 m->errorOut(e, "ShhherCommand", "writeQualities");
1692 /**************************************************************************************************/
1694 void ShhherCommand::writeSequences(vector<int> otuCounts){
1696 string thisOutputDir = outputDir;
1697 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1698 string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.fasta";
1700 m->openOutputFile(fastaFileName, fastaFile);
1702 vector<string> names(numOTUs, "");
1704 for(int i=0;i<numOTUs;i++){
1706 if (m->control_pressed) { break; }
1708 int index = centroids[i];
1710 if(otuCounts[i] > 0){
1711 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
1715 for(int j=0;j<numFlowCells;j++){
1717 char base = flowOrder[j % 4];
1718 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
1723 fastaFile << newSeq.substr(4) << endl;
1728 outputNames.push_back(fastaFileName);
1730 if(compositeFASTAFileName != ""){
1731 m->appendFiles(fastaFileName, compositeFASTAFileName);
1734 catch(exception& e) {
1735 m->errorOut(e, "ShhherCommand", "writeSequences");
1740 /**************************************************************************************************/
1742 void ShhherCommand::writeNames(vector<int> otuCounts){
1744 string thisOutputDir = outputDir;
1745 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1746 string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.names";
1748 m->openOutputFile(nameFileName, nameFile);
1750 for(int i=0;i<numOTUs;i++){
1752 if (m->control_pressed) { break; }
1754 if(otuCounts[i] > 0){
1755 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
1757 for(int j=1;j<nSeqsPerOTU[i];j++){
1758 nameFile << ',' << seqNameVector[aaI[i][j]];
1765 outputNames.push_back(nameFileName);
1768 if(compositeNamesFileName != ""){
1769 m->appendFiles(nameFileName, compositeNamesFileName);
1772 catch(exception& e) {
1773 m->errorOut(e, "ShhherCommand", "writeNames");
1778 /**************************************************************************************************/
1780 void ShhherCommand::writeGroups(){
1782 string thisOutputDir = outputDir;
1783 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1784 string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1785 string groupFileName = fileRoot + "shhh.groups";
1787 m->openOutputFile(groupFileName, groupFile);
1789 for(int i=0;i<numSeqs;i++){
1790 if (m->control_pressed) { break; }
1791 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
1794 outputNames.push_back(groupFileName);
1797 catch(exception& e) {
1798 m->errorOut(e, "ShhherCommand", "writeGroups");
1803 /**************************************************************************************************/
1805 void ShhherCommand::writeClusters(vector<int> otuCounts){
1807 string thisOutputDir = outputDir;
1808 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1809 string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.counts";
1810 ofstream otuCountsFile;
1811 m->openOutputFile(otuCountsFileName, otuCountsFile);
1813 string bases = flowOrder;
1815 for(int i=0;i<numOTUs;i++){
1817 if (m->control_pressed) {
1820 //output the translated version of the centroid sequence for the otu
1821 if(otuCounts[i] > 0){
1822 int index = centroids[i];
1824 otuCountsFile << "ideal\t";
1825 for(int j=8;j<numFlowCells;j++){
1826 char base = bases[j % 4];
1827 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
1828 otuCountsFile << base;
1831 otuCountsFile << endl;
1833 for(int j=0;j<nSeqsPerOTU[i];j++){
1834 int sequence = aaI[i][j];
1835 otuCountsFile << seqNameVector[sequence] << '\t';
1839 for(int k=0;k<lengths[sequence];k++){
1840 char base = bases[k % 4];
1841 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
1843 for(int s=0;s<freq;s++){
1845 //otuCountsFile << base;
1848 otuCountsFile << newSeq.substr(4) << endl;
1850 otuCountsFile << endl;
1853 otuCountsFile.close();
1854 outputNames.push_back(otuCountsFileName);
1857 catch(exception& e) {
1858 m->errorOut(e, "ShhherCommand", "writeClusters");
1864 //**********************************************************************************************************************
1866 int ShhherCommand::execute(){
1868 if (abort == true) { return 0; }
1870 getSingleLookUp(); if (m->control_pressed) { return 0; }
1871 getJointLookUp(); if (m->control_pressed) { return 0; }
1873 int numFiles = flowFileVector.size();
1875 if (numFiles < processors) { processors = numFiles; }
1877 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1878 if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size()); }
1879 else { createProcesses(flowFileVector); } //each processor processes one file
1881 driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size());
1884 if(compositeFASTAFileName != ""){
1885 outputNames.push_back(compositeFASTAFileName);
1886 outputNames.push_back(compositeNamesFileName);
1889 m->mothurOutEndLine();
1890 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
1891 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
1892 m->mothurOutEndLine();
1896 catch(exception& e) {
1897 m->errorOut(e, "ShhherCommand", "execute");
1902 /**************************************************************************************************/
1904 int ShhherCommand::createProcesses(vector<string> filenames){
1906 vector<int> processIDS;
1911 if (filenames.size() < processors) { processors = filenames.size(); }
1913 //divide the groups between the processors
1914 vector<linePair> lines;
1915 int numFilesPerProcessor = filenames.size() / processors;
1916 for (int i = 0; i < processors; i++) {
1917 int startIndex = i * numFilesPerProcessor;
1918 int endIndex = (i+1) * numFilesPerProcessor;
1919 if(i == (processors - 1)){ endIndex = filenames.size(); }
1920 lines.push_back(linePair(startIndex, endIndex));
1923 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1925 //loop through and create all the processes you want
1926 while (process != processors) {
1930 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1932 }else if (pid == 0){
1933 num = driver(filenames, compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp", lines[process].start, lines[process].end);
1936 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1937 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1943 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[0].start, lines[0].end);
1945 //force parent to wait until all the processes are done
1946 for (int i=0;i<processIDS.size();i++) {
1947 int temp = processIDS[i];
1953 //////////////////////////////////////////////////////////////////////////////////////////////////////
1955 /////////////////////// NOT WORKING, ACCESS VIOLATION ON READ OF FLOWGRAMS IN THREAD /////////////////
1957 //////////////////////////////////////////////////////////////////////////////////////////////////////
1958 //Windows version shared memory, so be careful when passing variables through the shhhFlowsData struct.
1959 //Above fork() will clone, so memory is separate, but that's not the case with windows,
1960 //////////////////////////////////////////////////////////////////////////////////////////////////////
1962 vector<shhhFlowsData*> pDataArray;
1963 DWORD dwThreadIdArray[processors-1];
1964 HANDLE hThreadArray[processors-1];
1966 //Create processor worker threads.
1967 for( int i=0; i<processors-1; i++ ){
1968 // Allocate memory for thread data.
1969 string extension = "";
1970 if (i != 0) { extension = toString(i) + ".temp"; }
1972 shhhFlowsData* tempFlow = new shhhFlowsData(filenames, (compositeFASTAFileName + extension), (compositeNamesFileName + extension), outputDir, flowOrder, jointLookUp, singleLookUp, m, lines[i].start, lines[i].end, cutoff, sigma, minDelta, maxIters, i);
1973 pDataArray.push_back(tempFlow);
1974 processIDS.push_back(i);
1976 hThreadArray[i] = CreateThread(NULL, 0, ShhhFlowsThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
1979 //using the main process as a worker saves time and memory
1981 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[processors-1].start, lines[processors-1].end);
1983 //Wait until all threads have terminated.
1984 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
1986 //Close all thread handles and free memory allocations.
1987 for(int i=0; i < pDataArray.size(); i++){
1988 for(int j=0; j < pDataArray[i]->outputNames.size(); j++){ outputNames.push_back(pDataArray[i]->outputNames[j]); }
1989 CloseHandle(hThreadArray[i]);
1990 delete pDataArray[i];
1995 for (int i=0;i<processIDS.size();i++) {
1996 if (compositeFASTAFileName != "") {
1997 m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName);
1998 m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName);
1999 m->mothurRemove((compositeFASTAFileName + toString(processIDS[i]) + ".temp"));
2000 m->mothurRemove((compositeNamesFileName + toString(processIDS[i]) + ".temp"));
2007 catch(exception& e) {
2008 m->errorOut(e, "ShhherCommand", "createProcesses");
2012 /**************************************************************************************************/
2014 int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){
2017 for(int i=start;i<end;i++){
2019 if (m->control_pressed) { break; }
2021 string flowFileName = filenames[i];
2023 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n");
2024 m->mothurOut("Reading flowgrams...\n");
2026 vector<string> seqNameVector;
2027 vector<int> lengths;
2028 vector<short> flowDataIntI;
2029 vector<double> flowDataPrI;
2030 map<string, int> nameMap;
2031 vector<short> uniqueFlowgrams;
2032 vector<int> uniqueCount;
2033 vector<int> mapSeqToUnique;
2034 vector<int> mapUniqueToSeq;
2035 vector<int> uniqueLengths;
2038 int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells);
2040 if (m->control_pressed) { break; }
2042 m->mothurOut("Identifying unique flowgrams...\n");
2043 int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI);
2045 if (m->control_pressed) { break; }
2047 m->mothurOut("Calculating distances between flowgrams...\n");
2048 string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
2049 unsigned long long begTime = time(NULL);
2050 double begClock = clock();
2052 flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2054 m->mothurOutEndLine();
2055 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
2058 string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
2059 createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq);
2061 if (m->control_pressed) { break; }
2063 m->mothurOut("\nClustering flowgrams...\n");
2064 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
2065 cluster(listFileName, distFileName, namesFileName);
2067 if (m->control_pressed) { break; }
2069 vector<int> otuData;
2070 vector<int> cumNumSeqs;
2071 vector<int> nSeqsPerOTU;
2072 vector<vector<int> > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2073 vector<vector<int> > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2074 vector<int> seqNumber; //tMaster->anP: the sequence id number sorted by OTU
2075 vector<int> seqIndex; //tMaster->anI; the index that corresponds to seqNumber
2078 int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap);
2080 if (m->control_pressed) { break; }
2082 m->mothurRemove(distFileName);
2083 m->mothurRemove(namesFileName);
2084 m->mothurRemove(listFileName);
2086 vector<double> dist; //adDist - distance of sequences to centroids
2087 vector<short> change; //did the centroid sequence change? 0 = no; 1 = yes
2088 vector<int> centroids; //the representative flowgram for each cluster m
2089 vector<double> weight;
2090 vector<double> singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2091 vector<int> nSeqsBreaks;
2092 vector<int> nOTUsBreaks;
2094 dist.assign(numSeqs * numOTUs, 0);
2095 change.assign(numOTUs, 1);
2096 centroids.assign(numOTUs, -1);
2097 weight.assign(numOTUs, 0);
2098 singleTau.assign(numSeqs, 1.0);
2100 nSeqsBreaks.assign(2, 0);
2101 nOTUsBreaks.assign(2, 0);
2104 nSeqsBreaks[1] = numSeqs;
2105 nOTUsBreaks[1] = numOTUs;
2107 if (m->control_pressed) { break; }
2109 double maxDelta = 0;
2113 begTime = time(NULL);
2115 m->mothurOut("\nDenoising flowgrams...\n");
2116 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
2118 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
2120 if (m->control_pressed) { break; }
2122 double cycClock = clock();
2123 unsigned long long cycTime = time(NULL);
2124 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2126 if (m->control_pressed) { break; }
2128 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2130 if (m->control_pressed) { break; }
2132 maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight);
2134 if (m->control_pressed) { break; }
2136 double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight);
2138 if (m->control_pressed) { break; }
2140 checkCentroids(numOTUs, centroids, weight);
2142 if (m->control_pressed) { break; }
2144 calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths);
2146 if (m->control_pressed) { break; }
2150 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
2154 if (m->control_pressed) { break; }
2156 m->mothurOut("\nFinalizing...\n");
2157 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2159 if (m->control_pressed) { break; }
2161 setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI);
2163 if (m->control_pressed) { break; }
2165 vector<int> otuCounts(numOTUs, 0);
2166 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
2168 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2170 if (m->control_pressed) { break; }
2172 writeQualities(numOTUs, numFlowCells, flowFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; }
2173 writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, flowFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; }
2174 writeNames(thisCompositeNamesFileName, numOTUs, flowFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; }
2175 writeClusters(flowFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; }
2176 writeGroups(flowFileName, numSeqs, seqNameVector); if (m->control_pressed) { break; }
2178 m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
2181 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
2185 }catch(exception& e) {
2186 m->errorOut(e, "ShhherCommand", "driver");
2191 /**************************************************************************************************/
2192 int ShhherCommand::getFlowData(string filename, vector<string>& thisSeqNameVector, vector<int>& thisLengths, vector<short>& thisFlowDataIntI, map<string, int>& thisNameMap, int& numFlowCells){
2197 m->openInputFile(filename, flowFile);
2200 int currentNumFlowCells;
2202 thisSeqNameVector.clear();
2203 thisLengths.clear();
2204 thisFlowDataIntI.clear();
2205 thisNameMap.clear();
2207 flowFile >> numFlowCells;
2208 int index = 0;//pcluster
2209 while(!flowFile.eof()){
2211 if (m->control_pressed) { break; }
2213 flowFile >> seqName >> currentNumFlowCells;
2214 thisLengths.push_back(currentNumFlowCells);
2216 thisSeqNameVector.push_back(seqName);
2217 thisNameMap[seqName] = index++;//pcluster
2219 for(int i=0;i<numFlowCells;i++){
2220 flowFile >> intensity;
2221 if(intensity > 9.99) { intensity = 9.99; }
2222 int intI = int(100 * intensity + 0.0001);
2223 thisFlowDataIntI.push_back(intI);
2225 m->gobble(flowFile);
2229 int numSeqs = thisSeqNameVector.size();
2231 for(int i=0;i<numSeqs;i++){
2233 if (m->control_pressed) { break; }
2235 int iNumFlowCells = i * numFlowCells;
2236 for(int j=thisLengths[i];j<numFlowCells;j++){
2237 thisFlowDataIntI[iNumFlowCells + j] = 0;
2244 catch(exception& e) {
2245 m->errorOut(e, "ShhherCommand", "getFlowData");
2249 /**************************************************************************************************/
2251 int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int stopSeq, vector<int>& mapUniqueToSeq, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2254 ostringstream outStream;
2255 outStream.setf(ios::fixed, ios::floatfield);
2256 outStream.setf(ios::dec, ios::basefield);
2257 outStream.setf(ios::showpoint);
2258 outStream.precision(6);
2260 int begTime = time(NULL);
2261 double begClock = clock();
2263 for(int i=0;i<stopSeq;i++){
2265 if (m->control_pressed) { break; }
2267 for(int j=0;j<i;j++){
2268 float flowDistance = calcPairwiseDist(numFlowCells, mapUniqueToSeq[i], mapUniqueToSeq[j], mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2270 if(flowDistance < 1e-6){
2271 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
2273 else if(flowDistance <= cutoff){
2274 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
2278 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
2279 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2280 m->mothurOutEndLine();
2284 ofstream distFile(distFileName.c_str());
2285 distFile << outStream.str();
2288 if (m->control_pressed) {}
2290 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
2291 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2292 m->mothurOutEndLine();
2297 catch(exception& e) {
2298 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
2302 /**************************************************************************************************/
2304 float ShhherCommand::calcPairwiseDist(int numFlowCells, int seqA, int seqB, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2306 int minLength = lengths[mapSeqToUnique[seqA]];
2307 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
2309 int ANumFlowCells = seqA * numFlowCells;
2310 int BNumFlowCells = seqB * numFlowCells;
2314 for(int i=0;i<minLength;i++){
2316 if (m->control_pressed) { break; }
2318 int flowAIntI = flowDataIntI[ANumFlowCells + i];
2319 float flowAPrI = flowDataPrI[ANumFlowCells + i];
2321 int flowBIntI = flowDataIntI[BNumFlowCells + i];
2322 float flowBPrI = flowDataPrI[BNumFlowCells + i];
2323 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
2326 dist /= (float) minLength;
2329 catch(exception& e) {
2330 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
2335 /**************************************************************************************************/
2337 int ShhherCommand::getUniques(int numSeqs, int numFlowCells, vector<short>& uniqueFlowgrams, vector<int>& uniqueCount, vector<int>& uniqueLengths, vector<int>& mapSeqToUnique, vector<int>& mapUniqueToSeq, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2340 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
2341 uniqueCount.assign(numSeqs, 0); // anWeights
2342 uniqueLengths.assign(numSeqs, 0);
2343 mapSeqToUnique.assign(numSeqs, -1);
2344 mapUniqueToSeq.assign(numSeqs, -1);
2346 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
2348 for(int i=0;i<numSeqs;i++){
2350 if (m->control_pressed) { break; }
2354 vector<short> current(numFlowCells);
2355 for(int j=0;j<numFlowCells;j++){
2356 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
2359 for(int j=0;j<numUniques;j++){
2360 int offset = j * numFlowCells;
2364 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
2365 else { shorterLength = uniqueLengths[j]; }
2367 for(int k=0;k<shorterLength;k++){
2368 if(current[k] != uniqueFlowgrams[offset + k]){
2375 mapSeqToUnique[i] = j;
2378 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
2384 if(index == numUniques){
2385 uniqueLengths[numUniques] = lengths[i];
2386 uniqueCount[numUniques] = 1;
2387 mapSeqToUnique[i] = numUniques;//anMap
2388 mapUniqueToSeq[numUniques] = i;//anF
2390 for(int k=0;k<numFlowCells;k++){
2391 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
2392 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
2398 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
2399 uniqueLengths.resize(numUniques);
2401 flowDataPrI.resize(numSeqs * numFlowCells, 0);
2402 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
2406 catch(exception& e) {
2407 m->errorOut(e, "ShhherCommand", "getUniques");
2411 /**************************************************************************************************/
2412 int ShhherCommand::createNamesFile(int numSeqs, int numUniques, string filename, vector<string>& seqNameVector, vector<int>& mapSeqToUnique, vector<int>& mapUniqueToSeq){
2415 vector<string> duplicateNames(numUniques, "");
2416 for(int i=0;i<numSeqs;i++){
2417 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
2421 m->openOutputFile(filename, nameFile);
2423 for(int i=0;i<numUniques;i++){
2425 if (m->control_pressed) { break; }
2427 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2428 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2435 catch(exception& e) {
2436 m->errorOut(e, "ShhherCommand", "createNamesFile");
2440 //**********************************************************************************************************************
2442 int ShhherCommand::cluster(string filename, string distFileName, string namesFileName){
2445 ReadMatrix* read = new ReadColumnMatrix(distFileName);
2446 read->setCutoff(cutoff);
2448 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
2449 clusterNameMap->readMap();
2450 read->read(clusterNameMap);
2452 ListVector* list = read->getListVector();
2453 SparseMatrix* matrix = read->getMatrix();
2456 delete clusterNameMap;
2458 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
2460 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
2461 string tag = cluster->getTag();
2463 double clusterCutoff = cutoff;
2464 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
2466 if (m->control_pressed) { break; }
2468 cluster->update(clusterCutoff);
2471 list->setLabel(toString(cutoff));
2474 m->openOutputFile(filename, listFile);
2475 list->print(listFile);
2478 delete matrix; delete cluster; delete rabund; delete list;
2482 catch(exception& e) {
2483 m->errorOut(e, "ShhherCommand", "cluster");
2487 /**************************************************************************************************/
2489 int ShhherCommand::getOTUData(int numSeqs, string fileName, vector<int>& otuData,
2490 vector<int>& cumNumSeqs,
2491 vector<int>& nSeqsPerOTU,
2492 vector<vector<int> >& aaP, //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2493 vector<vector<int> >& aaI, //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2494 vector<int>& seqNumber, //tMaster->anP: the sequence id number sorted by OTU
2495 vector<int>& seqIndex,
2496 map<string, int>& nameMap){
2500 m->openInputFile(fileName, listFile);
2504 listFile >> label >> numOTUs;
2506 otuData.assign(numSeqs, 0);
2507 cumNumSeqs.assign(numOTUs, 0);
2508 nSeqsPerOTU.assign(numOTUs, 0);
2509 aaP.clear();aaP.resize(numOTUs);
2515 string singleOTU = "";
2517 for(int i=0;i<numOTUs;i++){
2519 if (m->control_pressed) { break; }
2521 listFile >> singleOTU;
2523 istringstream otuString(singleOTU);
2527 string seqName = "";
2529 for(int j=0;j<singleOTU.length();j++){
2530 char letter = otuString.get();
2536 map<string,int>::iterator nmIt = nameMap.find(seqName);
2537 int index = nmIt->second;
2539 nameMap.erase(nmIt);
2543 aaP[i].push_back(index);
2548 map<string,int>::iterator nmIt = nameMap.find(seqName);
2550 int index = nmIt->second;
2551 nameMap.erase(nmIt);
2555 aaP[i].push_back(index);
2560 sort(aaP[i].begin(), aaP[i].end());
2561 for(int j=0;j<nSeqsPerOTU[i];j++){
2562 seqNumber.push_back(aaP[i][j]);
2564 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
2565 aaP[i].push_back(0);
2571 for(int i=1;i<numOTUs;i++){
2572 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
2575 seqIndex = seqNumber;
2582 catch(exception& e) {
2583 m->errorOut(e, "ShhherCommand", "getOTUData");
2587 /**************************************************************************************************/
2589 int ShhherCommand::calcCentroidsDriver(int numOTUs,
2590 vector<int>& cumNumSeqs,
2591 vector<int>& nSeqsPerOTU,
2592 vector<int>& seqIndex,
2593 vector<short>& change, //did the centroid sequence change? 0 = no; 1 = yes
2594 vector<int>& centroids, //the representative flowgram for each cluster m
2595 vector<double>& singleTau, //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2596 vector<int>& mapSeqToUnique,
2597 vector<short>& uniqueFlowgrams,
2598 vector<short>& flowDataIntI,
2599 vector<int>& lengths,
2601 vector<int>& seqNumber){
2603 //this function gets the most likely homopolymer length at a flow position for a group of sequences
2608 for(int i=0;i<numOTUs;i++){
2610 if (m->control_pressed) { break; }
2614 int minFlowGram = 100000000;
2615 double minFlowValue = 1e8;
2616 change[i] = 0; //FALSE
2618 for(int j=0;j<nSeqsPerOTU[i];j++){
2619 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
2622 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
2623 vector<double> adF(nSeqsPerOTU[i]);
2624 vector<int> anL(nSeqsPerOTU[i]);
2626 for(int j=0;j<nSeqsPerOTU[i];j++){
2627 int index = cumNumSeqs[i] + j;
2628 int nI = seqIndex[index];
2629 int nIU = mapSeqToUnique[nI];
2632 for(k=0;k<position;k++){
2638 anL[position] = nIU;
2639 adF[position] = 0.0000;
2644 for(int j=0;j<nSeqsPerOTU[i];j++){
2645 int index = cumNumSeqs[i] + j;
2646 int nI = seqIndex[index];
2648 double tauValue = singleTau[seqNumber[index]];
2650 for(int k=0;k<position;k++){
2651 double dist = getDistToCentroid(anL[k], nI, lengths[nI], uniqueFlowgrams, flowDataIntI, numFlowCells);
2652 adF[k] += dist * tauValue;
2656 for(int j=0;j<position;j++){
2657 if(adF[j] < minFlowValue){
2659 minFlowValue = adF[j];
2663 if(centroids[i] != anL[minFlowGram]){
2665 centroids[i] = anL[minFlowGram];
2668 else if(centroids[i] != -1){
2676 catch(exception& e) {
2677 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
2681 /**************************************************************************************************/
2683 double ShhherCommand::getDistToCentroid(int cent, int flow, int length, vector<short>& uniqueFlowgrams,
2684 vector<short>& flowDataIntI, int numFlowCells){
2687 int flowAValue = cent * numFlowCells;
2688 int flowBValue = flow * numFlowCells;
2692 for(int i=0;i<length;i++){
2693 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
2698 return dist / (double)length;
2700 catch(exception& e) {
2701 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
2705 /**************************************************************************************************/
2707 double ShhherCommand::getNewWeights(int numOTUs, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU, vector<double>& singleTau, vector<int>& seqNumber, vector<double>& weight){
2710 double maxChange = 0;
2712 for(int i=0;i<numOTUs;i++){
2714 if (m->control_pressed) { break; }
2716 double difference = weight[i];
2719 for(int j=0;j<nSeqsPerOTU[i];j++){
2720 int index = cumNumSeqs[i] + j;
2721 double tauValue = singleTau[seqNumber[index]];
2722 weight[i] += tauValue;
2725 difference = fabs(weight[i] - difference);
2726 if(difference > maxChange){ maxChange = difference; }
2730 catch(exception& e) {
2731 m->errorOut(e, "ShhherCommand", "getNewWeights");
2736 /**************************************************************************************************/
2738 double ShhherCommand::getLikelihood(int numSeqs, int numOTUs, vector<int>& nSeqsPerOTU, vector<int>& seqNumber, vector<int>& cumNumSeqs, vector<int>& seqIndex, vector<double>& dist, vector<double>& weight){
2742 vector<long double> P(numSeqs, 0);
2745 for(int i=0;i<numOTUs;i++){
2746 if(weight[i] > MIN_WEIGHT){
2752 for(int i=0;i<numOTUs;i++){
2754 if (m->control_pressed) { break; }
2756 for(int j=0;j<nSeqsPerOTU[i];j++){
2757 int index = cumNumSeqs[i] + j;
2758 int nI = seqIndex[index];
2759 double singleDist = dist[seqNumber[index]];
2761 P[nI] += weight[i] * exp(-singleDist * sigma);
2765 for(int i=0;i<numSeqs;i++){
2766 if(P[i] == 0){ P[i] = DBL_EPSILON; }
2771 nLL = nLL -(double)numSeqs * log(sigma);
2775 catch(exception& e) {
2776 m->errorOut(e, "ShhherCommand", "getNewWeights");
2781 /**************************************************************************************************/
2783 int ShhherCommand::checkCentroids(int numOTUs, vector<int>& centroids, vector<double>& weight){
2785 vector<int> unique(numOTUs, 1);
2787 for(int i=0;i<numOTUs;i++){
2788 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
2793 for(int i=0;i<numOTUs;i++){
2795 if (m->control_pressed) { break; }
2798 for(int j=i+1;j<numOTUs;j++){
2801 if(centroids[j] == centroids[i]){
2805 weight[i] += weight[j];
2815 catch(exception& e) {
2816 m->errorOut(e, "ShhherCommand", "checkCentroids");
2820 /**************************************************************************************************/
2822 void ShhherCommand::calcNewDistances(int numSeqs, int numOTUs, vector<int>& nSeqsPerOTU, vector<double>& dist,
2823 vector<double>& weight, vector<short>& change, vector<int>& centroids,
2824 vector<vector<int> >& aaP, vector<double>& singleTau, vector<vector<int> >& aaI,
2825 vector<int>& seqNumber, vector<int>& seqIndex,
2826 vector<short>& uniqueFlowgrams,
2827 vector<short>& flowDataIntI, int numFlowCells, vector<int>& lengths){
2832 vector<double> newTau(numOTUs,0);
2833 vector<double> norms(numSeqs, 0);
2834 nSeqsPerOTU.assign(numOTUs, 0);
2836 for(int i=0;i<numSeqs;i++){
2838 if (m->control_pressed) { break; }
2840 int indexOffset = i * numOTUs;
2842 double offset = 1e8;
2844 for(int j=0;j<numOTUs;j++){
2846 if(weight[j] > MIN_WEIGHT && change[j] == 1){
2847 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i], uniqueFlowgrams, flowDataIntI, numFlowCells);
2850 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
2851 offset = dist[indexOffset + j];
2855 for(int j=0;j<numOTUs;j++){
2856 if(weight[j] > MIN_WEIGHT){
2857 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
2858 norms[i] += newTau[j];
2865 for(int j=0;j<numOTUs;j++){
2866 newTau[j] /= norms[i];
2869 for(int j=0;j<numOTUs;j++){
2870 if(newTau[j] > MIN_TAU){
2872 int oldTotal = total;
2876 singleTau.resize(total, 0);
2877 seqNumber.resize(total, 0);
2878 seqIndex.resize(total, 0);
2880 singleTau[oldTotal] = newTau[j];
2882 aaP[j][nSeqsPerOTU[j]] = oldTotal;
2883 aaI[j][nSeqsPerOTU[j]] = i;
2891 catch(exception& e) {
2892 m->errorOut(e, "ShhherCommand", "calcNewDistances");
2896 /**************************************************************************************************/
2898 int ShhherCommand::fill(int numOTUs, vector<int>& seqNumber, vector<int>& seqIndex, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU, vector<vector<int> >& aaP, vector<vector<int> >& aaI){
2901 for(int i=0;i<numOTUs;i++){
2903 if (m->control_pressed) { return 0; }
2905 cumNumSeqs[i] = index;
2906 for(int j=0;j<nSeqsPerOTU[i];j++){
2907 seqNumber[index] = aaP[i][j];
2908 seqIndex[index] = aaI[i][j];
2916 catch(exception& e) {
2917 m->errorOut(e, "ShhherCommand", "fill");
2921 /**************************************************************************************************/
2923 void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector<int>& seqNumber, vector<int>& seqIndex, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU,
2924 vector<int>& otuData, vector<double>& singleTau, vector<double>& dist, vector<vector<int> >& aaP, vector<vector<int> >& aaI){
2927 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
2929 for(int i=0;i<numOTUs;i++){
2931 if (m->control_pressed) { break; }
2933 for(int j=0;j<nSeqsPerOTU[i];j++){
2934 int index = cumNumSeqs[i] + j;
2935 double tauValue = singleTau[seqNumber[index]];
2936 int sIndex = seqIndex[index];
2937 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
2941 for(int i=0;i<numSeqs;i++){
2942 double maxTau = -1.0000;
2944 for(int j=0;j<numOTUs;j++){
2945 if(bigTauMatrix[i * numOTUs + j] > maxTau){
2946 maxTau = bigTauMatrix[i * numOTUs + j];
2951 otuData[i] = maxOTU;
2954 nSeqsPerOTU.assign(numOTUs, 0);
2956 for(int i=0;i<numSeqs;i++){
2957 int index = otuData[i];
2959 singleTau[i] = 1.0000;
2962 aaP[index][nSeqsPerOTU[index]] = i;
2963 aaI[index][nSeqsPerOTU[index]] = i;
2965 nSeqsPerOTU[index]++;
2968 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2970 catch(exception& e) {
2971 m->errorOut(e, "ShhherCommand", "setOTUs");
2975 /**************************************************************************************************/
2977 void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string filename, vector<int> otuCounts, vector<int>& nSeqsPerOTU, vector<int>& seqNumber,
2978 vector<double>& singleTau, vector<short>& flowDataIntI, vector<short>& uniqueFlowgrams, vector<int>& cumNumSeqs,
2979 vector<int>& mapUniqueToSeq, vector<string>& seqNameVector, vector<int>& centroids, vector<vector<int> >& aaI){
2982 string thisOutputDir = outputDir;
2983 if (outputDir == "") { thisOutputDir += m->hasPath(filename); }
2984 string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.qual";
2986 ofstream qualityFile;
2987 m->openOutputFile(qualityFileName, qualityFile);
2989 qualityFile.setf(ios::fixed, ios::floatfield);
2990 qualityFile.setf(ios::showpoint);
2991 qualityFile << setprecision(6);
2993 vector<vector<int> > qualities(numOTUs);
2994 vector<double> pr(HOMOPS, 0);
2997 for(int i=0;i<numOTUs;i++){
2999 if (m->control_pressed) { break; }
3004 if(nSeqsPerOTU[i] > 0){
3005 qualities[i].assign(1024, -1);
3007 while(index < numFlowCells){
3008 double maxPrValue = 1e8;
3009 short maxPrIndex = -1;
3010 double count = 0.0000;
3012 pr.assign(HOMOPS, 0);
3014 for(int j=0;j<nSeqsPerOTU[i];j++){
3015 int lIndex = cumNumSeqs[i] + j;
3016 double tauValue = singleTau[seqNumber[lIndex]];
3017 int sequenceIndex = aaI[i][j];
3018 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
3022 for(int s=0;s<HOMOPS;s++){
3023 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
3027 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
3028 maxPrValue = pr[maxPrIndex];
3030 if(count > MIN_COUNT){
3032 double norm = 0.0000;
3034 for(int s=0;s<HOMOPS;s++){
3035 norm += exp(-(pr[s] - maxPrValue));
3038 for(int s=1;s<=maxPrIndex;s++){
3040 double temp = 0.0000;
3042 U += exp(-(pr[s-1]-maxPrValue))/norm;
3050 temp = floor(-10 * temp);
3051 value = (int)floor(temp);
3052 if(value > 100){ value = 100; }
3054 qualities[i][base] = (int)value;
3064 if(otuCounts[i] > 0){
3065 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
3067 int j=4; //need to get past the first four bases
3068 while(qualities[i][j] != -1){
3069 qualityFile << qualities[i][j] << ' ';
3070 if (j > qualities[i].size()) { break; }
3073 qualityFile << endl;
3076 qualityFile.close();
3077 outputNames.push_back(qualityFileName);
3080 catch(exception& e) {
3081 m->errorOut(e, "ShhherCommand", "writeQualities");
3086 /**************************************************************************************************/
3088 void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTUs, int numFlowCells, string filename, vector<int> otuCounts, vector<short>& uniqueFlowgrams, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& centroids){
3090 string thisOutputDir = outputDir;
3091 if (outputDir == "") { thisOutputDir += m->hasPath(filename); }
3092 string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.fasta";
3094 m->openOutputFile(fastaFileName, fastaFile);
3096 vector<string> names(numOTUs, "");
3098 for(int i=0;i<numOTUs;i++){
3100 if (m->control_pressed) { break; }
3102 int index = centroids[i];
3104 if(otuCounts[i] > 0){
3105 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
3109 for(int j=0;j<numFlowCells;j++){
3111 char base = flowOrder[j % 4];
3112 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
3117 fastaFile << newSeq.substr(4) << endl;
3122 outputNames.push_back(fastaFileName);
3124 if(thisCompositeFASTAFileName != ""){
3125 m->appendFiles(fastaFileName, thisCompositeFASTAFileName);
3128 catch(exception& e) {
3129 m->errorOut(e, "ShhherCommand", "writeSequences");
3134 /**************************************************************************************************/
3136 void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string filename, vector<int> otuCounts, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& nSeqsPerOTU){
3138 string thisOutputDir = outputDir;
3139 if (outputDir == "") { thisOutputDir += m->hasPath(filename); }
3140 string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.names";
3142 m->openOutputFile(nameFileName, nameFile);
3144 for(int i=0;i<numOTUs;i++){
3146 if (m->control_pressed) { break; }
3148 if(otuCounts[i] > 0){
3149 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
3151 for(int j=1;j<nSeqsPerOTU[i];j++){
3152 nameFile << ',' << seqNameVector[aaI[i][j]];
3159 outputNames.push_back(nameFileName);
3162 if(thisCompositeNamesFileName != ""){
3163 m->appendFiles(nameFileName, thisCompositeNamesFileName);
3166 catch(exception& e) {
3167 m->errorOut(e, "ShhherCommand", "writeNames");
3172 /**************************************************************************************************/
3174 void ShhherCommand::writeGroups(string filename, int numSeqs, vector<string>& seqNameVector){
3176 string thisOutputDir = outputDir;
3177 if (outputDir == "") { thisOutputDir += m->hasPath(filename); }
3178 string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(filename));
3179 string groupFileName = fileRoot + "shhh.groups";
3181 m->openOutputFile(groupFileName, groupFile);
3183 for(int i=0;i<numSeqs;i++){
3184 if (m->control_pressed) { break; }
3185 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
3188 outputNames.push_back(groupFileName);
3191 catch(exception& e) {
3192 m->errorOut(e, "ShhherCommand", "writeGroups");
3197 /**************************************************************************************************/
3199 void ShhherCommand::writeClusters(string filename, int numOTUs, int numFlowCells, vector<int> otuCounts, vector<int>& centroids, vector<short>& uniqueFlowgrams, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& nSeqsPerOTU, vector<int>& lengths, vector<short>& flowDataIntI){
3201 string thisOutputDir = outputDir;
3202 if (outputDir == "") { thisOutputDir += m->hasPath(filename); }
3203 string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.counts";
3204 ofstream otuCountsFile;
3205 m->openOutputFile(otuCountsFileName, otuCountsFile);
3207 string bases = flowOrder;
3209 for(int i=0;i<numOTUs;i++){
3211 if (m->control_pressed) {
3214 //output the translated version of the centroid sequence for the otu
3215 if(otuCounts[i] > 0){
3216 int index = centroids[i];
3218 otuCountsFile << "ideal\t";
3219 for(int j=8;j<numFlowCells;j++){
3220 char base = bases[j % 4];
3221 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
3222 otuCountsFile << base;
3225 otuCountsFile << endl;
3227 for(int j=0;j<nSeqsPerOTU[i];j++){
3228 int sequence = aaI[i][j];
3229 otuCountsFile << seqNameVector[sequence] << '\t';
3233 for(int k=0;k<lengths[sequence];k++){
3234 char base = bases[k % 4];
3235 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
3237 for(int s=0;s<freq;s++){
3239 //otuCountsFile << base;
3242 otuCountsFile << newSeq.substr(4) << endl;
3244 otuCountsFile << endl;
3247 otuCountsFile.close();
3248 outputNames.push_back(otuCountsFileName);
3251 catch(exception& e) {
3252 m->errorOut(e, "ShhherCommand", "writeClusters");
3257 /**************************************************************************************************/
3259 void ShhherCommand::getSingleLookUp(){
3261 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
3262 singleLookUp.assign(HOMOPS * NUMBINS, 0);
3265 ifstream lookUpFile;
3266 m->openInputFile(lookupFileName, lookUpFile);
3268 for(int i=0;i<HOMOPS;i++){
3270 if (m->control_pressed) { break; }
3273 lookUpFile >> logFracFreq;
3275 for(int j=0;j<NUMBINS;j++) {
3276 lookUpFile >> singleLookUp[index];
3282 catch(exception& e) {
3283 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
3288 /**************************************************************************************************/
3290 void ShhherCommand::getJointLookUp(){
3293 // the most likely joint probability (-log) that two intenities have the same polymer length
3294 jointLookUp.resize(NUMBINS * NUMBINS, 0);
3296 for(int i=0;i<NUMBINS;i++){
3298 if (m->control_pressed) { break; }
3300 for(int j=0;j<NUMBINS;j++){
3302 double minSum = 100000000;
3304 for(int k=0;k<HOMOPS;k++){
3305 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
3307 if(sum < minSum) { minSum = sum; }
3309 jointLookUp[i * NUMBINS + j] = minSum;
3313 catch(exception& e) {
3314 m->errorOut(e, "ShhherCommand", "getJointLookUp");
3319 /**************************************************************************************************/
3321 double ShhherCommand::getProbIntensity(int intIntensity){
3324 double minNegLogProb = 100000000;
3327 for(int i=0;i<HOMOPS;i++){//loop signal strength
3329 if (m->control_pressed) { break; }
3331 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
3332 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
3335 return minNegLogProb;
3337 catch(exception& e) {
3338 m->errorOut(e, "ShhherCommand", "getProbIntensity");