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 plarge("large", "Number", "", "-1", "", "", "",false,false); parameters.push_back(plarge);
22 CommandParameter psigma("sigma", "Number", "", "60", "", "", "",false,false); parameters.push_back(psigma);
23 CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "",false,false); parameters.push_back(pmindelta);
24 CommandParameter porder("order", "String", "", "", "", "", "",false,false); parameters.push_back(porder);
25 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
26 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
28 vector<string> myArray;
29 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
33 m->errorOut(e, "ShhherCommand", "setParameters");
37 //**********************************************************************************************************************
38 string ShhherCommand::getHelpString(){
40 string helpString = "";
41 helpString += "The shhh.flows command reads a file containing flowgrams and creates a file of corrected sequences.\n";
45 m->errorOut(e, "ShhherCommand", "getHelpString");
49 //**********************************************************************************************************************
51 ShhherCommand::ShhherCommand(){
53 abort = true; calledHelp = true;
56 //initialize outputTypes
57 // vector<string> tempOutNames;
58 // outputTypes["pn.dist"] = tempOutNames;
62 m->errorOut(e, "ShhherCommand", "ShhherCommand");
67 //**********************************************************************************************************************
69 ShhherCommand::ShhherCommand(string option) {
73 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
74 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
78 abort = false; calledHelp = false;
80 //allow user to run help
81 if(option == "help") { help(); abort = true; calledHelp = true; }
82 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
85 vector<string> myArray = setParameters();
87 OptionParser parser(option);
88 map<string,string> parameters = parser.getParameters();
90 ValidParameters validParameter;
91 map<string,string>::iterator it;
93 //check to make sure all parameters are valid for command
94 for (it = parameters.begin(); it != parameters.end(); it++) {
95 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
98 //initialize outputTypes
99 vector<string> tempOutNames;
100 // outputTypes["pn.dist"] = tempOutNames;
101 // outputTypes["fasta"] = tempOutNames;
103 //if the user changes the input directory command factory will send this info to us in the output parameter
104 string inputDir = validParameter.validFile(parameters, "inputdir", false);
105 if (inputDir == "not found"){ inputDir = ""; }
108 it = parameters.find("flow");
109 //user has given a template file
110 if(it != parameters.end()){
111 path = m->hasPath(it->second);
112 //if the user has not given a path then, add inputdir. else leave path alone.
113 if (path == "") { parameters["flow"] = inputDir + it->second; }
116 it = parameters.find("lookup");
117 //user has given a template file
118 if(it != parameters.end()){
119 path = m->hasPath(it->second);
120 //if the user has not given a path then, add inputdir. else leave path alone.
121 if (path == "") { parameters["lookup"] = inputDir + it->second; }
124 it = parameters.find("file");
125 //user has given a template file
126 if(it != parameters.end()){
127 path = m->hasPath(it->second);
128 //if the user has not given a path then, add inputdir. else leave path alone.
129 if (path == "") { parameters["file"] = inputDir + it->second; }
133 //if the user changes the output directory command factory will send this info to us in the output parameter
134 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
136 //check for required parameters
137 flowFileName = validParameter.validFile(parameters, "flow", true);
138 flowFilesFileName = validParameter.validFile(parameters, "file", true);
139 if (flowFileName == "not found" && flowFilesFileName == "not found") {
140 m->mothurOut("values for either flow or file must be provided for the shhh.flows command.");
141 m->mothurOutEndLine();
144 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
146 if(flowFileName != "not found"){
147 compositeFASTAFileName = "";
148 compositeNamesFileName = "";
153 string thisoutputDir = outputDir;
154 if (outputDir == "") { thisoutputDir = m->hasPath(flowFilesFileName); } //if user entered a file with a path then preserve it
156 //we want to rip off .files, and also .flow if its there
157 string fileroot = m->getRootName(m->getSimpleName(flowFilesFileName));
158 if (fileroot[fileroot.length()-1] == '.') { fileroot = fileroot.substr(0, fileroot.length()-1); } //rip off dot
159 string extension = m->getExtension(fileroot);
160 if (extension == ".flow") { fileroot = m->getRootName(fileroot); }
161 else { fileroot += "."; } //add back if needed
163 compositeFASTAFileName = thisoutputDir + fileroot + "shhh.fasta";
164 m->openOutputFile(compositeFASTAFileName, temp);
167 compositeNamesFileName = thisoutputDir + fileroot + "shhh.names";
168 m->openOutputFile(compositeNamesFileName, temp);
172 if(flowFilesFileName != "not found"){
175 ifstream flowFilesFile;
176 m->openInputFile(flowFilesFileName, flowFilesFile);
177 while(flowFilesFile){
178 fName = m->getline(flowFilesFile);
180 //test if file is valid
182 int ableToOpen = m->openInputFile(fName, in, "noerror");
184 if (ableToOpen == 1) {
185 if (inputDir != "") { //default path is set
186 string tryPath = inputDir + fName;
187 m->mothurOut("Unable to open " + fName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
189 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
195 if (ableToOpen == 1) {
196 if (m->getDefaultPath() != "") { //default path is set
197 string tryPath = m->getDefaultPath() + m->getSimpleName(fName);
198 m->mothurOut("Unable to open " + fName + ". Trying default " + tryPath); m->mothurOutEndLine();
200 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
206 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
207 if (ableToOpen == 1) {
208 string exepath = m->argv;
209 string tempPath = exepath;
210 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
211 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
213 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(fName);
214 m->mothurOut("Unable to open " + fName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
216 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
221 if (ableToOpen == 1) { m->mothurOut("Unable to open " + fName + ". Disregarding. "); m->mothurOutEndLine(); }
222 else { flowFileVector.push_back(fName); }
223 m->gobble(flowFilesFile);
225 flowFilesFile.close();
226 if (flowFileVector.size() == 0) { m->mothurOut("[ERROR]: no valid files."); m->mothurOutEndLine(); abort = true; }
229 if (outputDir == "") { outputDir = m->hasPath(flowFileName); }
230 flowFileVector.push_back(flowFileName);
233 //check for optional parameter and set defaults
234 // ...at some point should added some additional type checking...
236 temp = validParameter.validFile(parameters, "lookup", true);
237 if (temp == "not found") {
238 lookupFileName = "LookUp_Titanium.pat";
242 ableToOpen = m->openInputFile(lookupFileName, in, "noerror");
245 //if you can't open it, try input location
246 if (ableToOpen == 1) {
247 if (inputDir != "") { //default path is set
248 string tryPath = inputDir + lookupFileName;
249 m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
251 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
253 lookupFileName = tryPath;
257 //if you can't open it, try default location
258 if (ableToOpen == 1) {
259 if (m->getDefaultPath() != "") { //default path is set
260 string tryPath = m->getDefaultPath() + m->getSimpleName(lookupFileName);
261 m->mothurOut("Unable to open " + lookupFileName + ". Trying default " + tryPath); m->mothurOutEndLine();
263 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
265 lookupFileName = tryPath;
269 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
270 if (ableToOpen == 1) {
271 string exepath = m->argv;
272 string tempPath = exepath;
273 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
274 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
276 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
277 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
279 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
281 lookupFileName = tryPath;
284 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
286 else if(temp == "not open") {
288 lookupFileName = validParameter.validFile(parameters, "lookup", false);
290 //if you can't open it its not inputDir, try mothur excutable location
291 string exepath = m->argv;
292 string tempPath = exepath;
293 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
294 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
296 string tryPath = m->getFullPathName(exepath) + lookupFileName;
297 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
299 int ableToOpen = m->openInputFile(tryPath, in2, "noerror");
301 lookupFileName = tryPath;
303 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
304 }else { lookupFileName = temp; }
306 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
307 m->setProcessors(temp);
308 m->mothurConvert(temp, processors);
310 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
311 m->mothurConvert(temp, cutoff);
313 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
314 m->mothurConvert(temp, minDelta);
316 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
317 m->mothurConvert(temp, maxIters);
319 temp = validParameter.validFile(parameters, "large", false); if (temp == "not found"){ temp = "0"; }
320 m->mothurConvert(temp, largeSize);
321 if (largeSize != 0) { large = true; }
322 else { large = false; }
323 if (largeSize < 0) { m->mothurOut("The value of the large cannot be negative.\n"); }
326 if (large) { m->mothurOut("The large parameter is not available with the MPI-Enabled version.\n"); large=false; }
330 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
331 m->mothurConvert(temp, sigma);
333 flowOrder = validParameter.validFile(parameters, "order", false);
334 if (flowOrder == "not found"){ flowOrder = "TACG"; }
335 else if(flowOrder.length() != 4){
336 m->mothurOut("The value of the order option must be four bases long\n");
344 catch(exception& e) {
345 m->errorOut(e, "ShhherCommand", "ShhherCommand");
349 //**********************************************************************************************************************
351 int ShhherCommand::execute(){
353 if (abort == true) { if (calledHelp) { return 0; } return 2; }
360 for(int i=1;i<ncpus;i++){
361 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
363 if(abort == 1){ return 0; }
367 m->mothurOut("\nGetting preliminary data...\n");
368 getSingleLookUp(); if (m->control_pressed) { return 0; }
369 getJointLookUp(); if (m->control_pressed) { return 0; }
371 vector<string> flowFileVector;
372 if(flowFilesFileName != "not found"){
375 ifstream flowFilesFile;
376 m->openInputFile(flowFilesFileName, flowFilesFile);
377 while(flowFilesFile){
378 fName = m->getline(flowFilesFile);
379 flowFileVector.push_back(fName);
380 m->gobble(flowFilesFile);
384 flowFileVector.push_back(flowFileName);
387 int numFiles = flowFileVector.size();
389 for(int i=1;i<ncpus;i++){
390 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
393 for(int i=0;i<numFiles;i++){
395 if (m->control_pressed) { break; }
397 double begClock = clock();
398 unsigned long long begTime = time(NULL);
400 flowFileName = flowFileVector[i];
402 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
403 m->mothurOut("Reading flowgrams...\n");
406 if (m->control_pressed) { break; }
408 m->mothurOut("Identifying unique flowgrams...\n");
411 if (m->control_pressed) { break; }
413 m->mothurOut("Calculating distances between flowgrams...\n");
415 strcpy(fileName, flowFileName.c_str());
417 for(int i=1;i<ncpus;i++){
418 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
420 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
421 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
422 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
423 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
424 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
425 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
426 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
427 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
428 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
429 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
432 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
434 if (m->control_pressed) { break; }
437 for(int i=1;i<ncpus;i++){
438 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
440 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
441 m->mothurRemove((distFileName + ".temp." + toString(i)));
444 string namesFileName = createNamesFile();
446 if (m->control_pressed) { break; }
448 m->mothurOut("\nClustering flowgrams...\n");
449 string listFileName = cluster(distFileName, namesFileName);
451 if (m->control_pressed) { break; }
455 getOTUData(listFileName);
457 m->mothurRemove(distFileName);
458 m->mothurRemove(namesFileName);
459 m->mothurRemove(listFileName);
461 if (m->control_pressed) { break; }
465 if (m->control_pressed) { break; }
468 for(int i=1;i<ncpus;i++){
469 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
470 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
471 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
472 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
475 if (m->control_pressed) { break; }
480 int numOTUsOnCPU = numOTUs / ncpus;
481 int numSeqsOnCPU = numSeqs / ncpus;
482 m->mothurOut("\nDenoising flowgrams...\n");
483 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
485 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
487 double cycClock = clock();
488 unsigned long long cycTime = time(NULL);
491 if (m->control_pressed) { break; }
493 int total = singleTau.size();
494 for(int i=1;i<ncpus;i++){
495 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
496 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
497 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
499 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
500 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
501 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
502 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
503 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
506 calcCentroidsDriver(0, numOTUsOnCPU);
508 for(int i=1;i<ncpus;i++){
509 int otuStart = i * numOTUs / ncpus;
510 int otuStop = (i + 1) * numOTUs / ncpus;
512 vector<int> tempCentroids(numOTUs, 0);
513 vector<short> tempChange(numOTUs, 0);
515 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
516 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
518 for(int j=otuStart;j<otuStop;j++){
519 centroids[j] = tempCentroids[j];
520 change[j] = tempChange[j];
524 maxDelta = getNewWeights(); if (m->control_pressed) { break; }
525 double nLL = getLikelihood(); if (m->control_pressed) { break; }
526 checkCentroids(); if (m->control_pressed) { break; }
528 for(int i=1;i<ncpus;i++){
529 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
530 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
531 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
534 calcNewDistancesParent(0, numSeqsOnCPU);
536 total = singleTau.size();
538 for(int i=1;i<ncpus;i++){
540 int seqStart = i * numSeqs / ncpus;
541 int seqStop = (i + 1) * numSeqs / ncpus;
543 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
545 vector<int> childSeqIndex(childTotal, 0);
546 vector<double> childSingleTau(childTotal, 0);
547 vector<double> childDist(numSeqs * numOTUs, 0);
548 vector<int> otuIndex(childTotal, 0);
550 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
551 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
552 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
553 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
555 int oldTotal = total;
557 singleTau.resize(total, 0);
558 seqIndex.resize(total, 0);
559 seqNumber.resize(total, 0);
563 for(int j=oldTotal;j<total;j++){
564 int otuI = otuIndex[childIndex];
565 int seqI = childSeqIndex[childIndex];
567 singleTau[j] = childSingleTau[childIndex];
569 aaP[otuI][nSeqsPerOTU[otuI]] = j;
570 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
575 int index = seqStart * numOTUs;
576 for(int j=seqStart;j<seqStop;j++){
577 for(int k=0;k<numOTUs;k++){
578 dist[index] = childDist[index];
586 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
588 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
590 for(int i=1;i<ncpus;i++){
591 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
596 for(int i=1;i<ncpus;i++){
597 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
603 if (m->control_pressed) { break; }
605 m->mothurOut("\nFinalizing...\n");
608 if (m->control_pressed) { break; }
612 vector<int> otuCounts(numOTUs, 0);
613 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
614 calcCentroidsDriver(0, numOTUs);
616 if (m->control_pressed) { break; }
618 writeQualities(otuCounts); if (m->control_pressed) { break; }
619 writeSequences(otuCounts); if (m->control_pressed) { break; }
620 writeNames(otuCounts); if (m->control_pressed) { break; }
621 writeClusters(otuCounts); if (m->control_pressed) { break; }
622 writeGroups(); if (m->control_pressed) { break; }
625 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
631 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
632 if(abort){ return 0; }
635 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
637 for(int i=0;i<numFiles;i++){
639 if (m->control_pressed) { break; }
641 //Now into the pyrodist part
645 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
646 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
647 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
648 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
650 flowDataIntI.resize(numSeqs * numFlowCells);
651 flowDataPrI.resize(numSeqs * numFlowCells);
652 mapUniqueToSeq.resize(numSeqs);
653 mapSeqToUnique.resize(numSeqs);
654 lengths.resize(numSeqs);
655 jointLookUp.resize(NUMBINS * NUMBINS);
657 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
658 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
659 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
660 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
661 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
662 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
663 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
665 flowFileName = string(fileName);
666 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
667 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
669 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
671 if (m->control_pressed) { break; }
674 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
676 //Now into the pyronoise part
677 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
679 singleLookUp.resize(HOMOPS * NUMBINS);
680 uniqueFlowgrams.resize(numUniques * numFlowCells);
681 weight.resize(numOTUs);
682 centroids.resize(numOTUs);
683 change.resize(numOTUs);
684 dist.assign(numOTUs * numSeqs, 0);
685 nSeqsPerOTU.resize(numOTUs);
686 cumNumSeqs.resize(numOTUs);
688 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
689 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
690 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
692 int startOTU = pid * numOTUs / ncpus;
693 int endOTU = (pid + 1) * numOTUs / ncpus;
695 int startSeq = pid * numSeqs / ncpus;
696 int endSeq = (pid + 1) * numSeqs /ncpus;
702 if (m->control_pressed) { break; }
704 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
705 singleTau.assign(total, 0.0000);
706 seqNumber.assign(total, 0);
707 seqIndex.assign(total, 0);
709 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
710 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
711 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
712 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
713 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
714 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
715 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
717 calcCentroidsDriver(startOTU, endOTU);
719 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
720 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
722 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
723 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
724 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
726 vector<int> otuIndex(total, 0);
727 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
728 total = otuIndex.size();
730 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
731 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
732 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
733 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
734 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
736 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
741 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
743 MPI_Barrier(MPI_COMM_WORLD);
746 if(compositeFASTAFileName != ""){
747 outputNames.push_back(compositeFASTAFileName);
748 outputNames.push_back(compositeNamesFileName);
751 m->mothurOutEndLine();
752 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
753 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
754 m->mothurOutEndLine();
759 catch(exception& e) {
760 m->errorOut(e, "ShhherCommand", "execute");
764 /**************************************************************************************************/
765 string ShhherCommand::createNamesFile(){
768 vector<string> duplicateNames(numUniques, "");
769 for(int i=0;i<numSeqs;i++){
770 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
773 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
776 m->openOutputFile(nameFileName, nameFile);
778 for(int i=0;i<numUniques;i++){
780 if (m->control_pressed) { break; }
782 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
783 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
789 catch(exception& e) {
790 m->errorOut(e, "ShhherCommand", "createNamesFile");
794 /**************************************************************************************************/
796 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
798 ostringstream outStream;
799 outStream.setf(ios::fixed, ios::floatfield);
800 outStream.setf(ios::dec, ios::basefield);
801 outStream.setf(ios::showpoint);
802 outStream.precision(6);
804 int begTime = time(NULL);
805 double begClock = clock();
807 for(int i=startSeq;i<stopSeq;i++){
809 if (m->control_pressed) { break; }
811 for(int j=0;j<i;j++){
812 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
814 if(flowDistance < 1e-6){
815 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
817 else if(flowDistance <= cutoff){
818 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
822 m->mothurOut(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
826 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
827 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
829 if (m->control_pressed) { return fDistFileName; }
831 m->mothurOut(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
833 ofstream distFile(fDistFileName.c_str());
834 distFile << outStream.str();
837 return fDistFileName;
839 catch(exception& e) {
840 m->errorOut(e, "ShhherCommand", "flowDistMPI");
844 /**************************************************************************************************/
846 void ShhherCommand::getOTUData(string listFileName){
850 m->openInputFile(listFileName, listFile);
853 listFile >> label >> numOTUs;
855 otuData.assign(numSeqs, 0);
856 cumNumSeqs.assign(numOTUs, 0);
857 nSeqsPerOTU.assign(numOTUs, 0);
858 aaP.clear();aaP.resize(numOTUs);
864 string singleOTU = "";
866 for(int i=0;i<numOTUs;i++){
868 if (m->control_pressed) { break; }
870 listFile >> singleOTU;
872 istringstream otuString(singleOTU);
878 for(int j=0;j<singleOTU.length();j++){
879 char letter = otuString.get();
885 map<string,int>::iterator nmIt = nameMap.find(seqName);
886 int index = nmIt->second;
892 aaP[i].push_back(index);
897 map<string,int>::iterator nmIt = nameMap.find(seqName);
899 int index = nmIt->second;
904 aaP[i].push_back(index);
909 sort(aaP[i].begin(), aaP[i].end());
910 for(int j=0;j<nSeqsPerOTU[i];j++){
911 seqNumber.push_back(aaP[i][j]);
913 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
920 for(int i=1;i<numOTUs;i++){
921 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
924 seqIndex = seqNumber;
929 catch(exception& e) {
930 m->errorOut(e, "ShhherCommand", "getOTUData");
935 /**************************************************************************************************/
937 void ShhherCommand::initPyroCluster(){
939 if (numOTUs < processors) { processors = 1; }
941 dist.assign(numSeqs * numOTUs, 0);
942 change.assign(numOTUs, 1);
943 centroids.assign(numOTUs, -1);
944 weight.assign(numOTUs, 0);
945 singleTau.assign(numSeqs, 1.0);
947 nSeqsBreaks.assign(processors+1, 0);
948 nOTUsBreaks.assign(processors+1, 0);
951 for(int i=0;i<processors;i++){
952 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
953 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
955 nSeqsBreaks[processors] = numSeqs;
956 nOTUsBreaks[processors] = numOTUs;
958 catch(exception& e) {
959 m->errorOut(e, "ShhherCommand", "initPyroCluster");
964 /**************************************************************************************************/
966 void ShhherCommand::fill(){
969 for(int i=0;i<numOTUs;i++){
971 if (m->control_pressed) { break; }
973 cumNumSeqs[i] = index;
974 for(int j=0;j<nSeqsPerOTU[i];j++){
975 seqNumber[index] = aaP[i][j];
976 seqIndex[index] = aaI[i][j];
982 catch(exception& e) {
983 m->errorOut(e, "ShhherCommand", "fill");
988 /**************************************************************************************************/
990 void ShhherCommand::getFlowData(){
993 m->openInputFile(flowFileName, flowFile);
996 seqNameVector.clear();
998 flowDataIntI.clear();
1002 int currentNumFlowCells;
1006 flowFile >> numFlowCells;
1007 int index = 0;//pcluster
1008 while(!flowFile.eof()){
1010 if (m->control_pressed) { break; }
1012 flowFile >> seqName >> currentNumFlowCells;
1013 lengths.push_back(currentNumFlowCells);
1015 seqNameVector.push_back(seqName);
1016 nameMap[seqName] = index++;//pcluster
1018 for(int i=0;i<numFlowCells;i++){
1019 flowFile >> intensity;
1020 if(intensity > 9.99) { intensity = 9.99; }
1021 int intI = int(100 * intensity + 0.0001);
1022 flowDataIntI.push_back(intI);
1024 m->gobble(flowFile);
1028 numSeqs = seqNameVector.size();
1030 for(int i=0;i<numSeqs;i++){
1032 if (m->control_pressed) { break; }
1034 int iNumFlowCells = i * numFlowCells;
1035 for(int j=lengths[i];j<numFlowCells;j++){
1036 flowDataIntI[iNumFlowCells + j] = 0;
1041 catch(exception& e) {
1042 m->errorOut(e, "ShhherCommand", "getFlowData");
1046 /**************************************************************************************************/
1047 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1050 vector<double> newTau(numOTUs,0);
1051 vector<double> norms(numSeqs, 0);
1056 for(int i=startSeq;i<stopSeq;i++){
1058 if (m->control_pressed) { break; }
1060 double offset = 1e8;
1061 int indexOffset = i * numOTUs;
1063 for(int j=0;j<numOTUs;j++){
1065 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1066 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1068 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1069 offset = dist[indexOffset + j];
1073 for(int j=0;j<numOTUs;j++){
1074 if(weight[j] > MIN_WEIGHT){
1075 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1076 norms[i] += newTau[j];
1083 for(int j=0;j<numOTUs;j++){
1085 newTau[j] /= norms[i];
1087 if(newTau[j] > MIN_TAU){
1088 otuIndex.push_back(j);
1089 seqIndex.push_back(i);
1090 singleTau.push_back(newTau[j]);
1096 catch(exception& e) {
1097 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1102 /**************************************************************************************************/
1104 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1109 vector<double> newTau(numOTUs,0);
1110 vector<double> norms(numSeqs, 0);
1111 nSeqsPerOTU.assign(numOTUs, 0);
1113 for(int i=startSeq;i<stopSeq;i++){
1115 if (m->control_pressed) { break; }
1117 int indexOffset = i * numOTUs;
1119 double offset = 1e8;
1121 for(int j=0;j<numOTUs;j++){
1123 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1124 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1127 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1128 offset = dist[indexOffset + j];
1132 for(int j=0;j<numOTUs;j++){
1133 if(weight[j] > MIN_WEIGHT){
1134 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1135 norms[i] += newTau[j];
1142 for(int j=0;j<numOTUs;j++){
1143 newTau[j] /= norms[i];
1146 for(int j=0;j<numOTUs;j++){
1147 if(newTau[j] > MIN_TAU){
1149 int oldTotal = total;
1153 singleTau.resize(total, 0);
1154 seqNumber.resize(total, 0);
1155 seqIndex.resize(total, 0);
1157 singleTau[oldTotal] = newTau[j];
1159 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1160 aaI[j][nSeqsPerOTU[j]] = i;
1168 catch(exception& e) {
1169 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1174 /**************************************************************************************************/
1176 void ShhherCommand::setOTUs(){
1179 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1181 for(int i=0;i<numOTUs;i++){
1183 if (m->control_pressed) { break; }
1185 for(int j=0;j<nSeqsPerOTU[i];j++){
1186 int index = cumNumSeqs[i] + j;
1187 double tauValue = singleTau[seqNumber[index]];
1188 int sIndex = seqIndex[index];
1189 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1193 for(int i=0;i<numSeqs;i++){
1194 double maxTau = -1.0000;
1196 for(int j=0;j<numOTUs;j++){
1197 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1198 maxTau = bigTauMatrix[i * numOTUs + j];
1203 otuData[i] = maxOTU;
1206 nSeqsPerOTU.assign(numOTUs, 0);
1208 for(int i=0;i<numSeqs;i++){
1209 int index = otuData[i];
1211 singleTau[i] = 1.0000;
1214 aaP[index][nSeqsPerOTU[index]] = i;
1215 aaI[index][nSeqsPerOTU[index]] = i;
1217 nSeqsPerOTU[index]++;
1221 catch(exception& e) {
1222 m->errorOut(e, "ShhherCommand", "setOTUs");
1227 /**************************************************************************************************/
1229 void ShhherCommand::getUniques(){
1234 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
1235 uniqueCount.assign(numSeqs, 0); // anWeights
1236 uniqueLengths.assign(numSeqs, 0);
1237 mapSeqToUnique.assign(numSeqs, -1);
1238 mapUniqueToSeq.assign(numSeqs, -1);
1240 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
1242 for(int i=0;i<numSeqs;i++){
1244 if (m->control_pressed) { break; }
1248 vector<short> current(numFlowCells);
1249 for(int j=0;j<numFlowCells;j++){
1250 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
1253 for(int j=0;j<numUniques;j++){
1254 int offset = j * numFlowCells;
1258 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
1259 else { shorterLength = uniqueLengths[j]; }
1261 for(int k=0;k<shorterLength;k++){
1262 if(current[k] != uniqueFlowgrams[offset + k]){
1269 mapSeqToUnique[i] = j;
1272 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
1278 if(index == numUniques){
1279 uniqueLengths[numUniques] = lengths[i];
1280 uniqueCount[numUniques] = 1;
1281 mapSeqToUnique[i] = numUniques;//anMap
1282 mapUniqueToSeq[numUniques] = i;//anF
1284 for(int k=0;k<numFlowCells;k++){
1285 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
1286 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
1292 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
1293 uniqueLengths.resize(numUniques);
1295 flowDataPrI.resize(numSeqs * numFlowCells, 0);
1296 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
1298 catch(exception& e) {
1299 m->errorOut(e, "ShhherCommand", "getUniques");
1304 /**************************************************************************************************/
1306 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
1308 int minLength = lengths[mapSeqToUnique[seqA]];
1309 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
1311 int ANumFlowCells = seqA * numFlowCells;
1312 int BNumFlowCells = seqB * numFlowCells;
1316 for(int i=0;i<minLength;i++){
1318 if (m->control_pressed) { break; }
1320 int flowAIntI = flowDataIntI[ANumFlowCells + i];
1321 float flowAPrI = flowDataPrI[ANumFlowCells + i];
1323 int flowBIntI = flowDataIntI[BNumFlowCells + i];
1324 float flowBPrI = flowDataPrI[BNumFlowCells + i];
1325 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
1328 dist /= (float) minLength;
1331 catch(exception& e) {
1332 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
1337 //**********************************************************************************************************************/
1339 string ShhherCommand::cluster(string distFileName, string namesFileName){
1342 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1343 read->setCutoff(cutoff);
1345 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1346 clusterNameMap->readMap();
1347 read->read(clusterNameMap);
1349 ListVector* list = read->getListVector();
1350 SparseMatrix* matrix = read->getMatrix();
1353 delete clusterNameMap;
1355 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1357 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1358 string tag = cluster->getTag();
1360 double clusterCutoff = cutoff;
1361 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1363 if (m->control_pressed) { break; }
1365 cluster->update(clusterCutoff);
1368 list->setLabel(toString(cutoff));
1370 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
1372 m->openOutputFile(listFileName, listFile);
1373 list->print(listFile);
1376 delete matrix; delete cluster; delete rabund; delete list;
1378 return listFileName;
1380 catch(exception& e) {
1381 m->errorOut(e, "ShhherCommand", "cluster");
1386 /**************************************************************************************************/
1388 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1390 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1395 for(int i=start;i<finish;i++){
1397 if (m->control_pressed) { break; }
1401 int minFlowGram = 100000000;
1402 double minFlowValue = 1e8;
1403 change[i] = 0; //FALSE
1405 for(int j=0;j<nSeqsPerOTU[i];j++){
1406 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1409 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1410 vector<double> adF(nSeqsPerOTU[i]);
1411 vector<int> anL(nSeqsPerOTU[i]);
1413 for(int j=0;j<nSeqsPerOTU[i];j++){
1414 int index = cumNumSeqs[i] + j;
1415 int nI = seqIndex[index];
1416 int nIU = mapSeqToUnique[nI];
1419 for(k=0;k<position;k++){
1425 anL[position] = nIU;
1426 adF[position] = 0.0000;
1431 for(int j=0;j<nSeqsPerOTU[i];j++){
1432 int index = cumNumSeqs[i] + j;
1433 int nI = seqIndex[index];
1435 double tauValue = singleTau[seqNumber[index]];
1437 for(int k=0;k<position;k++){
1438 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1439 adF[k] += dist * tauValue;
1443 for(int j=0;j<position;j++){
1444 if(adF[j] < minFlowValue){
1446 minFlowValue = adF[j];
1450 if(centroids[i] != anL[minFlowGram]){
1452 centroids[i] = anL[minFlowGram];
1455 else if(centroids[i] != -1){
1461 catch(exception& e) {
1462 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1467 /**************************************************************************************************/
1469 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1472 int flowAValue = cent * numFlowCells;
1473 int flowBValue = flow * numFlowCells;
1477 for(int i=0;i<length;i++){
1478 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1483 return dist / (double)length;
1485 catch(exception& e) {
1486 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1491 /**************************************************************************************************/
1493 double ShhherCommand::getNewWeights(){
1496 double maxChange = 0;
1498 for(int i=0;i<numOTUs;i++){
1500 if (m->control_pressed) { break; }
1502 double difference = weight[i];
1505 for(int j=0;j<nSeqsPerOTU[i];j++){
1506 int index = cumNumSeqs[i] + j;
1507 double tauValue = singleTau[seqNumber[index]];
1508 weight[i] += tauValue;
1511 difference = fabs(weight[i] - difference);
1512 if(difference > maxChange){ maxChange = difference; }
1516 catch(exception& e) {
1517 m->errorOut(e, "ShhherCommand", "getNewWeights");
1522 /**************************************************************************************************/
1524 double ShhherCommand::getLikelihood(){
1528 vector<long double> P(numSeqs, 0);
1531 for(int i=0;i<numOTUs;i++){
1532 if(weight[i] > MIN_WEIGHT){
1538 for(int i=0;i<numOTUs;i++){
1540 if (m->control_pressed) { break; }
1542 for(int j=0;j<nSeqsPerOTU[i];j++){
1543 int index = cumNumSeqs[i] + j;
1544 int nI = seqIndex[index];
1545 double singleDist = dist[seqNumber[index]];
1547 P[nI] += weight[i] * exp(-singleDist * sigma);
1551 for(int i=0;i<numSeqs;i++){
1552 if(P[i] == 0){ P[i] = DBL_EPSILON; }
1557 nLL = nLL -(double)numSeqs * log(sigma);
1561 catch(exception& e) {
1562 m->errorOut(e, "ShhherCommand", "getNewWeights");
1567 /**************************************************************************************************/
1569 void ShhherCommand::checkCentroids(){
1571 vector<int> unique(numOTUs, 1);
1573 for(int i=0;i<numOTUs;i++){
1574 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1579 for(int i=0;i<numOTUs;i++){
1581 if (m->control_pressed) { break; }
1584 for(int j=i+1;j<numOTUs;j++){
1587 if(centroids[j] == centroids[i]){
1591 weight[i] += weight[j];
1599 catch(exception& e) {
1600 m->errorOut(e, "ShhherCommand", "checkCentroids");
1604 /**************************************************************************************************/
1608 void ShhherCommand::writeQualities(vector<int> otuCounts){
1611 string thisOutputDir = outputDir;
1612 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1613 string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.qual";
1615 ofstream qualityFile;
1616 m->openOutputFile(qualityFileName, qualityFile);
1618 qualityFile.setf(ios::fixed, ios::floatfield);
1619 qualityFile.setf(ios::showpoint);
1620 qualityFile << setprecision(6);
1622 vector<vector<int> > qualities(numOTUs);
1623 vector<double> pr(HOMOPS, 0);
1626 for(int i=0;i<numOTUs;i++){
1628 if (m->control_pressed) { break; }
1633 if(nSeqsPerOTU[i] > 0){
1634 qualities[i].assign(1024, -1);
1636 while(index < numFlowCells){
1637 double maxPrValue = 1e8;
1638 short maxPrIndex = -1;
1639 double count = 0.0000;
1641 pr.assign(HOMOPS, 0);
1643 for(int j=0;j<nSeqsPerOTU[i];j++){
1644 int lIndex = cumNumSeqs[i] + j;
1645 double tauValue = singleTau[seqNumber[lIndex]];
1646 int sequenceIndex = aaI[i][j];
1647 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1651 for(int s=0;s<HOMOPS;s++){
1652 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1656 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1657 maxPrValue = pr[maxPrIndex];
1659 if(count > MIN_COUNT){
1661 double norm = 0.0000;
1663 for(int s=0;s<HOMOPS;s++){
1664 norm += exp(-(pr[s] - maxPrValue));
1667 for(int s=1;s<=maxPrIndex;s++){
1669 double temp = 0.0000;
1671 U += exp(-(pr[s-1]-maxPrValue))/norm;
1679 temp = floor(-10 * temp);
1680 value = (int)floor(temp);
1681 if(value > 100){ value = 100; }
1683 qualities[i][base] = (int)value;
1693 if(otuCounts[i] > 0){
1694 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1696 int j=4; //need to get past the first four bases
1697 while(qualities[i][j] != -1){
1698 qualityFile << qualities[i][j] << ' ';
1701 qualityFile << endl;
1704 qualityFile.close();
1705 outputNames.push_back(qualityFileName);
1708 catch(exception& e) {
1709 m->errorOut(e, "ShhherCommand", "writeQualities");
1714 /**************************************************************************************************/
1716 void ShhherCommand::writeSequences(vector<int> otuCounts){
1718 string thisOutputDir = outputDir;
1719 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1720 string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.fasta";
1722 m->openOutputFile(fastaFileName, fastaFile);
1724 vector<string> names(numOTUs, "");
1726 for(int i=0;i<numOTUs;i++){
1728 if (m->control_pressed) { break; }
1730 int index = centroids[i];
1732 if(otuCounts[i] > 0){
1733 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
1737 for(int j=0;j<numFlowCells;j++){
1739 char base = flowOrder[j % 4];
1740 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
1745 fastaFile << newSeq.substr(4) << endl;
1750 outputNames.push_back(fastaFileName);
1752 if(compositeFASTAFileName != ""){
1753 m->appendFiles(fastaFileName, compositeFASTAFileName);
1756 catch(exception& e) {
1757 m->errorOut(e, "ShhherCommand", "writeSequences");
1762 /**************************************************************************************************/
1764 void ShhherCommand::writeNames(vector<int> otuCounts){
1766 string thisOutputDir = outputDir;
1767 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1768 string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.names";
1770 m->openOutputFile(nameFileName, nameFile);
1772 for(int i=0;i<numOTUs;i++){
1774 if (m->control_pressed) { break; }
1776 if(otuCounts[i] > 0){
1777 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
1779 for(int j=1;j<nSeqsPerOTU[i];j++){
1780 nameFile << ',' << seqNameVector[aaI[i][j]];
1787 outputNames.push_back(nameFileName);
1790 if(compositeNamesFileName != ""){
1791 m->appendFiles(nameFileName, compositeNamesFileName);
1794 catch(exception& e) {
1795 m->errorOut(e, "ShhherCommand", "writeNames");
1800 /**************************************************************************************************/
1802 void ShhherCommand::writeGroups(){
1804 string thisOutputDir = outputDir;
1805 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1806 string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1807 string groupFileName = fileRoot + "shhh.groups";
1809 m->openOutputFile(groupFileName, groupFile);
1811 for(int i=0;i<numSeqs;i++){
1812 if (m->control_pressed) { break; }
1813 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
1816 outputNames.push_back(groupFileName);
1819 catch(exception& e) {
1820 m->errorOut(e, "ShhherCommand", "writeGroups");
1825 /**************************************************************************************************/
1827 void ShhherCommand::writeClusters(vector<int> otuCounts){
1829 string thisOutputDir = outputDir;
1830 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1831 string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.counts";
1832 ofstream otuCountsFile;
1833 m->openOutputFile(otuCountsFileName, otuCountsFile);
1835 string bases = flowOrder;
1837 for(int i=0;i<numOTUs;i++){
1839 if (m->control_pressed) {
1842 //output the translated version of the centroid sequence for the otu
1843 if(otuCounts[i] > 0){
1844 int index = centroids[i];
1846 otuCountsFile << "ideal\t";
1847 for(int j=8;j<numFlowCells;j++){
1848 char base = bases[j % 4];
1849 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
1850 otuCountsFile << base;
1853 otuCountsFile << endl;
1855 for(int j=0;j<nSeqsPerOTU[i];j++){
1856 int sequence = aaI[i][j];
1857 otuCountsFile << seqNameVector[sequence] << '\t';
1861 for(int k=0;k<lengths[sequence];k++){
1862 char base = bases[k % 4];
1863 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
1865 for(int s=0;s<freq;s++){
1867 //otuCountsFile << base;
1870 otuCountsFile << newSeq.substr(4) << endl;
1872 otuCountsFile << endl;
1875 otuCountsFile.close();
1876 outputNames.push_back(otuCountsFileName);
1879 catch(exception& e) {
1880 m->errorOut(e, "ShhherCommand", "writeClusters");
1886 //**********************************************************************************************************************
1888 int ShhherCommand::execute(){
1890 if (abort == true) { return 0; }
1892 getSingleLookUp(); if (m->control_pressed) { return 0; }
1893 getJointLookUp(); if (m->control_pressed) { return 0; }
1895 int numFiles = flowFileVector.size();
1897 if (numFiles < processors) { processors = numFiles; }
1899 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1900 if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size()); }
1901 else { createProcesses(flowFileVector); } //each processor processes one file
1903 driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size());
1906 if(compositeFASTAFileName != ""){
1907 outputNames.push_back(compositeFASTAFileName);
1908 outputNames.push_back(compositeNamesFileName);
1911 m->mothurOutEndLine();
1912 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
1913 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
1914 m->mothurOutEndLine();
1918 catch(exception& e) {
1919 m->errorOut(e, "ShhherCommand", "execute");
1924 /**************************************************************************************************/
1926 int ShhherCommand::createProcesses(vector<string> filenames){
1928 vector<int> processIDS;
1933 if (filenames.size() < processors) { processors = filenames.size(); }
1935 //divide the groups between the processors
1936 vector<linePair> lines;
1937 vector<int> numFilesToComplete;
1938 int numFilesPerProcessor = filenames.size() / processors;
1939 for (int i = 0; i < processors; i++) {
1940 int startIndex = i * numFilesPerProcessor;
1941 int endIndex = (i+1) * numFilesPerProcessor;
1942 if(i == (processors - 1)){ endIndex = filenames.size(); }
1943 lines.push_back(linePair(startIndex, endIndex));
1944 numFilesToComplete.push_back((endIndex-startIndex));
1947 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1949 //loop through and create all the processes you want
1950 while (process != processors) {
1954 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1956 }else if (pid == 0){
1957 num = driver(filenames, compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp", lines[process].start, lines[process].end);
1959 //pass numSeqs to parent
1961 string tempFile = compositeFASTAFileName + toString(getpid()) + ".num.temp";
1962 m->openOutputFile(tempFile, out);
1968 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1969 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1975 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[0].start, lines[0].end);
1977 //force parent to wait until all the processes are done
1978 for (int i=0;i<processIDS.size();i++) {
1979 int temp = processIDS[i];
1985 //////////////////////////////////////////////////////////////////////////////////////////////////////
1987 /////////////////////// NOT WORKING, ACCESS VIOLATION ON READ OF FLOWGRAMS IN THREAD /////////////////
1989 //////////////////////////////////////////////////////////////////////////////////////////////////////
1990 //Windows version shared memory, so be careful when passing variables through the shhhFlowsData struct.
1991 //Above fork() will clone, so memory is separate, but that's not the case with windows,
1992 //////////////////////////////////////////////////////////////////////////////////////////////////////
1994 vector<shhhFlowsData*> pDataArray;
1995 DWORD dwThreadIdArray[processors-1];
1996 HANDLE hThreadArray[processors-1];
1998 //Create processor worker threads.
1999 for( int i=0; i<processors-1; i++ ){
2000 // Allocate memory for thread data.
2001 string extension = "";
2002 if (i != 0) { extension = toString(i) + ".temp"; }
2004 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);
2005 pDataArray.push_back(tempFlow);
2006 processIDS.push_back(i);
2008 hThreadArray[i] = CreateThread(NULL, 0, ShhhFlowsThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
2011 //using the main process as a worker saves time and memory
2013 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[processors-1].start, lines[processors-1].end);
2015 //Wait until all threads have terminated.
2016 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
2018 //Close all thread handles and free memory allocations.
2019 for(int i=0; i < pDataArray.size(); i++){
2020 for(int j=0; j < pDataArray[i]->outputNames.size(); j++){ outputNames.push_back(pDataArray[i]->outputNames[j]); }
2021 CloseHandle(hThreadArray[i]);
2022 delete pDataArray[i];
2027 for (int i=0;i<processIDS.size();i++) {
2029 string tempFile = compositeFASTAFileName + toString(processIDS[i]) + ".num.temp";
2030 m->openInputFile(tempFile, in);
2034 if (tempNum != numFilesToComplete[i+1]) {
2035 m->mothurOut("[ERROR]: main process expected " + toString(processIDS[i]) + " to complete " + toString(numFilesToComplete[i+1]) + " files, and it only reported completing " + toString(tempNum) + ". This will cause file mismatches. The flow files may be too large to process with multiple processors. \n");
2038 in.close(); m->mothurRemove(tempFile);
2040 if (compositeFASTAFileName != "") {
2041 m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName);
2042 m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName);
2043 m->mothurRemove((compositeFASTAFileName + toString(processIDS[i]) + ".temp"));
2044 m->mothurRemove((compositeNamesFileName + toString(processIDS[i]) + ".temp"));
2051 catch(exception& e) {
2052 m->errorOut(e, "ShhherCommand", "createProcesses");
2056 /**************************************************************************************************/
2058 vector<string> ShhherCommand::parseFlowFiles(string filename){
2060 vector<string> files;
2064 m->openInputFile(filename, in);
2066 int thisNumFLows = 0;
2067 in >> thisNumFLows; m->gobble(in);
2070 if (m->control_pressed) { break; }
2073 string outputFileName = filename + toString(count) + ".temp";
2074 m->openOutputFile(outputFileName, out);
2075 out << thisNumFLows << endl;
2076 files.push_back(outputFileName);
2078 int numLinesWrote = 0;
2079 for (int i = 0; i < largeSize; i++) {
2080 if (in.eof()) { break; }
2081 string line = m->getline(in); m->gobble(in);
2082 out << line << endl;
2087 if (numLinesWrote == 0) { m->mothurRemove(outputFileName); files.pop_back(); }
2092 if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } files.clear(); }
2094 m->mothurOut("\nDivided " + filename + " into " + toString(files.size()) + " files.\n\n");
2098 catch(exception& e) {
2099 m->errorOut(e, "ShhherCommand", "parseFlowFiles");
2103 /**************************************************************************************************/
2105 int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){
2108 int numCompleted = 0;
2110 for(int i=start;i<end;i++){
2112 if (m->control_pressed) { break; }
2114 vector<string> theseFlowFileNames; theseFlowFileNames.push_back(filenames[i]);
2115 if (large) { theseFlowFileNames = parseFlowFiles(filenames[i]); }
2117 if (m->control_pressed) { break; }
2119 double begClock = clock();
2120 unsigned long long begTime;
2122 for (int g = 0; g < theseFlowFileNames.size(); g++) {
2124 string flowFileName = theseFlowFileNames[g];
2125 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n");
2126 m->mothurOut("Reading flowgrams...\n");
2128 vector<string> seqNameVector;
2129 vector<int> lengths;
2130 vector<short> flowDataIntI;
2131 vector<double> flowDataPrI;
2132 map<string, int> nameMap;
2133 vector<short> uniqueFlowgrams;
2134 vector<int> uniqueCount;
2135 vector<int> mapSeqToUnique;
2136 vector<int> mapUniqueToSeq;
2137 vector<int> uniqueLengths;
2140 int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells);
2142 if (m->control_pressed) { break; }
2144 m->mothurOut("Identifying unique flowgrams...\n");
2145 int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI);
2147 if (m->control_pressed) { break; }
2149 m->mothurOut("Calculating distances between flowgrams...\n");
2150 string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
2151 begTime = time(NULL);
2154 flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2156 m->mothurOutEndLine();
2157 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
2160 string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
2161 createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq);
2163 if (m->control_pressed) { break; }
2165 m->mothurOut("\nClustering flowgrams...\n");
2166 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
2167 cluster(listFileName, distFileName, namesFileName);
2169 if (m->control_pressed) { break; }
2171 vector<int> otuData;
2172 vector<int> cumNumSeqs;
2173 vector<int> nSeqsPerOTU;
2174 vector<vector<int> > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2175 vector<vector<int> > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2176 vector<int> seqNumber; //tMaster->anP: the sequence id number sorted by OTU
2177 vector<int> seqIndex; //tMaster->anI; the index that corresponds to seqNumber
2180 int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap);
2182 if (m->control_pressed) { break; }
2184 m->mothurRemove(distFileName);
2185 m->mothurRemove(namesFileName);
2186 m->mothurRemove(listFileName);
2188 vector<double> dist; //adDist - distance of sequences to centroids
2189 vector<short> change; //did the centroid sequence change? 0 = no; 1 = yes
2190 vector<int> centroids; //the representative flowgram for each cluster m
2191 vector<double> weight;
2192 vector<double> singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2193 vector<int> nSeqsBreaks;
2194 vector<int> nOTUsBreaks;
2196 dist.assign(numSeqs * numOTUs, 0);
2197 change.assign(numOTUs, 1);
2198 centroids.assign(numOTUs, -1);
2199 weight.assign(numOTUs, 0);
2200 singleTau.assign(numSeqs, 1.0);
2202 nSeqsBreaks.assign(2, 0);
2203 nOTUsBreaks.assign(2, 0);
2206 nSeqsBreaks[1] = numSeqs;
2207 nOTUsBreaks[1] = numOTUs;
2209 if (m->control_pressed) { break; }
2211 double maxDelta = 0;
2215 begTime = time(NULL);
2217 m->mothurOut("\nDenoising flowgrams...\n");
2218 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
2220 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
2222 if (m->control_pressed) { break; }
2224 double cycClock = clock();
2225 unsigned long long cycTime = time(NULL);
2226 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2228 if (m->control_pressed) { break; }
2230 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2232 if (m->control_pressed) { break; }
2234 maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight);
2236 if (m->control_pressed) { break; }
2238 double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight);
2240 if (m->control_pressed) { break; }
2242 checkCentroids(numOTUs, centroids, weight);
2244 if (m->control_pressed) { break; }
2246 calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths);
2248 if (m->control_pressed) { break; }
2252 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
2256 if (m->control_pressed) { break; }
2258 m->mothurOut("\nFinalizing...\n");
2259 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2261 if (m->control_pressed) { break; }
2263 setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI);
2265 if (m->control_pressed) { break; }
2267 vector<int> otuCounts(numOTUs, 0);
2268 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
2270 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2272 if (m->control_pressed) { break; }
2274 if ((large) && (g == 0)) { flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; }
2275 string thisOutputDir = outputDir;
2276 if (outputDir == "") { thisOutputDir = m->hasPath(flowFileName); }
2277 string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.qual";
2278 string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.fasta";
2279 string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.names";
2280 string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.counts";
2281 string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
2282 string groupFileName = fileRoot + "shhh.groups";
2285 writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; }
2286 writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; }
2287 writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; }
2288 writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; }
2289 writeGroups(groupFileName, fileRoot, numSeqs, seqNameVector); if (m->control_pressed) { break; }
2293 m->appendFiles(qualityFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.qual"));
2294 m->mothurRemove(qualityFileName);
2295 m->appendFiles(fastaFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.fasta"));
2296 m->mothurRemove(fastaFileName);
2297 m->appendFiles(nameFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.names"));
2298 m->mothurRemove(nameFileName);
2299 m->appendFiles(otuCountsFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.counts"));
2300 m->mothurRemove(otuCountsFileName);
2301 m->appendFiles(groupFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.groups"));
2302 m->mothurRemove(groupFileName);
2304 m->mothurRemove(theseFlowFileNames[g]);
2309 m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
2312 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
2314 return numCompleted;
2316 }catch(exception& e) {
2317 m->errorOut(e, "ShhherCommand", "driver");
2322 /**************************************************************************************************/
2323 int ShhherCommand::getFlowData(string filename, vector<string>& thisSeqNameVector, vector<int>& thisLengths, vector<short>& thisFlowDataIntI, map<string, int>& thisNameMap, int& numFlowCells){
2328 m->openInputFile(filename, flowFile);
2331 int currentNumFlowCells;
2333 thisSeqNameVector.clear();
2334 thisLengths.clear();
2335 thisFlowDataIntI.clear();
2336 thisNameMap.clear();
2338 flowFile >> numFlowCells;
2339 int index = 0;//pcluster
2340 while(!flowFile.eof()){
2342 if (m->control_pressed) { break; }
2344 flowFile >> seqName >> currentNumFlowCells;
2345 thisLengths.push_back(currentNumFlowCells);
2347 thisSeqNameVector.push_back(seqName);
2348 thisNameMap[seqName] = index++;//pcluster
2350 for(int i=0;i<numFlowCells;i++){
2351 flowFile >> intensity;
2352 if(intensity > 9.99) { intensity = 9.99; }
2353 int intI = int(100 * intensity + 0.0001);
2354 thisFlowDataIntI.push_back(intI);
2356 m->gobble(flowFile);
2360 int numSeqs = thisSeqNameVector.size();
2362 for(int i=0;i<numSeqs;i++){
2364 if (m->control_pressed) { break; }
2366 int iNumFlowCells = i * numFlowCells;
2367 for(int j=thisLengths[i];j<numFlowCells;j++){
2368 thisFlowDataIntI[iNumFlowCells + j] = 0;
2375 catch(exception& e) {
2376 m->errorOut(e, "ShhherCommand", "getFlowData");
2380 /**************************************************************************************************/
2382 int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int stopSeq, vector<int>& mapUniqueToSeq, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2385 ostringstream outStream;
2386 outStream.setf(ios::fixed, ios::floatfield);
2387 outStream.setf(ios::dec, ios::basefield);
2388 outStream.setf(ios::showpoint);
2389 outStream.precision(6);
2391 int begTime = time(NULL);
2392 double begClock = clock();
2394 for(int i=0;i<stopSeq;i++){
2396 if (m->control_pressed) { break; }
2398 for(int j=0;j<i;j++){
2399 float flowDistance = calcPairwiseDist(numFlowCells, mapUniqueToSeq[i], mapUniqueToSeq[j], mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2401 if(flowDistance < 1e-6){
2402 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
2404 else if(flowDistance <= cutoff){
2405 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
2409 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
2410 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2411 m->mothurOutEndLine();
2415 ofstream distFile(distFileName.c_str());
2416 distFile << outStream.str();
2419 if (m->control_pressed) {}
2421 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
2422 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2423 m->mothurOutEndLine();
2428 catch(exception& e) {
2429 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
2433 /**************************************************************************************************/
2435 float ShhherCommand::calcPairwiseDist(int numFlowCells, int seqA, int seqB, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2437 int minLength = lengths[mapSeqToUnique[seqA]];
2438 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
2440 int ANumFlowCells = seqA * numFlowCells;
2441 int BNumFlowCells = seqB * numFlowCells;
2445 for(int i=0;i<minLength;i++){
2447 if (m->control_pressed) { break; }
2449 int flowAIntI = flowDataIntI[ANumFlowCells + i];
2450 float flowAPrI = flowDataPrI[ANumFlowCells + i];
2452 int flowBIntI = flowDataIntI[BNumFlowCells + i];
2453 float flowBPrI = flowDataPrI[BNumFlowCells + i];
2454 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
2457 dist /= (float) minLength;
2460 catch(exception& e) {
2461 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
2466 /**************************************************************************************************/
2468 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){
2471 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
2472 uniqueCount.assign(numSeqs, 0); // anWeights
2473 uniqueLengths.assign(numSeqs, 0);
2474 mapSeqToUnique.assign(numSeqs, -1);
2475 mapUniqueToSeq.assign(numSeqs, -1);
2477 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
2479 for(int i=0;i<numSeqs;i++){
2481 if (m->control_pressed) { break; }
2485 vector<short> current(numFlowCells);
2486 for(int j=0;j<numFlowCells;j++){
2487 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
2490 for(int j=0;j<numUniques;j++){
2491 int offset = j * numFlowCells;
2495 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
2496 else { shorterLength = uniqueLengths[j]; }
2498 for(int k=0;k<shorterLength;k++){
2499 if(current[k] != uniqueFlowgrams[offset + k]){
2506 mapSeqToUnique[i] = j;
2509 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
2515 if(index == numUniques){
2516 uniqueLengths[numUniques] = lengths[i];
2517 uniqueCount[numUniques] = 1;
2518 mapSeqToUnique[i] = numUniques;//anMap
2519 mapUniqueToSeq[numUniques] = i;//anF
2521 for(int k=0;k<numFlowCells;k++){
2522 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
2523 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
2529 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
2530 uniqueLengths.resize(numUniques);
2532 flowDataPrI.resize(numSeqs * numFlowCells, 0);
2533 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
2537 catch(exception& e) {
2538 m->errorOut(e, "ShhherCommand", "getUniques");
2542 /**************************************************************************************************/
2543 int ShhherCommand::createNamesFile(int numSeqs, int numUniques, string filename, vector<string>& seqNameVector, vector<int>& mapSeqToUnique, vector<int>& mapUniqueToSeq){
2546 vector<string> duplicateNames(numUniques, "");
2547 for(int i=0;i<numSeqs;i++){
2548 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
2552 m->openOutputFile(filename, nameFile);
2554 for(int i=0;i<numUniques;i++){
2556 if (m->control_pressed) { break; }
2558 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2559 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2566 catch(exception& e) {
2567 m->errorOut(e, "ShhherCommand", "createNamesFile");
2571 //**********************************************************************************************************************
2573 int ShhherCommand::cluster(string filename, string distFileName, string namesFileName){
2576 ReadMatrix* read = new ReadColumnMatrix(distFileName);
2577 read->setCutoff(cutoff);
2579 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
2580 clusterNameMap->readMap();
2581 read->read(clusterNameMap);
2583 ListVector* list = read->getListVector();
2584 SparseMatrix* matrix = read->getMatrix();
2587 delete clusterNameMap;
2589 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
2591 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
2592 string tag = cluster->getTag();
2594 double clusterCutoff = cutoff;
2595 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
2597 if (m->control_pressed) { break; }
2599 cluster->update(clusterCutoff);
2602 list->setLabel(toString(cutoff));
2605 m->openOutputFile(filename, listFile);
2606 list->print(listFile);
2609 delete matrix; delete cluster; delete rabund; delete list;
2613 catch(exception& e) {
2614 m->errorOut(e, "ShhherCommand", "cluster");
2618 /**************************************************************************************************/
2620 int ShhherCommand::getOTUData(int numSeqs, string fileName, vector<int>& otuData,
2621 vector<int>& cumNumSeqs,
2622 vector<int>& nSeqsPerOTU,
2623 vector<vector<int> >& aaP, //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2624 vector<vector<int> >& aaI, //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2625 vector<int>& seqNumber, //tMaster->anP: the sequence id number sorted by OTU
2626 vector<int>& seqIndex,
2627 map<string, int>& nameMap){
2631 m->openInputFile(fileName, listFile);
2635 listFile >> label >> numOTUs;
2637 otuData.assign(numSeqs, 0);
2638 cumNumSeqs.assign(numOTUs, 0);
2639 nSeqsPerOTU.assign(numOTUs, 0);
2640 aaP.clear();aaP.resize(numOTUs);
2646 string singleOTU = "";
2648 for(int i=0;i<numOTUs;i++){
2650 if (m->control_pressed) { break; }
2652 listFile >> singleOTU;
2654 istringstream otuString(singleOTU);
2658 string seqName = "";
2660 for(int j=0;j<singleOTU.length();j++){
2661 char letter = otuString.get();
2667 map<string,int>::iterator nmIt = nameMap.find(seqName);
2668 int index = nmIt->second;
2670 nameMap.erase(nmIt);
2674 aaP[i].push_back(index);
2679 map<string,int>::iterator nmIt = nameMap.find(seqName);
2681 int index = nmIt->second;
2682 nameMap.erase(nmIt);
2686 aaP[i].push_back(index);
2691 sort(aaP[i].begin(), aaP[i].end());
2692 for(int j=0;j<nSeqsPerOTU[i];j++){
2693 seqNumber.push_back(aaP[i][j]);
2695 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
2696 aaP[i].push_back(0);
2702 for(int i=1;i<numOTUs;i++){
2703 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
2706 seqIndex = seqNumber;
2713 catch(exception& e) {
2714 m->errorOut(e, "ShhherCommand", "getOTUData");
2718 /**************************************************************************************************/
2720 int ShhherCommand::calcCentroidsDriver(int numOTUs,
2721 vector<int>& cumNumSeqs,
2722 vector<int>& nSeqsPerOTU,
2723 vector<int>& seqIndex,
2724 vector<short>& change, //did the centroid sequence change? 0 = no; 1 = yes
2725 vector<int>& centroids, //the representative flowgram for each cluster m
2726 vector<double>& singleTau, //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2727 vector<int>& mapSeqToUnique,
2728 vector<short>& uniqueFlowgrams,
2729 vector<short>& flowDataIntI,
2730 vector<int>& lengths,
2732 vector<int>& seqNumber){
2734 //this function gets the most likely homopolymer length at a flow position for a group of sequences
2739 for(int i=0;i<numOTUs;i++){
2741 if (m->control_pressed) { break; }
2745 int minFlowGram = 100000000;
2746 double minFlowValue = 1e8;
2747 change[i] = 0; //FALSE
2749 for(int j=0;j<nSeqsPerOTU[i];j++){
2750 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
2753 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
2754 vector<double> adF(nSeqsPerOTU[i]);
2755 vector<int> anL(nSeqsPerOTU[i]);
2757 for(int j=0;j<nSeqsPerOTU[i];j++){
2758 int index = cumNumSeqs[i] + j;
2759 int nI = seqIndex[index];
2760 int nIU = mapSeqToUnique[nI];
2763 for(k=0;k<position;k++){
2769 anL[position] = nIU;
2770 adF[position] = 0.0000;
2775 for(int j=0;j<nSeqsPerOTU[i];j++){
2776 int index = cumNumSeqs[i] + j;
2777 int nI = seqIndex[index];
2779 double tauValue = singleTau[seqNumber[index]];
2781 for(int k=0;k<position;k++){
2782 double dist = getDistToCentroid(anL[k], nI, lengths[nI], uniqueFlowgrams, flowDataIntI, numFlowCells);
2783 adF[k] += dist * tauValue;
2787 for(int j=0;j<position;j++){
2788 if(adF[j] < minFlowValue){
2790 minFlowValue = adF[j];
2794 if(centroids[i] != anL[minFlowGram]){
2796 centroids[i] = anL[minFlowGram];
2799 else if(centroids[i] != -1){
2807 catch(exception& e) {
2808 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
2812 /**************************************************************************************************/
2814 double ShhherCommand::getDistToCentroid(int cent, int flow, int length, vector<short>& uniqueFlowgrams,
2815 vector<short>& flowDataIntI, int numFlowCells){
2818 int flowAValue = cent * numFlowCells;
2819 int flowBValue = flow * numFlowCells;
2823 for(int i=0;i<length;i++){
2824 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
2829 return dist / (double)length;
2831 catch(exception& e) {
2832 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
2836 /**************************************************************************************************/
2838 double ShhherCommand::getNewWeights(int numOTUs, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU, vector<double>& singleTau, vector<int>& seqNumber, vector<double>& weight){
2841 double maxChange = 0;
2843 for(int i=0;i<numOTUs;i++){
2845 if (m->control_pressed) { break; }
2847 double difference = weight[i];
2850 for(int j=0;j<nSeqsPerOTU[i];j++){
2851 int index = cumNumSeqs[i] + j;
2852 double tauValue = singleTau[seqNumber[index]];
2853 weight[i] += tauValue;
2856 difference = fabs(weight[i] - difference);
2857 if(difference > maxChange){ maxChange = difference; }
2861 catch(exception& e) {
2862 m->errorOut(e, "ShhherCommand", "getNewWeights");
2867 /**************************************************************************************************/
2869 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){
2873 vector<long double> P(numSeqs, 0);
2876 for(int i=0;i<numOTUs;i++){
2877 if(weight[i] > MIN_WEIGHT){
2883 for(int i=0;i<numOTUs;i++){
2885 if (m->control_pressed) { break; }
2887 for(int j=0;j<nSeqsPerOTU[i];j++){
2888 int index = cumNumSeqs[i] + j;
2889 int nI = seqIndex[index];
2890 double singleDist = dist[seqNumber[index]];
2892 P[nI] += weight[i] * exp(-singleDist * sigma);
2896 for(int i=0;i<numSeqs;i++){
2897 if(P[i] == 0){ P[i] = DBL_EPSILON; }
2902 nLL = nLL -(double)numSeqs * log(sigma);
2906 catch(exception& e) {
2907 m->errorOut(e, "ShhherCommand", "getNewWeights");
2912 /**************************************************************************************************/
2914 int ShhherCommand::checkCentroids(int numOTUs, vector<int>& centroids, vector<double>& weight){
2916 vector<int> unique(numOTUs, 1);
2918 for(int i=0;i<numOTUs;i++){
2919 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
2924 for(int i=0;i<numOTUs;i++){
2926 if (m->control_pressed) { break; }
2929 for(int j=i+1;j<numOTUs;j++){
2932 if(centroids[j] == centroids[i]){
2936 weight[i] += weight[j];
2946 catch(exception& e) {
2947 m->errorOut(e, "ShhherCommand", "checkCentroids");
2951 /**************************************************************************************************/
2953 void ShhherCommand::calcNewDistances(int numSeqs, int numOTUs, vector<int>& nSeqsPerOTU, vector<double>& dist,
2954 vector<double>& weight, vector<short>& change, vector<int>& centroids,
2955 vector<vector<int> >& aaP, vector<double>& singleTau, vector<vector<int> >& aaI,
2956 vector<int>& seqNumber, vector<int>& seqIndex,
2957 vector<short>& uniqueFlowgrams,
2958 vector<short>& flowDataIntI, int numFlowCells, vector<int>& lengths){
2963 vector<double> newTau(numOTUs,0);
2964 vector<double> norms(numSeqs, 0);
2965 nSeqsPerOTU.assign(numOTUs, 0);
2967 for(int i=0;i<numSeqs;i++){
2969 if (m->control_pressed) { break; }
2971 int indexOffset = i * numOTUs;
2973 double offset = 1e8;
2975 for(int j=0;j<numOTUs;j++){
2977 if(weight[j] > MIN_WEIGHT && change[j] == 1){
2978 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i], uniqueFlowgrams, flowDataIntI, numFlowCells);
2981 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
2982 offset = dist[indexOffset + j];
2986 for(int j=0;j<numOTUs;j++){
2987 if(weight[j] > MIN_WEIGHT){
2988 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
2989 norms[i] += newTau[j];
2996 for(int j=0;j<numOTUs;j++){
2997 newTau[j] /= norms[i];
3000 for(int j=0;j<numOTUs;j++){
3001 if(newTau[j] > MIN_TAU){
3003 int oldTotal = total;
3007 singleTau.resize(total, 0);
3008 seqNumber.resize(total, 0);
3009 seqIndex.resize(total, 0);
3011 singleTau[oldTotal] = newTau[j];
3013 aaP[j][nSeqsPerOTU[j]] = oldTotal;
3014 aaI[j][nSeqsPerOTU[j]] = i;
3022 catch(exception& e) {
3023 m->errorOut(e, "ShhherCommand", "calcNewDistances");
3027 /**************************************************************************************************/
3029 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){
3032 for(int i=0;i<numOTUs;i++){
3034 if (m->control_pressed) { return 0; }
3036 cumNumSeqs[i] = index;
3037 for(int j=0;j<nSeqsPerOTU[i];j++){
3038 seqNumber[index] = aaP[i][j];
3039 seqIndex[index] = aaI[i][j];
3047 catch(exception& e) {
3048 m->errorOut(e, "ShhherCommand", "fill");
3052 /**************************************************************************************************/
3054 void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector<int>& seqNumber, vector<int>& seqIndex, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU,
3055 vector<int>& otuData, vector<double>& singleTau, vector<double>& dist, vector<vector<int> >& aaP, vector<vector<int> >& aaI){
3058 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
3060 for(int i=0;i<numOTUs;i++){
3062 if (m->control_pressed) { break; }
3064 for(int j=0;j<nSeqsPerOTU[i];j++){
3065 int index = cumNumSeqs[i] + j;
3066 double tauValue = singleTau[seqNumber[index]];
3067 int sIndex = seqIndex[index];
3068 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
3072 for(int i=0;i<numSeqs;i++){
3073 double maxTau = -1.0000;
3075 for(int j=0;j<numOTUs;j++){
3076 if(bigTauMatrix[i * numOTUs + j] > maxTau){
3077 maxTau = bigTauMatrix[i * numOTUs + j];
3082 otuData[i] = maxOTU;
3085 nSeqsPerOTU.assign(numOTUs, 0);
3087 for(int i=0;i<numSeqs;i++){
3088 int index = otuData[i];
3090 singleTau[i] = 1.0000;
3093 aaP[index][nSeqsPerOTU[index]] = i;
3094 aaI[index][nSeqsPerOTU[index]] = i;
3096 nSeqsPerOTU[index]++;
3099 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
3101 catch(exception& e) {
3102 m->errorOut(e, "ShhherCommand", "setOTUs");
3106 /**************************************************************************************************/
3108 void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string qualityFileName, vector<int> otuCounts, vector<int>& nSeqsPerOTU, vector<int>& seqNumber,
3109 vector<double>& singleTau, vector<short>& flowDataIntI, vector<short>& uniqueFlowgrams, vector<int>& cumNumSeqs,
3110 vector<int>& mapUniqueToSeq, vector<string>& seqNameVector, vector<int>& centroids, vector<vector<int> >& aaI){
3114 ofstream qualityFile;
3115 m->openOutputFile(qualityFileName, qualityFile);
3117 qualityFile.setf(ios::fixed, ios::floatfield);
3118 qualityFile.setf(ios::showpoint);
3119 qualityFile << setprecision(6);
3121 vector<vector<int> > qualities(numOTUs);
3122 vector<double> pr(HOMOPS, 0);
3125 for(int i=0;i<numOTUs;i++){
3127 if (m->control_pressed) { break; }
3132 if(nSeqsPerOTU[i] > 0){
3133 qualities[i].assign(1024, -1);
3135 while(index < numFlowCells){
3136 double maxPrValue = 1e8;
3137 short maxPrIndex = -1;
3138 double count = 0.0000;
3140 pr.assign(HOMOPS, 0);
3142 for(int j=0;j<nSeqsPerOTU[i];j++){
3143 int lIndex = cumNumSeqs[i] + j;
3144 double tauValue = singleTau[seqNumber[lIndex]];
3145 int sequenceIndex = aaI[i][j];
3146 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
3150 for(int s=0;s<HOMOPS;s++){
3151 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
3155 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
3156 maxPrValue = pr[maxPrIndex];
3158 if(count > MIN_COUNT){
3160 double norm = 0.0000;
3162 for(int s=0;s<HOMOPS;s++){
3163 norm += exp(-(pr[s] - maxPrValue));
3166 for(int s=1;s<=maxPrIndex;s++){
3168 double temp = 0.0000;
3170 U += exp(-(pr[s-1]-maxPrValue))/norm;
3178 temp = floor(-10 * temp);
3179 value = (int)floor(temp);
3180 if(value > 100){ value = 100; }
3182 qualities[i][base] = (int)value;
3192 if(otuCounts[i] > 0){
3193 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
3195 int j=4; //need to get past the first four bases
3196 while(qualities[i][j] != -1){
3197 qualityFile << qualities[i][j] << ' ';
3198 if (j > qualities[i].size()) { break; }
3201 qualityFile << endl;
3204 qualityFile.close();
3205 outputNames.push_back(qualityFileName);
3208 catch(exception& e) {
3209 m->errorOut(e, "ShhherCommand", "writeQualities");
3214 /**************************************************************************************************/
3216 void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTUs, int numFlowCells, string fastaFileName, vector<int> otuCounts, vector<short>& uniqueFlowgrams, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& centroids){
3220 m->openOutputFile(fastaFileName, fastaFile);
3222 vector<string> names(numOTUs, "");
3224 for(int i=0;i<numOTUs;i++){
3226 if (m->control_pressed) { break; }
3228 int index = centroids[i];
3230 if(otuCounts[i] > 0){
3231 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
3235 for(int j=0;j<numFlowCells;j++){
3237 char base = flowOrder[j % 4];
3238 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
3243 if (newSeq.length() >= 4) { fastaFile << newSeq.substr(4) << endl; }
3244 else { fastaFile << "NNNN" << endl; }
3249 outputNames.push_back(fastaFileName);
3251 if(thisCompositeFASTAFileName != ""){
3252 m->appendFiles(fastaFileName, thisCompositeFASTAFileName);
3255 catch(exception& e) {
3256 m->errorOut(e, "ShhherCommand", "writeSequences");
3261 /**************************************************************************************************/
3263 void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string nameFileName, vector<int> otuCounts, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& nSeqsPerOTU){
3267 m->openOutputFile(nameFileName, nameFile);
3269 for(int i=0;i<numOTUs;i++){
3271 if (m->control_pressed) { break; }
3273 if(otuCounts[i] > 0){
3274 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
3276 for(int j=1;j<nSeqsPerOTU[i];j++){
3277 nameFile << ',' << seqNameVector[aaI[i][j]];
3284 outputNames.push_back(nameFileName);
3287 if(thisCompositeNamesFileName != ""){
3288 m->appendFiles(nameFileName, thisCompositeNamesFileName);
3291 catch(exception& e) {
3292 m->errorOut(e, "ShhherCommand", "writeNames");
3297 /**************************************************************************************************/
3299 void ShhherCommand::writeGroups(string groupFileName, string fileRoot, int numSeqs, vector<string>& seqNameVector){
3302 m->openOutputFile(groupFileName, groupFile);
3304 for(int i=0;i<numSeqs;i++){
3305 if (m->control_pressed) { break; }
3306 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
3309 outputNames.push_back(groupFileName);
3312 catch(exception& e) {
3313 m->errorOut(e, "ShhherCommand", "writeGroups");
3318 /**************************************************************************************************/
3320 void ShhherCommand::writeClusters(string otuCountsFileName, 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){
3322 ofstream otuCountsFile;
3323 m->openOutputFile(otuCountsFileName, otuCountsFile);
3325 string bases = flowOrder;
3327 for(int i=0;i<numOTUs;i++){
3329 if (m->control_pressed) {
3332 //output the translated version of the centroid sequence for the otu
3333 if(otuCounts[i] > 0){
3334 int index = centroids[i];
3336 otuCountsFile << "ideal\t";
3337 for(int j=8;j<numFlowCells;j++){
3338 char base = bases[j % 4];
3339 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
3340 otuCountsFile << base;
3343 otuCountsFile << endl;
3345 for(int j=0;j<nSeqsPerOTU[i];j++){
3346 int sequence = aaI[i][j];
3347 otuCountsFile << seqNameVector[sequence] << '\t';
3351 for(int k=0;k<lengths[sequence];k++){
3352 char base = bases[k % 4];
3353 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
3355 for(int s=0;s<freq;s++){
3357 //otuCountsFile << base;
3361 if (newSeq.length() >= 4) { otuCountsFile << newSeq.substr(4) << endl; }
3362 else { otuCountsFile << "NNNN" << endl; }
3364 otuCountsFile << endl;
3367 otuCountsFile.close();
3368 outputNames.push_back(otuCountsFileName);
3371 catch(exception& e) {
3372 m->errorOut(e, "ShhherCommand", "writeClusters");
3377 /**************************************************************************************************/
3379 void ShhherCommand::getSingleLookUp(){
3381 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
3382 singleLookUp.assign(HOMOPS * NUMBINS, 0);
3385 ifstream lookUpFile;
3386 m->openInputFile(lookupFileName, lookUpFile);
3388 for(int i=0;i<HOMOPS;i++){
3390 if (m->control_pressed) { break; }
3393 lookUpFile >> logFracFreq;
3395 for(int j=0;j<NUMBINS;j++) {
3396 lookUpFile >> singleLookUp[index];
3402 catch(exception& e) {
3403 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
3408 /**************************************************************************************************/
3410 void ShhherCommand::getJointLookUp(){
3413 // the most likely joint probability (-log) that two intenities have the same polymer length
3414 jointLookUp.resize(NUMBINS * NUMBINS, 0);
3416 for(int i=0;i<NUMBINS;i++){
3418 if (m->control_pressed) { break; }
3420 for(int j=0;j<NUMBINS;j++){
3422 double minSum = 100000000;
3424 for(int k=0;k<HOMOPS;k++){
3425 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
3427 if(sum < minSum) { minSum = sum; }
3429 jointLookUp[i * NUMBINS + j] = minSum;
3433 catch(exception& e) {
3434 m->errorOut(e, "ShhherCommand", "getJointLookUp");
3439 /**************************************************************************************************/
3441 double ShhherCommand::getProbIntensity(int intIntensity){
3444 double minNegLogProb = 100000000;
3447 for(int i=0;i<HOMOPS;i++){//loop signal strength
3449 if (m->control_pressed) { break; }
3451 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
3452 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
3455 return minNegLogProb;
3457 catch(exception& e) {
3458 m->errorOut(e, "ShhherCommand", "getProbIntensity");