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","fasta-name-group-counts-qfile",false,false,true); parameters.push_back(pflow);
16 CommandParameter pfile("file", "InputTypes", "", "", "none", "fileflow", "none","fasta-name-group-counts-qfile",false,false,true); parameters.push_back(pfile);
17 CommandParameter plookup("lookup", "InputTypes", "", "", "none", "none", "none","",false,false,true); 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,true); 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 //**********************************************************************************************************************
50 string ShhherCommand::getOutputPattern(string type) {
54 if (type == "fasta") { pattern = "[filename],shhh.fasta"; }
55 else if (type == "name") { pattern = "[filename],shhh.names"; }
56 else if (type == "group") { pattern = "[filename],shhh.groups"; }
57 else if (type == "counts") { pattern = "[filename],shhh.counts"; }
58 else if (type == "qfile") { pattern = "[filename],shhh.qual"; }
59 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
64 m->errorOut(e, "ShhherCommand", "getOutputPattern");
68 //**********************************************************************************************************************
70 ShhherCommand::ShhherCommand(){
72 abort = true; calledHelp = true;
75 //initialize outputTypes
76 vector<string> tempOutNames;
77 outputTypes["fasta"] = tempOutNames;
78 outputTypes["name"] = tempOutNames;
79 outputTypes["group"] = tempOutNames;
80 outputTypes["counts"] = tempOutNames;
81 outputTypes["qfile"] = tempOutNames;
85 m->errorOut(e, "ShhherCommand", "ShhherCommand");
90 //**********************************************************************************************************************
92 ShhherCommand::ShhherCommand(string option) {
96 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
97 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
101 abort = false; calledHelp = false;
103 //allow user to run help
104 if(option == "help") { help(); abort = true; calledHelp = true; }
105 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
108 vector<string> myArray = setParameters();
110 OptionParser parser(option);
111 map<string,string> parameters = parser.getParameters();
113 ValidParameters validParameter;
114 map<string,string>::iterator it;
116 //check to make sure all parameters are valid for command
117 for (it = parameters.begin(); it != parameters.end(); it++) {
118 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
121 //initialize outputTypes
122 vector<string> tempOutNames;
123 outputTypes["fasta"] = tempOutNames;
124 outputTypes["name"] = tempOutNames;
125 outputTypes["group"] = tempOutNames;
126 outputTypes["counts"] = tempOutNames;
127 outputTypes["qfile"] = tempOutNames;
130 //if the user changes the input directory command factory will send this info to us in the output parameter
131 string inputDir = validParameter.validFile(parameters, "inputdir", false);
132 if (inputDir == "not found"){ inputDir = ""; }
135 it = parameters.find("flow");
136 //user has given a template file
137 if(it != parameters.end()){
138 path = m->hasPath(it->second);
139 //if the user has not given a path then, add inputdir. else leave path alone.
140 if (path == "") { parameters["flow"] = inputDir + it->second; }
143 it = parameters.find("lookup");
144 //user has given a template file
145 if(it != parameters.end()){
146 path = m->hasPath(it->second);
147 //if the user has not given a path then, add inputdir. else leave path alone.
148 if (path == "") { parameters["lookup"] = inputDir + it->second; }
151 it = parameters.find("file");
152 //user has given a template file
153 if(it != parameters.end()){
154 path = m->hasPath(it->second);
155 //if the user has not given a path then, add inputdir. else leave path alone.
156 if (path == "") { parameters["file"] = inputDir + it->second; }
160 //if the user changes the output directory command factory will send this info to us in the output parameter
161 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
163 //check for required parameters
164 flowFileName = validParameter.validFile(parameters, "flow", true);
165 flowFilesFileName = validParameter.validFile(parameters, "file", true);
166 if (flowFileName == "not found" && flowFilesFileName == "not found") {
167 m->mothurOut("values for either flow or file must be provided for the shhh.flows command.");
168 m->mothurOutEndLine();
171 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
173 if(flowFileName != "not found"){
174 compositeFASTAFileName = "";
175 compositeNamesFileName = "";
180 string thisoutputDir = outputDir;
181 if (outputDir == "") { thisoutputDir = m->hasPath(flowFilesFileName); } //if user entered a file with a path then preserve it
183 //we want to rip off .files, and also .flow if its there
184 string fileroot = m->getRootName(m->getSimpleName(flowFilesFileName));
185 if (fileroot[fileroot.length()-1] == '.') { fileroot = fileroot.substr(0, fileroot.length()-1); } //rip off dot
186 string extension = m->getExtension(fileroot);
187 if (extension == ".flow") { fileroot = m->getRootName(fileroot); }
188 else { fileroot += "."; } //add back if needed
190 compositeFASTAFileName = thisoutputDir + fileroot + "shhh.fasta";
191 m->openOutputFile(compositeFASTAFileName, temp);
194 compositeNamesFileName = thisoutputDir + fileroot + "shhh.names";
195 m->openOutputFile(compositeNamesFileName, temp);
199 if(flowFilesFileName != "not found"){
202 ifstream flowFilesFile;
203 m->openInputFile(flowFilesFileName, flowFilesFile);
204 while(flowFilesFile){
205 fName = m->getline(flowFilesFile);
207 //test if file is valid
209 int ableToOpen = m->openInputFile(fName, in, "noerror");
211 if (ableToOpen == 1) {
212 if (inputDir != "") { //default path is set
213 string tryPath = inputDir + fName;
214 m->mothurOut("Unable to open " + fName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
216 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
222 if (ableToOpen == 1) {
223 if (m->getDefaultPath() != "") { //default path is set
224 string tryPath = m->getDefaultPath() + m->getSimpleName(fName);
225 m->mothurOut("Unable to open " + fName + ". Trying default " + tryPath); m->mothurOutEndLine();
227 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
233 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
234 if (ableToOpen == 1) {
235 string exepath = m->argv;
236 string tempPath = exepath;
237 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
238 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
240 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(fName);
241 m->mothurOut("Unable to open " + fName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
243 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
248 if (ableToOpen == 1) { m->mothurOut("Unable to open " + fName + ". Disregarding. "); m->mothurOutEndLine(); }
249 else { flowFileVector.push_back(fName); }
250 m->gobble(flowFilesFile);
252 flowFilesFile.close();
253 if (flowFileVector.size() == 0) { m->mothurOut("[ERROR]: no valid files."); m->mothurOutEndLine(); abort = true; }
256 if (outputDir == "") { outputDir = m->hasPath(flowFileName); }
257 flowFileVector.push_back(flowFileName);
260 //check for optional parameter and set defaults
261 // ...at some point should added some additional type checking...
263 temp = validParameter.validFile(parameters, "lookup", true);
264 if (temp == "not found") {
265 lookupFileName = "LookUp_Titanium.pat";
269 ableToOpen = m->openInputFile(lookupFileName, in, "noerror");
272 //if you can't open it, try input location
273 if (ableToOpen == 1) {
274 if (inputDir != "") { //default path is set
275 string tryPath = inputDir + lookupFileName;
276 m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
278 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
280 lookupFileName = tryPath;
284 //if you can't open it, try default location
285 if (ableToOpen == 1) {
286 if (m->getDefaultPath() != "") { //default path is set
287 string tryPath = m->getDefaultPath() + m->getSimpleName(lookupFileName);
288 m->mothurOut("Unable to open " + lookupFileName + ". Trying default " + tryPath); m->mothurOutEndLine();
290 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
292 lookupFileName = tryPath;
296 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
297 if (ableToOpen == 1) {
298 string exepath = m->argv;
299 string tempPath = exepath;
300 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
301 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
303 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
304 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
306 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
308 lookupFileName = tryPath;
311 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
313 else if(temp == "not open") {
315 lookupFileName = validParameter.validFile(parameters, "lookup", false);
317 //if you can't open it its not inputDir, try mothur excutable location
318 string exepath = m->argv;
319 string tempPath = exepath;
320 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
321 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
323 string tryPath = m->getFullPathName(exepath) + lookupFileName;
324 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
326 int ableToOpen = m->openInputFile(tryPath, in2, "noerror");
328 lookupFileName = tryPath;
330 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
331 }else { lookupFileName = temp; }
333 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
334 m->setProcessors(temp);
335 m->mothurConvert(temp, processors);
337 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
338 m->mothurConvert(temp, cutoff);
340 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
341 m->mothurConvert(temp, minDelta);
343 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
344 m->mothurConvert(temp, maxIters);
346 temp = validParameter.validFile(parameters, "large", false); if (temp == "not found"){ temp = "0"; }
347 m->mothurConvert(temp, largeSize);
348 if (largeSize != 0) { large = true; }
349 else { large = false; }
350 if (largeSize < 0) { m->mothurOut("The value of the large cannot be negative.\n"); }
353 if (large) { m->mothurOut("The large parameter is not available with the MPI-Enabled version.\n"); large=false; }
357 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
358 m->mothurConvert(temp, sigma);
360 flowOrder = validParameter.validFile(parameters, "order", false);
361 if (flowOrder == "not found"){ flowOrder = "TACG"; }
362 else if(flowOrder.length() != 4){
363 m->mothurOut("The value of the order option must be four bases long\n");
371 catch(exception& e) {
372 m->errorOut(e, "ShhherCommand", "ShhherCommand");
376 //**********************************************************************************************************************
378 int ShhherCommand::execute(){
380 if (abort == true) { if (calledHelp) { return 0; } return 2; }
387 for(int i=1;i<ncpus;i++){
388 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
390 if(abort == 1){ return 0; }
394 m->mothurOut("\nGetting preliminary data...\n");
395 getSingleLookUp(); if (m->control_pressed) { return 0; }
396 getJointLookUp(); if (m->control_pressed) { return 0; }
398 vector<string> flowFileVector;
399 if(flowFilesFileName != "not found"){
402 ifstream flowFilesFile;
403 m->openInputFile(flowFilesFileName, flowFilesFile);
404 while(flowFilesFile){
405 fName = m->getline(flowFilesFile);
406 flowFileVector.push_back(fName);
407 m->gobble(flowFilesFile);
411 flowFileVector.push_back(flowFileName);
414 int numFiles = flowFileVector.size();
416 for(int i=1;i<ncpus;i++){
417 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
420 for(int i=0;i<numFiles;i++){
422 if (m->control_pressed) { break; }
424 double begClock = clock();
425 unsigned long long begTime = time(NULL);
427 flowFileName = flowFileVector[i];
429 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
430 m->mothurOut("Reading flowgrams...\n");
433 if (m->control_pressed) { break; }
435 m->mothurOut("Identifying unique flowgrams...\n");
438 if (m->control_pressed) { break; }
440 m->mothurOut("Calculating distances between flowgrams...\n");
442 strcpy(fileName, flowFileName.c_str());
444 for(int i=1;i<ncpus;i++){
445 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
447 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
448 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
449 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
450 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
451 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
452 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
453 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
454 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
455 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
456 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
459 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
461 if (m->control_pressed) { break; }
464 for(int i=1;i<ncpus;i++){
465 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
467 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
468 m->mothurRemove((distFileName + ".temp." + toString(i)));
471 string namesFileName = createNamesFile();
473 if (m->control_pressed) { break; }
475 m->mothurOut("\nClustering flowgrams...\n");
476 string listFileName = cluster(distFileName, namesFileName);
478 if (m->control_pressed) { break; }
482 getOTUData(listFileName);
484 m->mothurRemove(distFileName);
485 m->mothurRemove(namesFileName);
486 m->mothurRemove(listFileName);
488 if (m->control_pressed) { break; }
492 if (m->control_pressed) { break; }
495 for(int i=1;i<ncpus;i++){
496 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
497 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
498 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
499 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
502 if (m->control_pressed) { break; }
507 int numOTUsOnCPU = numOTUs / ncpus;
508 int numSeqsOnCPU = numSeqs / ncpus;
509 m->mothurOut("\nDenoising flowgrams...\n");
510 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
512 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
514 double cycClock = clock();
515 unsigned long long cycTime = time(NULL);
518 if (m->control_pressed) { break; }
520 int total = singleTau.size();
521 for(int i=1;i<ncpus;i++){
522 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
523 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
524 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
526 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
527 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
528 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
529 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
530 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
533 calcCentroidsDriver(0, numOTUsOnCPU);
535 for(int i=1;i<ncpus;i++){
536 int otuStart = i * numOTUs / ncpus;
537 int otuStop = (i + 1) * numOTUs / ncpus;
539 vector<int> tempCentroids(numOTUs, 0);
540 vector<short> tempChange(numOTUs, 0);
542 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
543 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
545 for(int j=otuStart;j<otuStop;j++){
546 centroids[j] = tempCentroids[j];
547 change[j] = tempChange[j];
551 maxDelta = getNewWeights(); if (m->control_pressed) { break; }
552 double nLL = getLikelihood(); if (m->control_pressed) { break; }
553 checkCentroids(); if (m->control_pressed) { break; }
555 for(int i=1;i<ncpus;i++){
556 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
557 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
558 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
561 calcNewDistancesParent(0, numSeqsOnCPU);
563 total = singleTau.size();
565 for(int i=1;i<ncpus;i++){
567 int seqStart = i * numSeqs / ncpus;
568 int seqStop = (i + 1) * numSeqs / ncpus;
570 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
572 vector<int> childSeqIndex(childTotal, 0);
573 vector<double> childSingleTau(childTotal, 0);
574 vector<double> childDist(numSeqs * numOTUs, 0);
575 vector<int> otuIndex(childTotal, 0);
577 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
578 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
579 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
580 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
582 int oldTotal = total;
584 singleTau.resize(total, 0);
585 seqIndex.resize(total, 0);
586 seqNumber.resize(total, 0);
590 for(int j=oldTotal;j<total;j++){
591 int otuI = otuIndex[childIndex];
592 int seqI = childSeqIndex[childIndex];
594 singleTau[j] = childSingleTau[childIndex];
596 aaP[otuI][nSeqsPerOTU[otuI]] = j;
597 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
602 int index = seqStart * numOTUs;
603 for(int j=seqStart;j<seqStop;j++){
604 for(int k=0;k<numOTUs;k++){
605 dist[index] = childDist[index];
613 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
615 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
617 for(int i=1;i<ncpus;i++){
618 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
623 for(int i=1;i<ncpus;i++){
624 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
630 if (m->control_pressed) { break; }
632 m->mothurOut("\nFinalizing...\n");
635 if (m->control_pressed) { break; }
639 vector<int> otuCounts(numOTUs, 0);
640 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
641 calcCentroidsDriver(0, numOTUs);
643 if (m->control_pressed) { break; }
645 writeQualities(otuCounts); if (m->control_pressed) { break; }
646 writeSequences(otuCounts); if (m->control_pressed) { break; }
647 writeNames(otuCounts); if (m->control_pressed) { break; }
648 writeClusters(otuCounts); if (m->control_pressed) { break; }
649 writeGroups(); if (m->control_pressed) { break; }
652 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
658 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
659 if(abort){ return 0; }
662 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
664 for(int i=0;i<numFiles;i++){
666 if (m->control_pressed) { break; }
668 //Now into the pyrodist part
672 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
673 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
674 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
675 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
677 flowDataIntI.resize(numSeqs * numFlowCells);
678 flowDataPrI.resize(numSeqs * numFlowCells);
679 mapUniqueToSeq.resize(numSeqs);
680 mapSeqToUnique.resize(numSeqs);
681 lengths.resize(numSeqs);
682 jointLookUp.resize(NUMBINS * NUMBINS);
684 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
685 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
686 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
687 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
688 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
689 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
690 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
692 flowFileName = string(fileName);
693 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
694 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
696 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
698 if (m->control_pressed) { break; }
701 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
703 //Now into the pyronoise part
704 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
706 singleLookUp.resize(HOMOPS * NUMBINS);
707 uniqueFlowgrams.resize(numUniques * numFlowCells);
708 weight.resize(numOTUs);
709 centroids.resize(numOTUs);
710 change.resize(numOTUs);
711 dist.assign(numOTUs * numSeqs, 0);
712 nSeqsPerOTU.resize(numOTUs);
713 cumNumSeqs.resize(numOTUs);
715 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
716 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
717 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
719 int startOTU = pid * numOTUs / ncpus;
720 int endOTU = (pid + 1) * numOTUs / ncpus;
722 int startSeq = pid * numSeqs / ncpus;
723 int endSeq = (pid + 1) * numSeqs /ncpus;
729 if (m->control_pressed) { break; }
731 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
732 singleTau.assign(total, 0.0000);
733 seqNumber.assign(total, 0);
734 seqIndex.assign(total, 0);
736 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
737 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
738 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
739 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
740 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
741 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
742 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
744 calcCentroidsDriver(startOTU, endOTU);
746 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
747 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
749 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
750 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
751 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
753 vector<int> otuIndex(total, 0);
754 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
755 total = otuIndex.size();
757 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
758 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
759 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
760 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
761 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
763 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
768 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
770 MPI_Barrier(MPI_COMM_WORLD);
773 if(compositeFASTAFileName != ""){
774 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
775 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
778 m->mothurOutEndLine();
779 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
780 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
781 m->mothurOutEndLine();
786 catch(exception& e) {
787 m->errorOut(e, "ShhherCommand", "execute");
791 /**************************************************************************************************/
792 string ShhherCommand::createNamesFile(){
795 vector<string> duplicateNames(numUniques, "");
796 for(int i=0;i<numSeqs;i++){
797 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
799 map<string, string> variables;
800 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName));
801 string nameFileName = getOutputFileName("name",variables);
804 m->openOutputFile(nameFileName, nameFile);
806 for(int i=0;i<numUniques;i++){
808 if (m->control_pressed) { break; }
810 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
811 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
817 catch(exception& e) {
818 m->errorOut(e, "ShhherCommand", "createNamesFile");
822 /**************************************************************************************************/
824 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
826 ostringstream outStream;
827 outStream.setf(ios::fixed, ios::floatfield);
828 outStream.setf(ios::dec, ios::basefield);
829 outStream.setf(ios::showpoint);
830 outStream.precision(6);
832 int begTime = time(NULL);
833 double begClock = clock();
835 for(int i=startSeq;i<stopSeq;i++){
837 if (m->control_pressed) { break; }
839 for(int j=0;j<i;j++){
840 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
842 if(flowDistance < 1e-6){
843 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
845 else if(flowDistance <= cutoff){
846 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
850 m->mothurOut(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
854 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
855 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
857 if (m->control_pressed) { return fDistFileName; }
859 m->mothurOut(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
861 ofstream distFile(fDistFileName.c_str());
862 distFile << outStream.str();
865 return fDistFileName;
867 catch(exception& e) {
868 m->errorOut(e, "ShhherCommand", "flowDistMPI");
872 /**************************************************************************************************/
874 void ShhherCommand::getOTUData(string listFileName){
878 m->openInputFile(listFileName, listFile);
881 listFile >> label >> numOTUs;
883 otuData.assign(numSeqs, 0);
884 cumNumSeqs.assign(numOTUs, 0);
885 nSeqsPerOTU.assign(numOTUs, 0);
886 aaP.clear();aaP.resize(numOTUs);
892 string singleOTU = "";
894 for(int i=0;i<numOTUs;i++){
896 if (m->control_pressed) { break; }
898 listFile >> singleOTU;
900 istringstream otuString(singleOTU);
906 for(int j=0;j<singleOTU.length();j++){
907 char letter = otuString.get();
913 map<string,int>::iterator nmIt = nameMap.find(seqName);
914 int index = nmIt->second;
920 aaP[i].push_back(index);
925 map<string,int>::iterator nmIt = nameMap.find(seqName);
927 int index = nmIt->second;
932 aaP[i].push_back(index);
937 sort(aaP[i].begin(), aaP[i].end());
938 for(int j=0;j<nSeqsPerOTU[i];j++){
939 seqNumber.push_back(aaP[i][j]);
941 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
948 for(int i=1;i<numOTUs;i++){
949 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
952 seqIndex = seqNumber;
957 catch(exception& e) {
958 m->errorOut(e, "ShhherCommand", "getOTUData");
963 /**************************************************************************************************/
965 void ShhherCommand::initPyroCluster(){
967 if (numOTUs < processors) { processors = 1; }
969 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
971 dist.assign(numSeqs * numOTUs, 0);
972 change.assign(numOTUs, 1);
973 centroids.assign(numOTUs, -1);
974 weight.assign(numOTUs, 0);
975 singleTau.assign(numSeqs, 1.0);
977 nSeqsBreaks.assign(processors+1, 0);
978 nOTUsBreaks.assign(processors+1, 0);
980 if (m->debug) { m->mothurOut("[DEBUG]: made it through the memory allocation.\n"); }
983 for(int i=0;i<processors;i++){
984 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
985 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
987 nSeqsBreaks[processors] = numSeqs;
988 nOTUsBreaks[processors] = numOTUs;
990 catch(exception& e) {
991 m->errorOut(e, "ShhherCommand", "initPyroCluster");
996 /**************************************************************************************************/
998 void ShhherCommand::fill(){
1001 for(int i=0;i<numOTUs;i++){
1003 if (m->control_pressed) { break; }
1005 cumNumSeqs[i] = index;
1006 for(int j=0;j<nSeqsPerOTU[i];j++){
1007 seqNumber[index] = aaP[i][j];
1008 seqIndex[index] = aaI[i][j];
1014 catch(exception& e) {
1015 m->errorOut(e, "ShhherCommand", "fill");
1020 /**************************************************************************************************/
1022 void ShhherCommand::getFlowData(){
1025 m->openInputFile(flowFileName, flowFile);
1028 seqNameVector.clear();
1030 flowDataIntI.clear();
1034 int currentNumFlowCells;
1039 flowFile >> numFlowTest;
1041 if (!m->isContainingOnlyDigits(numFlowTest)) { m->mothurOut("[ERROR]: expected a number and got " + numFlowTest + ", quitting. Did you use the flow parameter instead of the file parameter?"); m->mothurOutEndLine(); exit(1); }
1042 else { convert(numFlowTest, numFlowCells); }
1044 int index = 0;//pcluster
1045 while(!flowFile.eof()){
1047 if (m->control_pressed) { break; }
1049 flowFile >> seqName >> currentNumFlowCells;
1050 lengths.push_back(currentNumFlowCells);
1052 seqNameVector.push_back(seqName);
1053 nameMap[seqName] = index++;//pcluster
1055 for(int i=0;i<numFlowCells;i++){
1056 flowFile >> intensity;
1057 if(intensity > 9.99) { intensity = 9.99; }
1058 int intI = int(100 * intensity + 0.0001);
1059 flowDataIntI.push_back(intI);
1061 m->gobble(flowFile);
1065 numSeqs = seqNameVector.size();
1067 for(int i=0;i<numSeqs;i++){
1069 if (m->control_pressed) { break; }
1071 int iNumFlowCells = i * numFlowCells;
1072 for(int j=lengths[i];j<numFlowCells;j++){
1073 flowDataIntI[iNumFlowCells + j] = 0;
1078 catch(exception& e) {
1079 m->errorOut(e, "ShhherCommand", "getFlowData");
1083 /**************************************************************************************************/
1084 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1087 vector<double> newTau(numOTUs,0);
1088 vector<double> norms(numSeqs, 0);
1093 for(int i=startSeq;i<stopSeq;i++){
1095 if (m->control_pressed) { break; }
1097 double offset = 1e8;
1098 int indexOffset = i * numOTUs;
1100 for(int j=0;j<numOTUs;j++){
1102 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1103 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1105 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1106 offset = dist[indexOffset + j];
1110 for(int j=0;j<numOTUs;j++){
1111 if(weight[j] > MIN_WEIGHT){
1112 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1113 norms[i] += newTau[j];
1120 for(int j=0;j<numOTUs;j++){
1122 newTau[j] /= norms[i];
1124 if(newTau[j] > MIN_TAU){
1125 otuIndex.push_back(j);
1126 seqIndex.push_back(i);
1127 singleTau.push_back(newTau[j]);
1133 catch(exception& e) {
1134 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1139 /**************************************************************************************************/
1141 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1146 vector<double> newTau(numOTUs,0);
1147 vector<double> norms(numSeqs, 0);
1148 nSeqsPerOTU.assign(numOTUs, 0);
1150 for(int i=startSeq;i<stopSeq;i++){
1152 if (m->control_pressed) { break; }
1154 int indexOffset = i * numOTUs;
1156 double offset = 1e8;
1158 for(int j=0;j<numOTUs;j++){
1160 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1161 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1164 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1165 offset = dist[indexOffset + j];
1169 for(int j=0;j<numOTUs;j++){
1170 if(weight[j] > MIN_WEIGHT){
1171 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1172 norms[i] += newTau[j];
1179 for(int j=0;j<numOTUs;j++){
1180 newTau[j] /= norms[i];
1183 for(int j=0;j<numOTUs;j++){
1184 if(newTau[j] > MIN_TAU){
1186 int oldTotal = total;
1190 singleTau.resize(total, 0);
1191 seqNumber.resize(total, 0);
1192 seqIndex.resize(total, 0);
1194 singleTau[oldTotal] = newTau[j];
1196 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1197 aaI[j][nSeqsPerOTU[j]] = i;
1205 catch(exception& e) {
1206 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1211 /**************************************************************************************************/
1213 void ShhherCommand::setOTUs(){
1216 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1218 for(int i=0;i<numOTUs;i++){
1220 if (m->control_pressed) { break; }
1222 for(int j=0;j<nSeqsPerOTU[i];j++){
1223 int index = cumNumSeqs[i] + j;
1224 double tauValue = singleTau[seqNumber[index]];
1225 int sIndex = seqIndex[index];
1226 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1230 for(int i=0;i<numSeqs;i++){
1231 double maxTau = -1.0000;
1233 for(int j=0;j<numOTUs;j++){
1234 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1235 maxTau = bigTauMatrix[i * numOTUs + j];
1240 otuData[i] = maxOTU;
1243 nSeqsPerOTU.assign(numOTUs, 0);
1245 for(int i=0;i<numSeqs;i++){
1246 int index = otuData[i];
1248 singleTau[i] = 1.0000;
1251 aaP[index][nSeqsPerOTU[index]] = i;
1252 aaI[index][nSeqsPerOTU[index]] = i;
1254 nSeqsPerOTU[index]++;
1258 catch(exception& e) {
1259 m->errorOut(e, "ShhherCommand", "setOTUs");
1264 /**************************************************************************************************/
1266 void ShhherCommand::getUniques(){
1271 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
1272 uniqueCount.assign(numSeqs, 0); // anWeights
1273 uniqueLengths.assign(numSeqs, 0);
1274 mapSeqToUnique.assign(numSeqs, -1);
1275 mapUniqueToSeq.assign(numSeqs, -1);
1277 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
1279 for(int i=0;i<numSeqs;i++){
1281 if (m->control_pressed) { break; }
1285 vector<short> current(numFlowCells);
1286 for(int j=0;j<numFlowCells;j++){
1287 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
1290 for(int j=0;j<numUniques;j++){
1291 int offset = j * numFlowCells;
1295 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
1296 else { shorterLength = uniqueLengths[j]; }
1298 for(int k=0;k<shorterLength;k++){
1299 if(current[k] != uniqueFlowgrams[offset + k]){
1306 mapSeqToUnique[i] = j;
1309 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
1315 if(index == numUniques){
1316 uniqueLengths[numUniques] = lengths[i];
1317 uniqueCount[numUniques] = 1;
1318 mapSeqToUnique[i] = numUniques;//anMap
1319 mapUniqueToSeq[numUniques] = i;//anF
1321 for(int k=0;k<numFlowCells;k++){
1322 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
1323 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
1329 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
1330 uniqueLengths.resize(numUniques);
1332 flowDataPrI.resize(numSeqs * numFlowCells, 0);
1333 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
1335 catch(exception& e) {
1336 m->errorOut(e, "ShhherCommand", "getUniques");
1341 /**************************************************************************************************/
1343 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
1345 int minLength = lengths[mapSeqToUnique[seqA]];
1346 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
1348 int ANumFlowCells = seqA * numFlowCells;
1349 int BNumFlowCells = seqB * numFlowCells;
1353 for(int i=0;i<minLength;i++){
1355 if (m->control_pressed) { break; }
1357 int flowAIntI = flowDataIntI[ANumFlowCells + i];
1358 float flowAPrI = flowDataPrI[ANumFlowCells + i];
1360 int flowBIntI = flowDataIntI[BNumFlowCells + i];
1361 float flowBPrI = flowDataPrI[BNumFlowCells + i];
1362 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
1365 dist /= (float) minLength;
1368 catch(exception& e) {
1369 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
1374 //**********************************************************************************************************************/
1376 string ShhherCommand::cluster(string distFileName, string namesFileName){
1379 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1380 read->setCutoff(cutoff);
1382 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1383 clusterNameMap->readMap();
1384 read->read(clusterNameMap);
1386 ListVector* list = read->getListVector();
1387 SparseDistanceMatrix* matrix = read->getDMatrix();
1390 delete clusterNameMap;
1392 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1394 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1395 string tag = cluster->getTag();
1397 double clusterCutoff = cutoff;
1398 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1400 if (m->control_pressed) { break; }
1402 cluster->update(clusterCutoff);
1405 list->setLabel(toString(cutoff));
1407 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
1409 m->openOutputFile(listFileName, listFile);
1410 list->print(listFile);
1413 delete matrix; delete cluster; delete rabund; delete list;
1415 return listFileName;
1417 catch(exception& e) {
1418 m->errorOut(e, "ShhherCommand", "cluster");
1423 /**************************************************************************************************/
1425 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1427 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1432 for(int i=start;i<finish;i++){
1434 if (m->control_pressed) { break; }
1438 int minFlowGram = 100000000;
1439 double minFlowValue = 1e8;
1440 change[i] = 0; //FALSE
1442 for(int j=0;j<nSeqsPerOTU[i];j++){
1443 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1446 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1447 vector<double> adF(nSeqsPerOTU[i]);
1448 vector<int> anL(nSeqsPerOTU[i]);
1450 for(int j=0;j<nSeqsPerOTU[i];j++){
1451 int index = cumNumSeqs[i] + j;
1452 int nI = seqIndex[index];
1453 int nIU = mapSeqToUnique[nI];
1456 for(k=0;k<position;k++){
1462 anL[position] = nIU;
1463 adF[position] = 0.0000;
1468 for(int j=0;j<nSeqsPerOTU[i];j++){
1469 int index = cumNumSeqs[i] + j;
1470 int nI = seqIndex[index];
1472 double tauValue = singleTau[seqNumber[index]];
1474 for(int k=0;k<position;k++){
1475 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1476 adF[k] += dist * tauValue;
1480 for(int j=0;j<position;j++){
1481 if(adF[j] < minFlowValue){
1483 minFlowValue = adF[j];
1487 if(centroids[i] != anL[minFlowGram]){
1489 centroids[i] = anL[minFlowGram];
1492 else if(centroids[i] != -1){
1498 catch(exception& e) {
1499 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1504 /**************************************************************************************************/
1506 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1509 int flowAValue = cent * numFlowCells;
1510 int flowBValue = flow * numFlowCells;
1514 for(int i=0;i<length;i++){
1515 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1520 return dist / (double)length;
1522 catch(exception& e) {
1523 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1528 /**************************************************************************************************/
1530 double ShhherCommand::getNewWeights(){
1533 double maxChange = 0;
1535 for(int i=0;i<numOTUs;i++){
1537 if (m->control_pressed) { break; }
1539 double difference = weight[i];
1542 for(int j=0;j<nSeqsPerOTU[i];j++){
1543 int index = cumNumSeqs[i] + j;
1544 double tauValue = singleTau[seqNumber[index]];
1545 weight[i] += tauValue;
1548 difference = fabs(weight[i] - difference);
1549 if(difference > maxChange){ maxChange = difference; }
1553 catch(exception& e) {
1554 m->errorOut(e, "ShhherCommand", "getNewWeights");
1559 /**************************************************************************************************/
1561 double ShhherCommand::getLikelihood(){
1565 vector<long double> P(numSeqs, 0);
1568 for(int i=0;i<numOTUs;i++){
1569 if(weight[i] > MIN_WEIGHT){
1575 for(int i=0;i<numOTUs;i++){
1577 if (m->control_pressed) { break; }
1579 for(int j=0;j<nSeqsPerOTU[i];j++){
1580 int index = cumNumSeqs[i] + j;
1581 int nI = seqIndex[index];
1582 double singleDist = dist[seqNumber[index]];
1584 P[nI] += weight[i] * exp(-singleDist * sigma);
1588 for(int i=0;i<numSeqs;i++){
1589 if(P[i] == 0){ P[i] = DBL_EPSILON; }
1594 nLL = nLL -(double)numSeqs * log(sigma);
1598 catch(exception& e) {
1599 m->errorOut(e, "ShhherCommand", "getNewWeights");
1604 /**************************************************************************************************/
1606 void ShhherCommand::checkCentroids(){
1608 vector<int> unique(numOTUs, 1);
1610 for(int i=0;i<numOTUs;i++){
1611 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1616 for(int i=0;i<numOTUs;i++){
1618 if (m->control_pressed) { break; }
1621 for(int j=i+1;j<numOTUs;j++){
1624 if(centroids[j] == centroids[i]){
1628 weight[i] += weight[j];
1636 catch(exception& e) {
1637 m->errorOut(e, "ShhherCommand", "checkCentroids");
1641 /**************************************************************************************************/
1645 void ShhherCommand::writeQualities(vector<int> otuCounts){
1648 string thisOutputDir = outputDir;
1649 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1650 map<string, string> variables;
1651 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName));
1652 string qualityFileName = getOutputFileName("qfile",variables);
1654 ofstream qualityFile;
1655 m->openOutputFile(qualityFileName, qualityFile);
1657 qualityFile.setf(ios::fixed, ios::floatfield);
1658 qualityFile.setf(ios::showpoint);
1659 qualityFile << setprecision(6);
1661 vector<vector<int> > qualities(numOTUs);
1662 vector<double> pr(HOMOPS, 0);
1665 for(int i=0;i<numOTUs;i++){
1667 if (m->control_pressed) { break; }
1672 if(nSeqsPerOTU[i] > 0){
1673 qualities[i].assign(1024, -1);
1675 while(index < numFlowCells){
1676 double maxPrValue = 1e8;
1677 short maxPrIndex = -1;
1678 double count = 0.0000;
1680 pr.assign(HOMOPS, 0);
1682 for(int j=0;j<nSeqsPerOTU[i];j++){
1683 int lIndex = cumNumSeqs[i] + j;
1684 double tauValue = singleTau[seqNumber[lIndex]];
1685 int sequenceIndex = aaI[i][j];
1686 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1690 for(int s=0;s<HOMOPS;s++){
1691 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1695 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1696 maxPrValue = pr[maxPrIndex];
1698 if(count > MIN_COUNT){
1700 double norm = 0.0000;
1702 for(int s=0;s<HOMOPS;s++){
1703 norm += exp(-(pr[s] - maxPrValue));
1706 for(int s=1;s<=maxPrIndex;s++){
1708 double temp = 0.0000;
1710 U += exp(-(pr[s-1]-maxPrValue))/norm;
1718 temp = floor(-10 * temp);
1719 value = (int)floor(temp);
1720 if(value > 100){ value = 100; }
1722 qualities[i][base] = (int)value;
1732 if(otuCounts[i] > 0){
1733 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1735 int j=4; //need to get past the first four bases
1736 while(qualities[i][j] != -1){
1737 qualityFile << qualities[i][j] << ' ';
1740 qualityFile << endl;
1743 qualityFile.close();
1744 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
1747 catch(exception& e) {
1748 m->errorOut(e, "ShhherCommand", "writeQualities");
1753 /**************************************************************************************************/
1755 void ShhherCommand::writeSequences(vector<int> otuCounts){
1757 string thisOutputDir = outputDir;
1758 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1759 map<string, string> variables;
1760 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1761 string fastaFileName = getOutputFileName("fasta",variables);
1763 m->openOutputFile(fastaFileName, fastaFile);
1765 vector<string> names(numOTUs, "");
1767 for(int i=0;i<numOTUs;i++){
1769 if (m->control_pressed) { break; }
1771 int index = centroids[i];
1773 if(otuCounts[i] > 0){
1774 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
1778 for(int j=0;j<numFlowCells;j++){
1780 char base = flowOrder[j % 4];
1781 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
1786 fastaFile << newSeq.substr(4) << endl;
1791 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
1793 if(compositeFASTAFileName != ""){
1794 m->appendFiles(fastaFileName, compositeFASTAFileName);
1797 catch(exception& e) {
1798 m->errorOut(e, "ShhherCommand", "writeSequences");
1803 /**************************************************************************************************/
1805 void ShhherCommand::writeNames(vector<int> otuCounts){
1807 string thisOutputDir = outputDir;
1808 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1809 map<string, string> variables;
1810 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1811 string nameFileName = getOutputFileName("name",variables);
1813 m->openOutputFile(nameFileName, nameFile);
1815 for(int i=0;i<numOTUs;i++){
1817 if (m->control_pressed) { break; }
1819 if(otuCounts[i] > 0){
1820 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
1822 for(int j=1;j<nSeqsPerOTU[i];j++){
1823 nameFile << ',' << seqNameVector[aaI[i][j]];
1830 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
1833 if(compositeNamesFileName != ""){
1834 m->appendFiles(nameFileName, compositeNamesFileName);
1837 catch(exception& e) {
1838 m->errorOut(e, "ShhherCommand", "writeNames");
1843 /**************************************************************************************************/
1845 void ShhherCommand::writeGroups(){
1847 string thisOutputDir = outputDir;
1848 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1849 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
1850 int pos = fileRoot.find_first_of('.');
1851 string fileGroup = fileRoot;
1852 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
1853 map<string, string> variables;
1854 variables["[filename]"] = thisOutputDir + fileRoot;
1855 string groupFileName = getOutputFileName("group",variables);
1857 m->openOutputFile(groupFileName, groupFile);
1859 for(int i=0;i<numSeqs;i++){
1860 if (m->control_pressed) { break; }
1861 groupFile << seqNameVector[i] << '\t' << fileGroup << endl;
1864 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
1867 catch(exception& e) {
1868 m->errorOut(e, "ShhherCommand", "writeGroups");
1873 /**************************************************************************************************/
1875 void ShhherCommand::writeClusters(vector<int> otuCounts){
1877 string thisOutputDir = outputDir;
1878 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1879 map<string, string> variables;
1880 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1881 string otuCountsFileName = getOutputFileName("counts",variables);
1882 ofstream otuCountsFile;
1883 m->openOutputFile(otuCountsFileName, otuCountsFile);
1885 string bases = flowOrder;
1887 for(int i=0;i<numOTUs;i++){
1889 if (m->control_pressed) {
1892 //output the translated version of the centroid sequence for the otu
1893 if(otuCounts[i] > 0){
1894 int index = centroids[i];
1896 otuCountsFile << "ideal\t";
1897 for(int j=8;j<numFlowCells;j++){
1898 char base = bases[j % 4];
1899 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
1900 otuCountsFile << base;
1903 otuCountsFile << endl;
1905 for(int j=0;j<nSeqsPerOTU[i];j++){
1906 int sequence = aaI[i][j];
1907 otuCountsFile << seqNameVector[sequence] << '\t';
1911 for(int k=0;k<lengths[sequence];k++){
1912 char base = bases[k % 4];
1913 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
1915 for(int s=0;s<freq;s++){
1917 //otuCountsFile << base;
1920 otuCountsFile << newSeq.substr(4) << endl;
1922 otuCountsFile << endl;
1925 otuCountsFile.close();
1926 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
1929 catch(exception& e) {
1930 m->errorOut(e, "ShhherCommand", "writeClusters");
1936 //**********************************************************************************************************************
1938 int ShhherCommand::execute(){
1940 if (abort == true) { if (calledHelp) { return 0; } return 2; }
1942 getSingleLookUp(); if (m->control_pressed) { return 0; }
1943 getJointLookUp(); if (m->control_pressed) { return 0; }
1945 int numFiles = flowFileVector.size();
1947 if (numFiles < processors) { processors = numFiles; }
1949 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1950 if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size()); }
1951 else { createProcesses(flowFileVector); } //each processor processes one file
1953 driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size());
1956 if(compositeFASTAFileName != ""){
1957 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
1958 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
1961 m->mothurOutEndLine();
1962 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
1963 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
1964 m->mothurOutEndLine();
1968 catch(exception& e) {
1969 m->errorOut(e, "ShhherCommand", "execute");
1974 /**************************************************************************************************/
1976 int ShhherCommand::createProcesses(vector<string> filenames){
1978 vector<int> processIDS;
1983 if (filenames.size() < processors) { processors = filenames.size(); }
1985 //divide the groups between the processors
1986 vector<linePair> lines;
1987 vector<int> numFilesToComplete;
1988 int numFilesPerProcessor = filenames.size() / processors;
1989 for (int i = 0; i < processors; i++) {
1990 int startIndex = i * numFilesPerProcessor;
1991 int endIndex = (i+1) * numFilesPerProcessor;
1992 if(i == (processors - 1)){ endIndex = filenames.size(); }
1993 lines.push_back(linePair(startIndex, endIndex));
1994 numFilesToComplete.push_back((endIndex-startIndex));
1997 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1999 //loop through and create all the processes you want
2000 while (process != processors) {
2004 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
2006 }else if (pid == 0){
2007 num = driver(filenames, compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp", lines[process].start, lines[process].end);
2009 //pass numSeqs to parent
2011 string tempFile = compositeFASTAFileName + toString(getpid()) + ".num.temp";
2012 m->openOutputFile(tempFile, out);
2018 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
2019 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
2025 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[0].start, lines[0].end);
2027 //force parent to wait until all the processes are done
2028 for (int i=0;i<processIDS.size();i++) {
2029 int temp = processIDS[i];
2035 //////////////////////////////////////////////////////////////////////////////////////////////////////
2037 /////////////////////// NOT WORKING, ACCESS VIOLATION ON READ OF FLOWGRAMS IN THREAD /////////////////
2039 //////////////////////////////////////////////////////////////////////////////////////////////////////
2040 //Windows version shared memory, so be careful when passing variables through the shhhFlowsData struct.
2041 //Above fork() will clone, so memory is separate, but that's not the case with windows,
2042 //////////////////////////////////////////////////////////////////////////////////////////////////////
2044 vector<shhhFlowsData*> pDataArray;
2045 DWORD dwThreadIdArray[processors-1];
2046 HANDLE hThreadArray[processors-1];
2048 //Create processor worker threads.
2049 for( int i=0; i<processors-1; i++ ){
2050 // Allocate memory for thread data.
2051 string extension = "";
2052 if (i != 0) { extension = toString(i) + ".temp"; }
2054 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);
2055 pDataArray.push_back(tempFlow);
2056 processIDS.push_back(i);
2058 hThreadArray[i] = CreateThread(NULL, 0, ShhhFlowsThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
2061 //using the main process as a worker saves time and memory
2063 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[processors-1].start, lines[processors-1].end);
2065 //Wait until all threads have terminated.
2066 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
2068 //Close all thread handles and free memory allocations.
2069 for(int i=0; i < pDataArray.size(); i++){
2070 for(int j=0; j < pDataArray[i]->outputNames.size(); j++){ outputNames.push_back(pDataArray[i]->outputNames[j]); }
2071 CloseHandle(hThreadArray[i]);
2072 delete pDataArray[i];
2077 for (int i=0;i<processIDS.size();i++) {
2079 string tempFile = compositeFASTAFileName + toString(processIDS[i]) + ".num.temp";
2080 m->openInputFile(tempFile, in);
2084 if (tempNum != numFilesToComplete[i+1]) {
2085 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");
2088 in.close(); m->mothurRemove(tempFile);
2090 if (compositeFASTAFileName != "") {
2091 m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName);
2092 m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName);
2093 m->mothurRemove((compositeFASTAFileName + toString(processIDS[i]) + ".temp"));
2094 m->mothurRemove((compositeNamesFileName + toString(processIDS[i]) + ".temp"));
2101 catch(exception& e) {
2102 m->errorOut(e, "ShhherCommand", "createProcesses");
2106 /**************************************************************************************************/
2108 vector<string> ShhherCommand::parseFlowFiles(string filename){
2110 vector<string> files;
2114 m->openInputFile(filename, in);
2116 int thisNumFLows = 0;
2117 in >> thisNumFLows; m->gobble(in);
2120 if (m->control_pressed) { break; }
2123 string outputFileName = filename + toString(count) + ".temp";
2124 m->openOutputFile(outputFileName, out);
2125 out << thisNumFLows << endl;
2126 files.push_back(outputFileName);
2128 int numLinesWrote = 0;
2129 for (int i = 0; i < largeSize; i++) {
2130 if (in.eof()) { break; }
2131 string line = m->getline(in); m->gobble(in);
2132 out << line << endl;
2137 if (numLinesWrote == 0) { m->mothurRemove(outputFileName); files.pop_back(); }
2142 if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } files.clear(); }
2144 m->mothurOut("\nDivided " + filename + " into " + toString(files.size()) + " files.\n\n");
2148 catch(exception& e) {
2149 m->errorOut(e, "ShhherCommand", "parseFlowFiles");
2153 /**************************************************************************************************/
2155 int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){
2158 int numCompleted = 0;
2160 for(int i=start;i<end;i++){
2162 if (m->control_pressed) { break; }
2164 vector<string> theseFlowFileNames; theseFlowFileNames.push_back(filenames[i]);
2165 if (large) { theseFlowFileNames = parseFlowFiles(filenames[i]); }
2167 if (m->control_pressed) { break; }
2169 double begClock = clock();
2170 unsigned long long begTime;
2172 for (int g = 0; g < theseFlowFileNames.size(); g++) {
2174 string flowFileName = theseFlowFileNames[g];
2175 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n");
2176 m->mothurOut("Reading flowgrams...\n");
2178 vector<string> seqNameVector;
2179 vector<int> lengths;
2180 vector<short> flowDataIntI;
2181 vector<double> flowDataPrI;
2182 map<string, int> nameMap;
2183 vector<short> uniqueFlowgrams;
2184 vector<int> uniqueCount;
2185 vector<int> mapSeqToUnique;
2186 vector<int> mapUniqueToSeq;
2187 vector<int> uniqueLengths;
2190 if (m->debug) { m->mothurOut("[DEBUG]: About to read flowgrams.\n"); }
2191 int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells);
2193 if (m->control_pressed) { break; }
2195 m->mothurOut("Identifying unique flowgrams...\n");
2196 int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI);
2198 if (m->control_pressed) { break; }
2200 m->mothurOut("Calculating distances between flowgrams...\n");
2201 string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
2202 begTime = time(NULL);
2205 flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2207 m->mothurOutEndLine();
2208 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
2211 string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
2212 createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq);
2214 if (m->control_pressed) { break; }
2216 m->mothurOut("\nClustering flowgrams...\n");
2217 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
2218 cluster(listFileName, distFileName, namesFileName);
2220 if (m->control_pressed) { break; }
2222 vector<int> otuData;
2223 vector<int> cumNumSeqs;
2224 vector<int> nSeqsPerOTU;
2225 vector<vector<int> > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2226 vector<vector<int> > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2227 vector<int> seqNumber; //tMaster->anP: the sequence id number sorted by OTU
2228 vector<int> seqIndex; //tMaster->anI; the index that corresponds to seqNumber
2231 int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap);
2233 if (m->control_pressed) { break; }
2235 m->mothurRemove(distFileName);
2236 m->mothurRemove(namesFileName);
2237 m->mothurRemove(listFileName);
2239 vector<double> dist; //adDist - distance of sequences to centroids
2240 vector<short> change; //did the centroid sequence change? 0 = no; 1 = yes
2241 vector<int> centroids; //the representative flowgram for each cluster m
2242 vector<double> weight;
2243 vector<double> singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2244 vector<int> nSeqsBreaks;
2245 vector<int> nOTUsBreaks;
2247 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
2249 dist.assign(numSeqs * numOTUs, 0);
2250 change.assign(numOTUs, 1);
2251 centroids.assign(numOTUs, -1);
2252 weight.assign(numOTUs, 0);
2253 singleTau.assign(numSeqs, 1.0);
2255 nSeqsBreaks.assign(2, 0);
2256 nOTUsBreaks.assign(2, 0);
2259 nSeqsBreaks[1] = numSeqs;
2260 nOTUsBreaks[1] = numOTUs;
2262 if (m->debug) { m->mothurOut("[DEBUG]: done allocating memory, about to denoise.\n"); }
2264 if (m->control_pressed) { break; }
2266 double maxDelta = 0;
2270 begTime = time(NULL);
2272 m->mothurOut("\nDenoising flowgrams...\n");
2273 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
2275 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
2277 if (m->control_pressed) { break; }
2279 double cycClock = clock();
2280 unsigned long long cycTime = time(NULL);
2281 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2283 if (m->control_pressed) { break; }
2285 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2287 if (m->control_pressed) { break; }
2289 maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight);
2291 if (m->control_pressed) { break; }
2293 double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight);
2295 if (m->control_pressed) { break; }
2297 checkCentroids(numOTUs, centroids, weight);
2299 if (m->control_pressed) { break; }
2301 calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths);
2303 if (m->control_pressed) { break; }
2307 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
2311 if (m->control_pressed) { break; }
2313 m->mothurOut("\nFinalizing...\n");
2314 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2316 if (m->control_pressed) { break; }
2318 setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI);
2320 if (m->control_pressed) { break; }
2322 vector<int> otuCounts(numOTUs, 0);
2323 for(int j=0;j<numSeqs;j++) { otuCounts[otuData[j]]++; }
2325 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2327 if (m->control_pressed) { break; }
2329 if ((large) && (g == 0)) { flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; }
2330 string thisOutputDir = outputDir;
2331 if (outputDir == "") { thisOutputDir = m->hasPath(flowFileName); }
2332 map<string, string> variables;
2333 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
2334 string qualityFileName = getOutputFileName("qfile",variables);
2335 string fastaFileName = getOutputFileName("fasta",variables);
2336 string nameFileName = getOutputFileName("name",variables);
2337 string otuCountsFileName = getOutputFileName("counts",variables);
2338 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
2339 int pos = fileRoot.find_first_of('.');
2340 string fileGroup = fileRoot;
2341 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
2342 string groupFileName = getOutputFileName("group",variables);
2345 writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; }
2346 writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; }
2347 writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; }
2348 writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; }
2349 writeGroups(groupFileName, fileGroup, numSeqs, seqNameVector); if (m->control_pressed) { break; }
2353 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0]));
2354 m->appendFiles(qualityFileName, getOutputFileName("qfile",variables));
2355 m->mothurRemove(qualityFileName);
2356 m->appendFiles(fastaFileName, getOutputFileName("fasta",variables));
2357 m->mothurRemove(fastaFileName);
2358 m->appendFiles(nameFileName, getOutputFileName("name",variables));
2359 m->mothurRemove(nameFileName);
2360 m->appendFiles(otuCountsFileName, getOutputFileName("counts",variables));
2361 m->mothurRemove(otuCountsFileName);
2362 m->appendFiles(groupFileName, getOutputFileName("group",variables));
2363 m->mothurRemove(groupFileName);
2365 m->mothurRemove(theseFlowFileNames[g]);
2370 m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
2373 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
2375 return numCompleted;
2377 }catch(exception& e) {
2378 m->errorOut(e, "ShhherCommand", "driver");
2383 /**************************************************************************************************/
2384 int ShhherCommand::getFlowData(string filename, vector<string>& thisSeqNameVector, vector<int>& thisLengths, vector<short>& thisFlowDataIntI, map<string, int>& thisNameMap, int& numFlowCells){
2389 m->openInputFile(filename, flowFile);
2392 int currentNumFlowCells;
2394 thisSeqNameVector.clear();
2395 thisLengths.clear();
2396 thisFlowDataIntI.clear();
2397 thisNameMap.clear();
2400 flowFile >> numFlowTest;
2402 if (!m->isContainingOnlyDigits(numFlowTest)) { m->mothurOut("[ERROR]: expected a number and got " + numFlowTest + ", quitting. Did you use the flow parameter instead of the file parameter?"); m->mothurOutEndLine(); exit(1); }
2403 else { convert(numFlowTest, numFlowCells); }
2405 if (m->debug) { m->mothurOut("[DEBUG]: numFlowCells = " + toString(numFlowCells) + ".\n"); }
2406 int index = 0;//pcluster
2407 while(!flowFile.eof()){
2409 if (m->control_pressed) { break; }
2411 flowFile >> seqName >> currentNumFlowCells;
2413 thisLengths.push_back(currentNumFlowCells);
2415 thisSeqNameVector.push_back(seqName);
2416 thisNameMap[seqName] = index++;//pcluster
2418 if (m->debug) { m->mothurOut("[DEBUG]: seqName = " + seqName + " length = " + toString(currentNumFlowCells) + " index = " + toString(index) + "\n"); }
2420 for(int i=0;i<numFlowCells;i++){
2421 flowFile >> intensity;
2422 if(intensity > 9.99) { intensity = 9.99; }
2423 int intI = int(100 * intensity + 0.0001);
2424 thisFlowDataIntI.push_back(intI);
2426 m->gobble(flowFile);
2430 int numSeqs = thisSeqNameVector.size();
2432 for(int i=0;i<numSeqs;i++){
2434 if (m->control_pressed) { break; }
2436 int iNumFlowCells = i * numFlowCells;
2437 for(int j=thisLengths[i];j<numFlowCells;j++){
2438 thisFlowDataIntI[iNumFlowCells + j] = 0;
2445 catch(exception& e) {
2446 m->errorOut(e, "ShhherCommand", "getFlowData");
2450 /**************************************************************************************************/
2452 int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int stopSeq, vector<int>& mapUniqueToSeq, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2455 ostringstream outStream;
2456 outStream.setf(ios::fixed, ios::floatfield);
2457 outStream.setf(ios::dec, ios::basefield);
2458 outStream.setf(ios::showpoint);
2459 outStream.precision(6);
2461 int begTime = time(NULL);
2462 double begClock = clock();
2464 for(int i=0;i<stopSeq;i++){
2466 if (m->control_pressed) { break; }
2468 for(int j=0;j<i;j++){
2469 float flowDistance = calcPairwiseDist(numFlowCells, mapUniqueToSeq[i], mapUniqueToSeq[j], mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2471 if(flowDistance < 1e-6){
2472 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
2474 else if(flowDistance <= cutoff){
2475 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
2479 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
2480 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2481 m->mothurOutEndLine();
2485 ofstream distFile(distFileName.c_str());
2486 distFile << outStream.str();
2489 if (m->control_pressed) {}
2491 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
2492 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2493 m->mothurOutEndLine();
2498 catch(exception& e) {
2499 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
2503 /**************************************************************************************************/
2505 float ShhherCommand::calcPairwiseDist(int numFlowCells, int seqA, int seqB, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2507 int minLength = lengths[mapSeqToUnique[seqA]];
2508 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
2510 int ANumFlowCells = seqA * numFlowCells;
2511 int BNumFlowCells = seqB * numFlowCells;
2515 for(int i=0;i<minLength;i++){
2517 if (m->control_pressed) { break; }
2519 int flowAIntI = flowDataIntI[ANumFlowCells + i];
2520 float flowAPrI = flowDataPrI[ANumFlowCells + i];
2522 int flowBIntI = flowDataIntI[BNumFlowCells + i];
2523 float flowBPrI = flowDataPrI[BNumFlowCells + i];
2524 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
2527 dist /= (float) minLength;
2530 catch(exception& e) {
2531 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
2536 /**************************************************************************************************/
2538 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){
2541 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
2542 uniqueCount.assign(numSeqs, 0); // anWeights
2543 uniqueLengths.assign(numSeqs, 0);
2544 mapSeqToUnique.assign(numSeqs, -1);
2545 mapUniqueToSeq.assign(numSeqs, -1);
2547 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
2549 for(int i=0;i<numSeqs;i++){
2551 if (m->control_pressed) { break; }
2555 vector<short> current(numFlowCells);
2556 for(int j=0;j<numFlowCells;j++){
2557 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
2560 for(int j=0;j<numUniques;j++){
2561 int offset = j * numFlowCells;
2565 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
2566 else { shorterLength = uniqueLengths[j]; }
2568 for(int k=0;k<shorterLength;k++){
2569 if(current[k] != uniqueFlowgrams[offset + k]){
2576 mapSeqToUnique[i] = j;
2579 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
2585 if(index == numUniques){
2586 uniqueLengths[numUniques] = lengths[i];
2587 uniqueCount[numUniques] = 1;
2588 mapSeqToUnique[i] = numUniques;//anMap
2589 mapUniqueToSeq[numUniques] = i;//anF
2591 for(int k=0;k<numFlowCells;k++){
2592 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
2593 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
2599 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
2600 uniqueLengths.resize(numUniques);
2602 flowDataPrI.resize(numSeqs * numFlowCells, 0);
2603 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
2607 catch(exception& e) {
2608 m->errorOut(e, "ShhherCommand", "getUniques");
2612 /**************************************************************************************************/
2613 int ShhherCommand::createNamesFile(int numSeqs, int numUniques, string filename, vector<string>& seqNameVector, vector<int>& mapSeqToUnique, vector<int>& mapUniqueToSeq){
2616 vector<string> duplicateNames(numUniques, "");
2617 for(int i=0;i<numSeqs;i++){
2618 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
2622 m->openOutputFile(filename, nameFile);
2624 for(int i=0;i<numUniques;i++){
2626 if (m->control_pressed) { break; }
2628 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2629 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2636 catch(exception& e) {
2637 m->errorOut(e, "ShhherCommand", "createNamesFile");
2641 //**********************************************************************************************************************
2643 int ShhherCommand::cluster(string filename, string distFileName, string namesFileName){
2646 ReadMatrix* read = new ReadColumnMatrix(distFileName);
2647 read->setCutoff(cutoff);
2649 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
2650 clusterNameMap->readMap();
2651 read->read(clusterNameMap);
2653 ListVector* list = read->getListVector();
2654 SparseDistanceMatrix* matrix = read->getDMatrix();
2657 delete clusterNameMap;
2659 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
2661 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
2662 string tag = cluster->getTag();
2664 double clusterCutoff = cutoff;
2665 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
2667 if (m->control_pressed) { break; }
2669 cluster->update(clusterCutoff);
2672 list->setLabel(toString(cutoff));
2675 m->openOutputFile(filename, listFile);
2676 list->print(listFile);
2679 delete matrix; delete cluster; delete rabund; delete list;
2683 catch(exception& e) {
2684 m->errorOut(e, "ShhherCommand", "cluster");
2688 /**************************************************************************************************/
2690 int ShhherCommand::getOTUData(int numSeqs, string fileName, vector<int>& otuData,
2691 vector<int>& cumNumSeqs,
2692 vector<int>& nSeqsPerOTU,
2693 vector<vector<int> >& aaP, //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2694 vector<vector<int> >& aaI, //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2695 vector<int>& seqNumber, //tMaster->anP: the sequence id number sorted by OTU
2696 vector<int>& seqIndex,
2697 map<string, int>& nameMap){
2701 m->openInputFile(fileName, listFile);
2705 listFile >> label >> numOTUs;
2707 if (m->debug) { m->mothurOut("[DEBUG]: Getting OTU Data...\n"); }
2709 otuData.assign(numSeqs, 0);
2710 cumNumSeqs.assign(numOTUs, 0);
2711 nSeqsPerOTU.assign(numOTUs, 0);
2712 aaP.clear();aaP.resize(numOTUs);
2718 string singleOTU = "";
2720 for(int i=0;i<numOTUs;i++){
2722 if (m->control_pressed) { break; }
2723 if (m->debug) { m->mothurOut("[DEBUG]: processing OTU " + toString(i) + ".\n"); }
2725 listFile >> singleOTU;
2727 istringstream otuString(singleOTU);
2731 string seqName = "";
2733 for(int j=0;j<singleOTU.length();j++){
2734 char letter = otuString.get();
2740 map<string,int>::iterator nmIt = nameMap.find(seqName);
2741 int index = nmIt->second;
2743 nameMap.erase(nmIt);
2747 aaP[i].push_back(index);
2752 map<string,int>::iterator nmIt = nameMap.find(seqName);
2754 int index = nmIt->second;
2755 nameMap.erase(nmIt);
2759 aaP[i].push_back(index);
2764 sort(aaP[i].begin(), aaP[i].end());
2765 for(int j=0;j<nSeqsPerOTU[i];j++){
2766 seqNumber.push_back(aaP[i][j]);
2768 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
2769 aaP[i].push_back(0);
2775 for(int i=1;i<numOTUs;i++){
2776 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
2779 seqIndex = seqNumber;
2786 catch(exception& e) {
2787 m->errorOut(e, "ShhherCommand", "getOTUData");
2791 /**************************************************************************************************/
2793 int ShhherCommand::calcCentroidsDriver(int numOTUs,
2794 vector<int>& cumNumSeqs,
2795 vector<int>& nSeqsPerOTU,
2796 vector<int>& seqIndex,
2797 vector<short>& change, //did the centroid sequence change? 0 = no; 1 = yes
2798 vector<int>& centroids, //the representative flowgram for each cluster m
2799 vector<double>& singleTau, //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2800 vector<int>& mapSeqToUnique,
2801 vector<short>& uniqueFlowgrams,
2802 vector<short>& flowDataIntI,
2803 vector<int>& lengths,
2805 vector<int>& seqNumber){
2807 //this function gets the most likely homopolymer length at a flow position for a group of sequences
2812 for(int i=0;i<numOTUs;i++){
2814 if (m->control_pressed) { break; }
2818 int minFlowGram = 100000000;
2819 double minFlowValue = 1e8;
2820 change[i] = 0; //FALSE
2822 for(int j=0;j<nSeqsPerOTU[i];j++){
2823 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
2826 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
2827 vector<double> adF(nSeqsPerOTU[i]);
2828 vector<int> anL(nSeqsPerOTU[i]);
2830 for(int j=0;j<nSeqsPerOTU[i];j++){
2831 int index = cumNumSeqs[i] + j;
2832 int nI = seqIndex[index];
2833 int nIU = mapSeqToUnique[nI];
2836 for(k=0;k<position;k++){
2842 anL[position] = nIU;
2843 adF[position] = 0.0000;
2848 for(int j=0;j<nSeqsPerOTU[i];j++){
2849 int index = cumNumSeqs[i] + j;
2850 int nI = seqIndex[index];
2852 double tauValue = singleTau[seqNumber[index]];
2854 for(int k=0;k<position;k++){
2855 double dist = getDistToCentroid(anL[k], nI, lengths[nI], uniqueFlowgrams, flowDataIntI, numFlowCells);
2856 adF[k] += dist * tauValue;
2860 for(int j=0;j<position;j++){
2861 if(adF[j] < minFlowValue){
2863 minFlowValue = adF[j];
2867 if(centroids[i] != anL[minFlowGram]){
2869 centroids[i] = anL[minFlowGram];
2872 else if(centroids[i] != -1){
2880 catch(exception& e) {
2881 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
2885 /**************************************************************************************************/
2887 double ShhherCommand::getDistToCentroid(int cent, int flow, int length, vector<short>& uniqueFlowgrams,
2888 vector<short>& flowDataIntI, int numFlowCells){
2891 int flowAValue = cent * numFlowCells;
2892 int flowBValue = flow * numFlowCells;
2896 for(int i=0;i<length;i++){
2897 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
2902 return dist / (double)length;
2904 catch(exception& e) {
2905 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
2909 /**************************************************************************************************/
2911 double ShhherCommand::getNewWeights(int numOTUs, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU, vector<double>& singleTau, vector<int>& seqNumber, vector<double>& weight){
2914 double maxChange = 0;
2916 for(int i=0;i<numOTUs;i++){
2918 if (m->control_pressed) { break; }
2920 double difference = weight[i];
2923 for(int j=0;j<nSeqsPerOTU[i];j++){
2924 int index = cumNumSeqs[i] + j;
2925 double tauValue = singleTau[seqNumber[index]];
2926 weight[i] += tauValue;
2929 difference = fabs(weight[i] - difference);
2930 if(difference > maxChange){ maxChange = difference; }
2934 catch(exception& e) {
2935 m->errorOut(e, "ShhherCommand", "getNewWeights");
2940 /**************************************************************************************************/
2942 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){
2946 vector<long double> P(numSeqs, 0);
2949 for(int i=0;i<numOTUs;i++){
2950 if(weight[i] > MIN_WEIGHT){
2956 for(int i=0;i<numOTUs;i++){
2958 if (m->control_pressed) { break; }
2960 for(int j=0;j<nSeqsPerOTU[i];j++){
2961 int index = cumNumSeqs[i] + j;
2962 int nI = seqIndex[index];
2963 double singleDist = dist[seqNumber[index]];
2965 P[nI] += weight[i] * exp(-singleDist * sigma);
2969 for(int i=0;i<numSeqs;i++){
2970 if(P[i] == 0){ P[i] = DBL_EPSILON; }
2975 nLL = nLL -(double)numSeqs * log(sigma);
2979 catch(exception& e) {
2980 m->errorOut(e, "ShhherCommand", "getNewWeights");
2985 /**************************************************************************************************/
2987 int ShhherCommand::checkCentroids(int numOTUs, vector<int>& centroids, vector<double>& weight){
2989 vector<int> unique(numOTUs, 1);
2991 for(int i=0;i<numOTUs;i++){
2992 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
2997 for(int i=0;i<numOTUs;i++){
2999 if (m->control_pressed) { break; }
3002 for(int j=i+1;j<numOTUs;j++){
3005 if(centroids[j] == centroids[i]){
3009 weight[i] += weight[j];
3019 catch(exception& e) {
3020 m->errorOut(e, "ShhherCommand", "checkCentroids");
3024 /**************************************************************************************************/
3026 void ShhherCommand::calcNewDistances(int numSeqs, int numOTUs, vector<int>& nSeqsPerOTU, vector<double>& dist,
3027 vector<double>& weight, vector<short>& change, vector<int>& centroids,
3028 vector<vector<int> >& aaP, vector<double>& singleTau, vector<vector<int> >& aaI,
3029 vector<int>& seqNumber, vector<int>& seqIndex,
3030 vector<short>& uniqueFlowgrams,
3031 vector<short>& flowDataIntI, int numFlowCells, vector<int>& lengths){
3036 vector<double> newTau(numOTUs,0);
3037 vector<double> norms(numSeqs, 0);
3038 nSeqsPerOTU.assign(numOTUs, 0);
3040 for(int i=0;i<numSeqs;i++){
3042 if (m->control_pressed) { break; }
3044 int indexOffset = i * numOTUs;
3046 double offset = 1e8;
3048 for(int j=0;j<numOTUs;j++){
3050 if(weight[j] > MIN_WEIGHT && change[j] == 1){
3051 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i], uniqueFlowgrams, flowDataIntI, numFlowCells);
3054 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
3055 offset = dist[indexOffset + j];
3059 for(int j=0;j<numOTUs;j++){
3060 if(weight[j] > MIN_WEIGHT){
3061 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
3062 norms[i] += newTau[j];
3069 for(int j=0;j<numOTUs;j++){
3070 newTau[j] /= norms[i];
3073 for(int j=0;j<numOTUs;j++){
3074 if(newTau[j] > MIN_TAU){
3076 int oldTotal = total;
3080 singleTau.resize(total, 0);
3081 seqNumber.resize(total, 0);
3082 seqIndex.resize(total, 0);
3084 singleTau[oldTotal] = newTau[j];
3086 aaP[j][nSeqsPerOTU[j]] = oldTotal;
3087 aaI[j][nSeqsPerOTU[j]] = i;
3095 catch(exception& e) {
3096 m->errorOut(e, "ShhherCommand", "calcNewDistances");
3100 /**************************************************************************************************/
3102 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){
3105 for(int i=0;i<numOTUs;i++){
3107 if (m->control_pressed) { return 0; }
3109 cumNumSeqs[i] = index;
3110 for(int j=0;j<nSeqsPerOTU[i];j++){
3111 seqNumber[index] = aaP[i][j];
3112 seqIndex[index] = aaI[i][j];
3120 catch(exception& e) {
3121 m->errorOut(e, "ShhherCommand", "fill");
3125 /**************************************************************************************************/
3127 void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector<int>& seqNumber, vector<int>& seqIndex, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU,
3128 vector<int>& otuData, vector<double>& singleTau, vector<double>& dist, vector<vector<int> >& aaP, vector<vector<int> >& aaI){
3131 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
3133 for(int i=0;i<numOTUs;i++){
3135 if (m->control_pressed) { break; }
3137 for(int j=0;j<nSeqsPerOTU[i];j++){
3138 int index = cumNumSeqs[i] + j;
3139 double tauValue = singleTau[seqNumber[index]];
3140 int sIndex = seqIndex[index];
3141 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
3145 for(int i=0;i<numSeqs;i++){
3146 double maxTau = -1.0000;
3148 for(int j=0;j<numOTUs;j++){
3149 if(bigTauMatrix[i * numOTUs + j] > maxTau){
3150 maxTau = bigTauMatrix[i * numOTUs + j];
3155 otuData[i] = maxOTU;
3158 nSeqsPerOTU.assign(numOTUs, 0);
3160 for(int i=0;i<numSeqs;i++){
3161 int index = otuData[i];
3163 singleTau[i] = 1.0000;
3166 aaP[index][nSeqsPerOTU[index]] = i;
3167 aaI[index][nSeqsPerOTU[index]] = i;
3169 nSeqsPerOTU[index]++;
3172 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
3174 catch(exception& e) {
3175 m->errorOut(e, "ShhherCommand", "setOTUs");
3179 /**************************************************************************************************/
3181 void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string qualityFileName, vector<int> otuCounts, vector<int>& nSeqsPerOTU, vector<int>& seqNumber,
3182 vector<double>& singleTau, vector<short>& flowDataIntI, vector<short>& uniqueFlowgrams, vector<int>& cumNumSeqs,
3183 vector<int>& mapUniqueToSeq, vector<string>& seqNameVector, vector<int>& centroids, vector<vector<int> >& aaI){
3187 ofstream qualityFile;
3188 m->openOutputFile(qualityFileName, qualityFile);
3190 qualityFile.setf(ios::fixed, ios::floatfield);
3191 qualityFile.setf(ios::showpoint);
3192 qualityFile << setprecision(6);
3194 vector<vector<int> > qualities(numOTUs);
3195 vector<double> pr(HOMOPS, 0);
3198 for(int i=0;i<numOTUs;i++){
3200 if (m->control_pressed) { break; }
3205 if(nSeqsPerOTU[i] > 0){
3206 qualities[i].assign(1024, -1);
3208 while(index < numFlowCells){
3209 double maxPrValue = 1e8;
3210 short maxPrIndex = -1;
3211 double count = 0.0000;
3213 pr.assign(HOMOPS, 0);
3215 for(int j=0;j<nSeqsPerOTU[i];j++){
3216 int lIndex = cumNumSeqs[i] + j;
3217 double tauValue = singleTau[seqNumber[lIndex]];
3218 int sequenceIndex = aaI[i][j];
3219 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
3223 for(int s=0;s<HOMOPS;s++){
3224 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
3228 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
3229 maxPrValue = pr[maxPrIndex];
3231 if(count > MIN_COUNT){
3233 double norm = 0.0000;
3235 for(int s=0;s<HOMOPS;s++){
3236 norm += exp(-(pr[s] - maxPrValue));
3239 for(int s=1;s<=maxPrIndex;s++){
3241 double temp = 0.0000;
3243 U += exp(-(pr[s-1]-maxPrValue))/norm;
3251 temp = floor(-10 * temp);
3252 value = (int)floor(temp);
3253 if(value > 100){ value = 100; }
3255 qualities[i][base] = (int)value;
3265 if(otuCounts[i] > 0){
3266 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
3268 int j=4; //need to get past the first four bases
3269 while(qualities[i][j] != -1){
3270 qualityFile << qualities[i][j] << ' ';
3271 if (j > qualities[i].size()) { break; }
3274 qualityFile << endl;
3277 qualityFile.close();
3278 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
3281 catch(exception& e) {
3282 m->errorOut(e, "ShhherCommand", "writeQualities");
3287 /**************************************************************************************************/
3289 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){
3293 m->openOutputFile(fastaFileName, fastaFile);
3295 vector<string> names(numOTUs, "");
3297 for(int i=0;i<numOTUs;i++){
3299 if (m->control_pressed) { break; }
3301 int index = centroids[i];
3303 if(otuCounts[i] > 0){
3304 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
3308 for(int j=0;j<numFlowCells;j++){
3310 char base = flowOrder[j % 4];
3311 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
3316 if (newSeq.length() >= 4) { fastaFile << newSeq.substr(4) << endl; }
3317 else { fastaFile << "NNNN" << endl; }
3322 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
3324 if(thisCompositeFASTAFileName != ""){
3325 m->appendFiles(fastaFileName, thisCompositeFASTAFileName);
3328 catch(exception& e) {
3329 m->errorOut(e, "ShhherCommand", "writeSequences");
3334 /**************************************************************************************************/
3336 void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string nameFileName, vector<int> otuCounts, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& nSeqsPerOTU){
3340 m->openOutputFile(nameFileName, nameFile);
3342 for(int i=0;i<numOTUs;i++){
3344 if (m->control_pressed) { break; }
3346 if(otuCounts[i] > 0){
3347 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
3349 for(int j=1;j<nSeqsPerOTU[i];j++){
3350 nameFile << ',' << seqNameVector[aaI[i][j]];
3357 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
3360 if(thisCompositeNamesFileName != ""){
3361 m->appendFiles(nameFileName, thisCompositeNamesFileName);
3364 catch(exception& e) {
3365 m->errorOut(e, "ShhherCommand", "writeNames");
3370 /**************************************************************************************************/
3372 void ShhherCommand::writeGroups(string groupFileName, string fileRoot, int numSeqs, vector<string>& seqNameVector){
3375 m->openOutputFile(groupFileName, groupFile);
3377 for(int i=0;i<numSeqs;i++){
3378 if (m->control_pressed) { break; }
3379 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
3382 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
3385 catch(exception& e) {
3386 m->errorOut(e, "ShhherCommand", "writeGroups");
3391 /**************************************************************************************************/
3393 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){
3395 ofstream otuCountsFile;
3396 m->openOutputFile(otuCountsFileName, otuCountsFile);
3398 string bases = flowOrder;
3400 for(int i=0;i<numOTUs;i++){
3402 if (m->control_pressed) {
3405 //output the translated version of the centroid sequence for the otu
3406 if(otuCounts[i] > 0){
3407 int index = centroids[i];
3409 otuCountsFile << "ideal\t";
3410 for(int j=8;j<numFlowCells;j++){
3411 char base = bases[j % 4];
3412 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
3413 otuCountsFile << base;
3416 otuCountsFile << endl;
3418 for(int j=0;j<nSeqsPerOTU[i];j++){
3419 int sequence = aaI[i][j];
3420 otuCountsFile << seqNameVector[sequence] << '\t';
3424 for(int k=0;k<lengths[sequence];k++){
3425 char base = bases[k % 4];
3426 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
3428 for(int s=0;s<freq;s++){
3430 //otuCountsFile << base;
3434 if (newSeq.length() >= 4) { otuCountsFile << newSeq.substr(4) << endl; }
3435 else { otuCountsFile << "NNNN" << endl; }
3437 otuCountsFile << endl;
3440 otuCountsFile.close();
3441 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
3444 catch(exception& e) {
3445 m->errorOut(e, "ShhherCommand", "writeClusters");
3450 /**************************************************************************************************/
3452 void ShhherCommand::getSingleLookUp(){
3454 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
3455 singleLookUp.assign(HOMOPS * NUMBINS, 0);
3458 ifstream lookUpFile;
3459 m->openInputFile(lookupFileName, lookUpFile);
3461 for(int i=0;i<HOMOPS;i++){
3463 if (m->control_pressed) { break; }
3466 lookUpFile >> logFracFreq;
3468 for(int j=0;j<NUMBINS;j++) {
3469 lookUpFile >> singleLookUp[index];
3475 catch(exception& e) {
3476 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
3481 /**************************************************************************************************/
3483 void ShhherCommand::getJointLookUp(){
3486 // the most likely joint probability (-log) that two intenities have the same polymer length
3487 jointLookUp.resize(NUMBINS * NUMBINS, 0);
3489 for(int i=0;i<NUMBINS;i++){
3491 if (m->control_pressed) { break; }
3493 for(int j=0;j<NUMBINS;j++){
3495 double minSum = 100000000;
3497 for(int k=0;k<HOMOPS;k++){
3498 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
3500 if(sum < minSum) { minSum = sum; }
3502 jointLookUp[i * NUMBINS + j] = minSum;
3506 catch(exception& e) {
3507 m->errorOut(e, "ShhherCommand", "getJointLookUp");
3512 /**************************************************************************************************/
3514 double ShhherCommand::getProbIntensity(int intIntensity){
3517 double minNegLogProb = 100000000;
3520 for(int i=0;i<HOMOPS;i++){//loop signal strength
3522 if (m->control_pressed) { break; }
3524 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
3525 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
3528 return minNegLogProb;
3530 catch(exception& e) {
3531 m->errorOut(e, "ShhherCommand", "getProbIntensity");