5 * Created by Pat Schloss on 12/27/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "shhhercommand.h"
12 //**********************************************************************************************************************
13 vector<string> ShhherCommand::setParameters(){
15 CommandParameter pflow("flow", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pflow);
16 CommandParameter pfile("file", "InputTypes", "", "", "none", "fileflow", "none",false,false); parameters.push_back(pfile);
17 CommandParameter plookup("lookup", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(plookup);
18 CommandParameter pcutoff("cutoff", "Number", "", "0.01", "", "", "",false,false); parameters.push_back(pcutoff);
19 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
20 CommandParameter pmaxiter("maxiter", "Number", "", "1000", "", "", "",false,false); parameters.push_back(pmaxiter);
21 CommandParameter plarge("large", "Number", "", "-1", "", "", "",false,false); parameters.push_back(plarge);
22 CommandParameter psigma("sigma", "Number", "", "60", "", "", "",false,false); parameters.push_back(psigma);
23 CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "",false,false); parameters.push_back(pmindelta);
24 CommandParameter porder("order", "String", "", "", "", "", "",false,false); parameters.push_back(porder);
25 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
26 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
28 vector<string> myArray;
29 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
33 m->errorOut(e, "ShhherCommand", "setParameters");
37 //**********************************************************************************************************************
38 string ShhherCommand::getHelpString(){
40 string helpString = "";
41 helpString += "The shhh.flows command reads a file containing flowgrams and creates a file of corrected sequences.\n";
45 m->errorOut(e, "ShhherCommand", "getHelpString");
49 //**********************************************************************************************************************
50 string ShhherCommand::getOutputFileNameTag(string type, string inputName=""){
52 string outputFileName = "";
53 map<string, vector<string> >::iterator it;
55 //is this a type this command creates
56 it = outputTypes.find(type);
57 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
59 if (type == "fasta") { outputFileName = "shhh.fasta"; }
60 else if (type == "name") { outputFileName = "shhh.names"; }
61 else if (type == "group") { outputFileName = "shhh.groups"; }
62 else if (type == "counts") { outputFileName = "shhh.counts"; }
63 else if (type == "qfile") { outputFileName = "shhh.qual"; }
64 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
66 return outputFileName;
69 m->errorOut(e, "ShhherCommand", "getOutputFileNameTag");
73 //**********************************************************************************************************************
75 ShhherCommand::ShhherCommand(){
77 abort = true; calledHelp = true;
80 //initialize outputTypes
81 vector<string> tempOutNames;
82 outputTypes["fasta"] = tempOutNames;
83 outputTypes["name"] = tempOutNames;
84 outputTypes["group"] = tempOutNames;
85 outputTypes["counts"] = tempOutNames;
86 outputTypes["qfile"] = tempOutNames;
90 m->errorOut(e, "ShhherCommand", "ShhherCommand");
95 //**********************************************************************************************************************
97 ShhherCommand::ShhherCommand(string option) {
101 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
102 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
106 abort = false; calledHelp = false;
108 //allow user to run help
109 if(option == "help") { help(); abort = true; calledHelp = true; }
110 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
113 vector<string> myArray = setParameters();
115 OptionParser parser(option);
116 map<string,string> parameters = parser.getParameters();
118 ValidParameters validParameter;
119 map<string,string>::iterator it;
121 //check to make sure all parameters are valid for command
122 for (it = parameters.begin(); it != parameters.end(); it++) {
123 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
126 //initialize outputTypes
127 vector<string> tempOutNames;
128 outputTypes["fasta"] = tempOutNames;
129 outputTypes["name"] = tempOutNames;
130 outputTypes["group"] = tempOutNames;
131 outputTypes["counts"] = tempOutNames;
132 outputTypes["qfile"] = tempOutNames;
135 //if the user changes the input directory command factory will send this info to us in the output parameter
136 string inputDir = validParameter.validFile(parameters, "inputdir", false);
137 if (inputDir == "not found"){ inputDir = ""; }
140 it = parameters.find("flow");
141 //user has given a template file
142 if(it != parameters.end()){
143 path = m->hasPath(it->second);
144 //if the user has not given a path then, add inputdir. else leave path alone.
145 if (path == "") { parameters["flow"] = inputDir + it->second; }
148 it = parameters.find("lookup");
149 //user has given a template file
150 if(it != parameters.end()){
151 path = m->hasPath(it->second);
152 //if the user has not given a path then, add inputdir. else leave path alone.
153 if (path == "") { parameters["lookup"] = inputDir + it->second; }
156 it = parameters.find("file");
157 //user has given a template file
158 if(it != parameters.end()){
159 path = m->hasPath(it->second);
160 //if the user has not given a path then, add inputdir. else leave path alone.
161 if (path == "") { parameters["file"] = inputDir + it->second; }
165 //if the user changes the output directory command factory will send this info to us in the output parameter
166 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
168 //check for required parameters
169 flowFileName = validParameter.validFile(parameters, "flow", true);
170 flowFilesFileName = validParameter.validFile(parameters, "file", true);
171 if (flowFileName == "not found" && flowFilesFileName == "not found") {
172 m->mothurOut("values for either flow or file must be provided for the shhh.flows command.");
173 m->mothurOutEndLine();
176 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
178 if(flowFileName != "not found"){
179 compositeFASTAFileName = "";
180 compositeNamesFileName = "";
185 string thisoutputDir = outputDir;
186 if (outputDir == "") { thisoutputDir = m->hasPath(flowFilesFileName); } //if user entered a file with a path then preserve it
188 //we want to rip off .files, and also .flow if its there
189 string fileroot = m->getRootName(m->getSimpleName(flowFilesFileName));
190 if (fileroot[fileroot.length()-1] == '.') { fileroot = fileroot.substr(0, fileroot.length()-1); } //rip off dot
191 string extension = m->getExtension(fileroot);
192 if (extension == ".flow") { fileroot = m->getRootName(fileroot); }
193 else { fileroot += "."; } //add back if needed
195 compositeFASTAFileName = thisoutputDir + fileroot + "shhh.fasta";
196 m->openOutputFile(compositeFASTAFileName, temp);
199 compositeNamesFileName = thisoutputDir + fileroot + "shhh.names";
200 m->openOutputFile(compositeNamesFileName, temp);
204 if(flowFilesFileName != "not found"){
207 ifstream flowFilesFile;
208 m->openInputFile(flowFilesFileName, flowFilesFile);
209 while(flowFilesFile){
210 fName = m->getline(flowFilesFile);
212 //test if file is valid
214 int ableToOpen = m->openInputFile(fName, in, "noerror");
216 if (ableToOpen == 1) {
217 if (inputDir != "") { //default path is set
218 string tryPath = inputDir + fName;
219 m->mothurOut("Unable to open " + fName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
221 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
227 if (ableToOpen == 1) {
228 if (m->getDefaultPath() != "") { //default path is set
229 string tryPath = m->getDefaultPath() + m->getSimpleName(fName);
230 m->mothurOut("Unable to open " + fName + ". Trying default " + tryPath); m->mothurOutEndLine();
232 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
238 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
239 if (ableToOpen == 1) {
240 string exepath = m->argv;
241 string tempPath = exepath;
242 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
243 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
245 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(fName);
246 m->mothurOut("Unable to open " + fName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
248 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
253 if (ableToOpen == 1) { m->mothurOut("Unable to open " + fName + ". Disregarding. "); m->mothurOutEndLine(); }
254 else { flowFileVector.push_back(fName); }
255 m->gobble(flowFilesFile);
257 flowFilesFile.close();
258 if (flowFileVector.size() == 0) { m->mothurOut("[ERROR]: no valid files."); m->mothurOutEndLine(); abort = true; }
261 if (outputDir == "") { outputDir = m->hasPath(flowFileName); }
262 flowFileVector.push_back(flowFileName);
265 //check for optional parameter and set defaults
266 // ...at some point should added some additional type checking...
268 temp = validParameter.validFile(parameters, "lookup", true);
269 if (temp == "not found") {
270 lookupFileName = "LookUp_Titanium.pat";
274 ableToOpen = m->openInputFile(lookupFileName, in, "noerror");
277 //if you can't open it, try input location
278 if (ableToOpen == 1) {
279 if (inputDir != "") { //default path is set
280 string tryPath = inputDir + lookupFileName;
281 m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
283 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
285 lookupFileName = tryPath;
289 //if you can't open it, try default location
290 if (ableToOpen == 1) {
291 if (m->getDefaultPath() != "") { //default path is set
292 string tryPath = m->getDefaultPath() + m->getSimpleName(lookupFileName);
293 m->mothurOut("Unable to open " + lookupFileName + ". Trying default " + tryPath); m->mothurOutEndLine();
295 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
297 lookupFileName = tryPath;
301 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
302 if (ableToOpen == 1) {
303 string exepath = m->argv;
304 string tempPath = exepath;
305 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
306 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
308 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
309 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
311 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
313 lookupFileName = tryPath;
316 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
318 else if(temp == "not open") {
320 lookupFileName = validParameter.validFile(parameters, "lookup", false);
322 //if you can't open it its not inputDir, try mothur excutable location
323 string exepath = m->argv;
324 string tempPath = exepath;
325 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
326 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
328 string tryPath = m->getFullPathName(exepath) + lookupFileName;
329 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
331 int ableToOpen = m->openInputFile(tryPath, in2, "noerror");
333 lookupFileName = tryPath;
335 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
336 }else { lookupFileName = temp; }
338 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
339 m->setProcessors(temp);
340 m->mothurConvert(temp, processors);
342 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
343 m->mothurConvert(temp, cutoff);
345 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
346 m->mothurConvert(temp, minDelta);
348 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
349 m->mothurConvert(temp, maxIters);
351 temp = validParameter.validFile(parameters, "large", false); if (temp == "not found"){ temp = "0"; }
352 m->mothurConvert(temp, largeSize);
353 if (largeSize != 0) { large = true; }
354 else { large = false; }
355 if (largeSize < 0) { m->mothurOut("The value of the large cannot be negative.\n"); }
358 if (large) { m->mothurOut("The large parameter is not available with the MPI-Enabled version.\n"); large=false; }
362 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
363 m->mothurConvert(temp, sigma);
365 flowOrder = validParameter.validFile(parameters, "order", false);
366 if (flowOrder == "not found"){ flowOrder = "TACG"; }
367 else if(flowOrder.length() != 4){
368 m->mothurOut("The value of the order option must be four bases long\n");
376 catch(exception& e) {
377 m->errorOut(e, "ShhherCommand", "ShhherCommand");
381 //**********************************************************************************************************************
383 int ShhherCommand::execute(){
385 if (abort == true) { if (calledHelp) { return 0; } return 2; }
392 for(int i=1;i<ncpus;i++){
393 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
395 if(abort == 1){ return 0; }
399 m->mothurOut("\nGetting preliminary data...\n");
400 getSingleLookUp(); if (m->control_pressed) { return 0; }
401 getJointLookUp(); if (m->control_pressed) { return 0; }
403 vector<string> flowFileVector;
404 if(flowFilesFileName != "not found"){
407 ifstream flowFilesFile;
408 m->openInputFile(flowFilesFileName, flowFilesFile);
409 while(flowFilesFile){
410 fName = m->getline(flowFilesFile);
411 flowFileVector.push_back(fName);
412 m->gobble(flowFilesFile);
416 flowFileVector.push_back(flowFileName);
419 int numFiles = flowFileVector.size();
421 for(int i=1;i<ncpus;i++){
422 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
425 for(int i=0;i<numFiles;i++){
427 if (m->control_pressed) { break; }
429 double begClock = clock();
430 unsigned long long begTime = time(NULL);
432 flowFileName = flowFileVector[i];
434 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
435 m->mothurOut("Reading flowgrams...\n");
438 if (m->control_pressed) { break; }
440 m->mothurOut("Identifying unique flowgrams...\n");
443 if (m->control_pressed) { break; }
445 m->mothurOut("Calculating distances between flowgrams...\n");
447 strcpy(fileName, flowFileName.c_str());
449 for(int i=1;i<ncpus;i++){
450 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
452 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
453 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
454 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
455 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
456 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
457 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
458 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
459 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
460 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
461 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
464 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
466 if (m->control_pressed) { break; }
469 for(int i=1;i<ncpus;i++){
470 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
472 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
473 m->mothurRemove((distFileName + ".temp." + toString(i)));
476 string namesFileName = createNamesFile();
478 if (m->control_pressed) { break; }
480 m->mothurOut("\nClustering flowgrams...\n");
481 string listFileName = cluster(distFileName, namesFileName);
483 if (m->control_pressed) { break; }
487 getOTUData(listFileName);
489 m->mothurRemove(distFileName);
490 m->mothurRemove(namesFileName);
491 m->mothurRemove(listFileName);
493 if (m->control_pressed) { break; }
497 if (m->control_pressed) { break; }
500 for(int i=1;i<ncpus;i++){
501 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
502 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
503 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
504 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
507 if (m->control_pressed) { break; }
512 int numOTUsOnCPU = numOTUs / ncpus;
513 int numSeqsOnCPU = numSeqs / ncpus;
514 m->mothurOut("\nDenoising flowgrams...\n");
515 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
517 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
519 double cycClock = clock();
520 unsigned long long cycTime = time(NULL);
523 if (m->control_pressed) { break; }
525 int total = singleTau.size();
526 for(int i=1;i<ncpus;i++){
527 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
528 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
529 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
531 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
532 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
533 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
534 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
535 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
538 calcCentroidsDriver(0, numOTUsOnCPU);
540 for(int i=1;i<ncpus;i++){
541 int otuStart = i * numOTUs / ncpus;
542 int otuStop = (i + 1) * numOTUs / ncpus;
544 vector<int> tempCentroids(numOTUs, 0);
545 vector<short> tempChange(numOTUs, 0);
547 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
548 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
550 for(int j=otuStart;j<otuStop;j++){
551 centroids[j] = tempCentroids[j];
552 change[j] = tempChange[j];
556 maxDelta = getNewWeights(); if (m->control_pressed) { break; }
557 double nLL = getLikelihood(); if (m->control_pressed) { break; }
558 checkCentroids(); if (m->control_pressed) { break; }
560 for(int i=1;i<ncpus;i++){
561 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
562 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
563 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
566 calcNewDistancesParent(0, numSeqsOnCPU);
568 total = singleTau.size();
570 for(int i=1;i<ncpus;i++){
572 int seqStart = i * numSeqs / ncpus;
573 int seqStop = (i + 1) * numSeqs / ncpus;
575 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
577 vector<int> childSeqIndex(childTotal, 0);
578 vector<double> childSingleTau(childTotal, 0);
579 vector<double> childDist(numSeqs * numOTUs, 0);
580 vector<int> otuIndex(childTotal, 0);
582 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
583 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
584 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
585 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
587 int oldTotal = total;
589 singleTau.resize(total, 0);
590 seqIndex.resize(total, 0);
591 seqNumber.resize(total, 0);
595 for(int j=oldTotal;j<total;j++){
596 int otuI = otuIndex[childIndex];
597 int seqI = childSeqIndex[childIndex];
599 singleTau[j] = childSingleTau[childIndex];
601 aaP[otuI][nSeqsPerOTU[otuI]] = j;
602 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
607 int index = seqStart * numOTUs;
608 for(int j=seqStart;j<seqStop;j++){
609 for(int k=0;k<numOTUs;k++){
610 dist[index] = childDist[index];
618 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
620 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
622 for(int i=1;i<ncpus;i++){
623 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
628 for(int i=1;i<ncpus;i++){
629 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
635 if (m->control_pressed) { break; }
637 m->mothurOut("\nFinalizing...\n");
640 if (m->control_pressed) { break; }
644 vector<int> otuCounts(numOTUs, 0);
645 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
646 calcCentroidsDriver(0, numOTUs);
648 if (m->control_pressed) { break; }
650 writeQualities(otuCounts); if (m->control_pressed) { break; }
651 writeSequences(otuCounts); if (m->control_pressed) { break; }
652 writeNames(otuCounts); if (m->control_pressed) { break; }
653 writeClusters(otuCounts); if (m->control_pressed) { break; }
654 writeGroups(); if (m->control_pressed) { break; }
657 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
663 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
664 if(abort){ return 0; }
667 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
669 for(int i=0;i<numFiles;i++){
671 if (m->control_pressed) { break; }
673 //Now into the pyrodist part
677 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
678 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
679 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
680 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
682 flowDataIntI.resize(numSeqs * numFlowCells);
683 flowDataPrI.resize(numSeqs * numFlowCells);
684 mapUniqueToSeq.resize(numSeqs);
685 mapSeqToUnique.resize(numSeqs);
686 lengths.resize(numSeqs);
687 jointLookUp.resize(NUMBINS * NUMBINS);
689 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
690 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
691 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
692 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
693 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
694 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
695 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
697 flowFileName = string(fileName);
698 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
699 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
701 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
703 if (m->control_pressed) { break; }
706 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
708 //Now into the pyronoise part
709 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
711 singleLookUp.resize(HOMOPS * NUMBINS);
712 uniqueFlowgrams.resize(numUniques * numFlowCells);
713 weight.resize(numOTUs);
714 centroids.resize(numOTUs);
715 change.resize(numOTUs);
716 dist.assign(numOTUs * numSeqs, 0);
717 nSeqsPerOTU.resize(numOTUs);
718 cumNumSeqs.resize(numOTUs);
720 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
721 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
722 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
724 int startOTU = pid * numOTUs / ncpus;
725 int endOTU = (pid + 1) * numOTUs / ncpus;
727 int startSeq = pid * numSeqs / ncpus;
728 int endSeq = (pid + 1) * numSeqs /ncpus;
734 if (m->control_pressed) { break; }
736 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
737 singleTau.assign(total, 0.0000);
738 seqNumber.assign(total, 0);
739 seqIndex.assign(total, 0);
741 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
742 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
743 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
744 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
745 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
746 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
747 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
749 calcCentroidsDriver(startOTU, endOTU);
751 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
752 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
754 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
755 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
756 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
758 vector<int> otuIndex(total, 0);
759 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
760 total = otuIndex.size();
762 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
763 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
764 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
765 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
766 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
768 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
773 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
775 MPI_Barrier(MPI_COMM_WORLD);
778 if(compositeFASTAFileName != ""){
779 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
780 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
783 m->mothurOutEndLine();
784 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
785 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
786 m->mothurOutEndLine();
791 catch(exception& e) {
792 m->errorOut(e, "ShhherCommand", "execute");
796 /**************************************************************************************************/
797 string ShhherCommand::createNamesFile(){
800 vector<string> duplicateNames(numUniques, "");
801 for(int i=0;i<numSeqs;i++){
802 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
805 string nameFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("name");
808 m->openOutputFile(nameFileName, nameFile);
810 for(int i=0;i<numUniques;i++){
812 if (m->control_pressed) { break; }
814 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
815 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
821 catch(exception& e) {
822 m->errorOut(e, "ShhherCommand", "createNamesFile");
826 /**************************************************************************************************/
828 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
830 ostringstream outStream;
831 outStream.setf(ios::fixed, ios::floatfield);
832 outStream.setf(ios::dec, ios::basefield);
833 outStream.setf(ios::showpoint);
834 outStream.precision(6);
836 int begTime = time(NULL);
837 double begClock = clock();
839 for(int i=startSeq;i<stopSeq;i++){
841 if (m->control_pressed) { break; }
843 for(int j=0;j<i;j++){
844 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
846 if(flowDistance < 1e-6){
847 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
849 else if(flowDistance <= cutoff){
850 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
854 m->mothurOut(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
858 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
859 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
861 if (m->control_pressed) { return fDistFileName; }
863 m->mothurOut(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
865 ofstream distFile(fDistFileName.c_str());
866 distFile << outStream.str();
869 return fDistFileName;
871 catch(exception& e) {
872 m->errorOut(e, "ShhherCommand", "flowDistMPI");
876 /**************************************************************************************************/
878 void ShhherCommand::getOTUData(string listFileName){
882 m->openInputFile(listFileName, listFile);
885 listFile >> label >> numOTUs;
887 otuData.assign(numSeqs, 0);
888 cumNumSeqs.assign(numOTUs, 0);
889 nSeqsPerOTU.assign(numOTUs, 0);
890 aaP.clear();aaP.resize(numOTUs);
896 string singleOTU = "";
898 for(int i=0;i<numOTUs;i++){
900 if (m->control_pressed) { break; }
902 listFile >> singleOTU;
904 istringstream otuString(singleOTU);
910 for(int j=0;j<singleOTU.length();j++){
911 char letter = otuString.get();
917 map<string,int>::iterator nmIt = nameMap.find(seqName);
918 int index = nmIt->second;
924 aaP[i].push_back(index);
929 map<string,int>::iterator nmIt = nameMap.find(seqName);
931 int index = nmIt->second;
936 aaP[i].push_back(index);
941 sort(aaP[i].begin(), aaP[i].end());
942 for(int j=0;j<nSeqsPerOTU[i];j++){
943 seqNumber.push_back(aaP[i][j]);
945 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
952 for(int i=1;i<numOTUs;i++){
953 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
956 seqIndex = seqNumber;
961 catch(exception& e) {
962 m->errorOut(e, "ShhherCommand", "getOTUData");
967 /**************************************************************************************************/
969 void ShhherCommand::initPyroCluster(){
971 if (numOTUs < processors) { processors = 1; }
973 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
975 dist.assign(numSeqs * numOTUs, 0);
976 change.assign(numOTUs, 1);
977 centroids.assign(numOTUs, -1);
978 weight.assign(numOTUs, 0);
979 singleTau.assign(numSeqs, 1.0);
981 nSeqsBreaks.assign(processors+1, 0);
982 nOTUsBreaks.assign(processors+1, 0);
984 if (m->debug) { m->mothurOut("[DEBUG]: made it through the memory allocation.\n"); }
987 for(int i=0;i<processors;i++){
988 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
989 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
991 nSeqsBreaks[processors] = numSeqs;
992 nOTUsBreaks[processors] = numOTUs;
994 catch(exception& e) {
995 m->errorOut(e, "ShhherCommand", "initPyroCluster");
1000 /**************************************************************************************************/
1002 void ShhherCommand::fill(){
1005 for(int i=0;i<numOTUs;i++){
1007 if (m->control_pressed) { break; }
1009 cumNumSeqs[i] = index;
1010 for(int j=0;j<nSeqsPerOTU[i];j++){
1011 seqNumber[index] = aaP[i][j];
1012 seqIndex[index] = aaI[i][j];
1018 catch(exception& e) {
1019 m->errorOut(e, "ShhherCommand", "fill");
1024 /**************************************************************************************************/
1026 void ShhherCommand::getFlowData(){
1029 m->openInputFile(flowFileName, flowFile);
1032 seqNameVector.clear();
1034 flowDataIntI.clear();
1038 int currentNumFlowCells;
1043 flowFile >> numFlowTest;
1045 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); }
1046 else { convert(numFlowTest, numFlowCells); }
1048 int index = 0;//pcluster
1049 while(!flowFile.eof()){
1051 if (m->control_pressed) { break; }
1053 flowFile >> seqName >> currentNumFlowCells;
1054 lengths.push_back(currentNumFlowCells);
1056 seqNameVector.push_back(seqName);
1057 nameMap[seqName] = index++;//pcluster
1059 for(int i=0;i<numFlowCells;i++){
1060 flowFile >> intensity;
1061 if(intensity > 9.99) { intensity = 9.99; }
1062 int intI = int(100 * intensity + 0.0001);
1063 flowDataIntI.push_back(intI);
1065 m->gobble(flowFile);
1069 numSeqs = seqNameVector.size();
1071 for(int i=0;i<numSeqs;i++){
1073 if (m->control_pressed) { break; }
1075 int iNumFlowCells = i * numFlowCells;
1076 for(int j=lengths[i];j<numFlowCells;j++){
1077 flowDataIntI[iNumFlowCells + j] = 0;
1082 catch(exception& e) {
1083 m->errorOut(e, "ShhherCommand", "getFlowData");
1087 /**************************************************************************************************/
1088 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1091 vector<double> newTau(numOTUs,0);
1092 vector<double> norms(numSeqs, 0);
1097 for(int i=startSeq;i<stopSeq;i++){
1099 if (m->control_pressed) { break; }
1101 double offset = 1e8;
1102 int indexOffset = i * numOTUs;
1104 for(int j=0;j<numOTUs;j++){
1106 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1107 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1109 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1110 offset = dist[indexOffset + j];
1114 for(int j=0;j<numOTUs;j++){
1115 if(weight[j] > MIN_WEIGHT){
1116 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1117 norms[i] += newTau[j];
1124 for(int j=0;j<numOTUs;j++){
1126 newTau[j] /= norms[i];
1128 if(newTau[j] > MIN_TAU){
1129 otuIndex.push_back(j);
1130 seqIndex.push_back(i);
1131 singleTau.push_back(newTau[j]);
1137 catch(exception& e) {
1138 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1143 /**************************************************************************************************/
1145 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1150 vector<double> newTau(numOTUs,0);
1151 vector<double> norms(numSeqs, 0);
1152 nSeqsPerOTU.assign(numOTUs, 0);
1154 for(int i=startSeq;i<stopSeq;i++){
1156 if (m->control_pressed) { break; }
1158 int indexOffset = i * numOTUs;
1160 double offset = 1e8;
1162 for(int j=0;j<numOTUs;j++){
1164 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1165 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1168 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1169 offset = dist[indexOffset + j];
1173 for(int j=0;j<numOTUs;j++){
1174 if(weight[j] > MIN_WEIGHT){
1175 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1176 norms[i] += newTau[j];
1183 for(int j=0;j<numOTUs;j++){
1184 newTau[j] /= norms[i];
1187 for(int j=0;j<numOTUs;j++){
1188 if(newTau[j] > MIN_TAU){
1190 int oldTotal = total;
1194 singleTau.resize(total, 0);
1195 seqNumber.resize(total, 0);
1196 seqIndex.resize(total, 0);
1198 singleTau[oldTotal] = newTau[j];
1200 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1201 aaI[j][nSeqsPerOTU[j]] = i;
1209 catch(exception& e) {
1210 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1215 /**************************************************************************************************/
1217 void ShhherCommand::setOTUs(){
1220 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1222 for(int i=0;i<numOTUs;i++){
1224 if (m->control_pressed) { break; }
1226 for(int j=0;j<nSeqsPerOTU[i];j++){
1227 int index = cumNumSeqs[i] + j;
1228 double tauValue = singleTau[seqNumber[index]];
1229 int sIndex = seqIndex[index];
1230 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1234 for(int i=0;i<numSeqs;i++){
1235 double maxTau = -1.0000;
1237 for(int j=0;j<numOTUs;j++){
1238 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1239 maxTau = bigTauMatrix[i * numOTUs + j];
1244 otuData[i] = maxOTU;
1247 nSeqsPerOTU.assign(numOTUs, 0);
1249 for(int i=0;i<numSeqs;i++){
1250 int index = otuData[i];
1252 singleTau[i] = 1.0000;
1255 aaP[index][nSeqsPerOTU[index]] = i;
1256 aaI[index][nSeqsPerOTU[index]] = i;
1258 nSeqsPerOTU[index]++;
1262 catch(exception& e) {
1263 m->errorOut(e, "ShhherCommand", "setOTUs");
1268 /**************************************************************************************************/
1270 void ShhherCommand::getUniques(){
1275 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
1276 uniqueCount.assign(numSeqs, 0); // anWeights
1277 uniqueLengths.assign(numSeqs, 0);
1278 mapSeqToUnique.assign(numSeqs, -1);
1279 mapUniqueToSeq.assign(numSeqs, -1);
1281 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
1283 for(int i=0;i<numSeqs;i++){
1285 if (m->control_pressed) { break; }
1289 vector<short> current(numFlowCells);
1290 for(int j=0;j<numFlowCells;j++){
1291 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
1294 for(int j=0;j<numUniques;j++){
1295 int offset = j * numFlowCells;
1299 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
1300 else { shorterLength = uniqueLengths[j]; }
1302 for(int k=0;k<shorterLength;k++){
1303 if(current[k] != uniqueFlowgrams[offset + k]){
1310 mapSeqToUnique[i] = j;
1313 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
1319 if(index == numUniques){
1320 uniqueLengths[numUniques] = lengths[i];
1321 uniqueCount[numUniques] = 1;
1322 mapSeqToUnique[i] = numUniques;//anMap
1323 mapUniqueToSeq[numUniques] = i;//anF
1325 for(int k=0;k<numFlowCells;k++){
1326 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
1327 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
1333 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
1334 uniqueLengths.resize(numUniques);
1336 flowDataPrI.resize(numSeqs * numFlowCells, 0);
1337 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
1339 catch(exception& e) {
1340 m->errorOut(e, "ShhherCommand", "getUniques");
1345 /**************************************************************************************************/
1347 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
1349 int minLength = lengths[mapSeqToUnique[seqA]];
1350 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
1352 int ANumFlowCells = seqA * numFlowCells;
1353 int BNumFlowCells = seqB * numFlowCells;
1357 for(int i=0;i<minLength;i++){
1359 if (m->control_pressed) { break; }
1361 int flowAIntI = flowDataIntI[ANumFlowCells + i];
1362 float flowAPrI = flowDataPrI[ANumFlowCells + i];
1364 int flowBIntI = flowDataIntI[BNumFlowCells + i];
1365 float flowBPrI = flowDataPrI[BNumFlowCells + i];
1366 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
1369 dist /= (float) minLength;
1372 catch(exception& e) {
1373 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
1378 //**********************************************************************************************************************/
1380 string ShhherCommand::cluster(string distFileName, string namesFileName){
1383 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1384 read->setCutoff(cutoff);
1386 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1387 clusterNameMap->readMap();
1388 read->read(clusterNameMap);
1390 ListVector* list = read->getListVector();
1391 SparseDistanceMatrix* matrix = read->getDMatrix();
1394 delete clusterNameMap;
1396 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1398 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1399 string tag = cluster->getTag();
1401 double clusterCutoff = cutoff;
1402 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1404 if (m->control_pressed) { break; }
1406 cluster->update(clusterCutoff);
1409 list->setLabel(toString(cutoff));
1411 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
1413 m->openOutputFile(listFileName, listFile);
1414 list->print(listFile);
1417 delete matrix; delete cluster; delete rabund; delete list;
1419 return listFileName;
1421 catch(exception& e) {
1422 m->errorOut(e, "ShhherCommand", "cluster");
1427 /**************************************************************************************************/
1429 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1431 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1436 for(int i=start;i<finish;i++){
1438 if (m->control_pressed) { break; }
1442 int minFlowGram = 100000000;
1443 double minFlowValue = 1e8;
1444 change[i] = 0; //FALSE
1446 for(int j=0;j<nSeqsPerOTU[i];j++){
1447 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1450 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1451 vector<double> adF(nSeqsPerOTU[i]);
1452 vector<int> anL(nSeqsPerOTU[i]);
1454 for(int j=0;j<nSeqsPerOTU[i];j++){
1455 int index = cumNumSeqs[i] + j;
1456 int nI = seqIndex[index];
1457 int nIU = mapSeqToUnique[nI];
1460 for(k=0;k<position;k++){
1466 anL[position] = nIU;
1467 adF[position] = 0.0000;
1472 for(int j=0;j<nSeqsPerOTU[i];j++){
1473 int index = cumNumSeqs[i] + j;
1474 int nI = seqIndex[index];
1476 double tauValue = singleTau[seqNumber[index]];
1478 for(int k=0;k<position;k++){
1479 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1480 adF[k] += dist * tauValue;
1484 for(int j=0;j<position;j++){
1485 if(adF[j] < minFlowValue){
1487 minFlowValue = adF[j];
1491 if(centroids[i] != anL[minFlowGram]){
1493 centroids[i] = anL[minFlowGram];
1496 else if(centroids[i] != -1){
1502 catch(exception& e) {
1503 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1508 /**************************************************************************************************/
1510 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1513 int flowAValue = cent * numFlowCells;
1514 int flowBValue = flow * numFlowCells;
1518 for(int i=0;i<length;i++){
1519 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1524 return dist / (double)length;
1526 catch(exception& e) {
1527 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1532 /**************************************************************************************************/
1534 double ShhherCommand::getNewWeights(){
1537 double maxChange = 0;
1539 for(int i=0;i<numOTUs;i++){
1541 if (m->control_pressed) { break; }
1543 double difference = weight[i];
1546 for(int j=0;j<nSeqsPerOTU[i];j++){
1547 int index = cumNumSeqs[i] + j;
1548 double tauValue = singleTau[seqNumber[index]];
1549 weight[i] += tauValue;
1552 difference = fabs(weight[i] - difference);
1553 if(difference > maxChange){ maxChange = difference; }
1557 catch(exception& e) {
1558 m->errorOut(e, "ShhherCommand", "getNewWeights");
1563 /**************************************************************************************************/
1565 double ShhherCommand::getLikelihood(){
1569 vector<long double> P(numSeqs, 0);
1572 for(int i=0;i<numOTUs;i++){
1573 if(weight[i] > MIN_WEIGHT){
1579 for(int i=0;i<numOTUs;i++){
1581 if (m->control_pressed) { break; }
1583 for(int j=0;j<nSeqsPerOTU[i];j++){
1584 int index = cumNumSeqs[i] + j;
1585 int nI = seqIndex[index];
1586 double singleDist = dist[seqNumber[index]];
1588 P[nI] += weight[i] * exp(-singleDist * sigma);
1592 for(int i=0;i<numSeqs;i++){
1593 if(P[i] == 0){ P[i] = DBL_EPSILON; }
1598 nLL = nLL -(double)numSeqs * log(sigma);
1602 catch(exception& e) {
1603 m->errorOut(e, "ShhherCommand", "getNewWeights");
1608 /**************************************************************************************************/
1610 void ShhherCommand::checkCentroids(){
1612 vector<int> unique(numOTUs, 1);
1614 for(int i=0;i<numOTUs;i++){
1615 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1620 for(int i=0;i<numOTUs;i++){
1622 if (m->control_pressed) { break; }
1625 for(int j=i+1;j<numOTUs;j++){
1628 if(centroids[j] == centroids[i]){
1632 weight[i] += weight[j];
1640 catch(exception& e) {
1641 m->errorOut(e, "ShhherCommand", "checkCentroids");
1645 /**************************************************************************************************/
1649 void ShhherCommand::writeQualities(vector<int> otuCounts){
1652 string thisOutputDir = outputDir;
1653 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1654 string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("qfile");
1656 ofstream qualityFile;
1657 m->openOutputFile(qualityFileName, qualityFile);
1659 qualityFile.setf(ios::fixed, ios::floatfield);
1660 qualityFile.setf(ios::showpoint);
1661 qualityFile << setprecision(6);
1663 vector<vector<int> > qualities(numOTUs);
1664 vector<double> pr(HOMOPS, 0);
1667 for(int i=0;i<numOTUs;i++){
1669 if (m->control_pressed) { break; }
1674 if(nSeqsPerOTU[i] > 0){
1675 qualities[i].assign(1024, -1);
1677 while(index < numFlowCells){
1678 double maxPrValue = 1e8;
1679 short maxPrIndex = -1;
1680 double count = 0.0000;
1682 pr.assign(HOMOPS, 0);
1684 for(int j=0;j<nSeqsPerOTU[i];j++){
1685 int lIndex = cumNumSeqs[i] + j;
1686 double tauValue = singleTau[seqNumber[lIndex]];
1687 int sequenceIndex = aaI[i][j];
1688 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1692 for(int s=0;s<HOMOPS;s++){
1693 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1697 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1698 maxPrValue = pr[maxPrIndex];
1700 if(count > MIN_COUNT){
1702 double norm = 0.0000;
1704 for(int s=0;s<HOMOPS;s++){
1705 norm += exp(-(pr[s] - maxPrValue));
1708 for(int s=1;s<=maxPrIndex;s++){
1710 double temp = 0.0000;
1712 U += exp(-(pr[s-1]-maxPrValue))/norm;
1720 temp = floor(-10 * temp);
1721 value = (int)floor(temp);
1722 if(value > 100){ value = 100; }
1724 qualities[i][base] = (int)value;
1734 if(otuCounts[i] > 0){
1735 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1737 int j=4; //need to get past the first four bases
1738 while(qualities[i][j] != -1){
1739 qualityFile << qualities[i][j] << ' ';
1742 qualityFile << endl;
1745 qualityFile.close();
1746 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
1749 catch(exception& e) {
1750 m->errorOut(e, "ShhherCommand", "writeQualities");
1755 /**************************************************************************************************/
1757 void ShhherCommand::writeSequences(vector<int> otuCounts){
1759 string thisOutputDir = outputDir;
1760 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1761 string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("fasta");
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 string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("name");
1811 m->openOutputFile(nameFileName, nameFile);
1813 for(int i=0;i<numOTUs;i++){
1815 if (m->control_pressed) { break; }
1817 if(otuCounts[i] > 0){
1818 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
1820 for(int j=1;j<nSeqsPerOTU[i];j++){
1821 nameFile << ',' << seqNameVector[aaI[i][j]];
1828 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
1831 if(compositeNamesFileName != ""){
1832 m->appendFiles(nameFileName, compositeNamesFileName);
1835 catch(exception& e) {
1836 m->errorOut(e, "ShhherCommand", "writeNames");
1841 /**************************************************************************************************/
1843 void ShhherCommand::writeGroups(){
1845 string thisOutputDir = outputDir;
1846 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1847 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
1848 int pos = fileRoot.find_first_of('.');
1849 string fileGroup = fileRoot;
1850 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
1851 string groupFileName = thisOutputDir + fileRoot + getOutputFileNameTag("group");
1853 m->openOutputFile(groupFileName, groupFile);
1855 for(int i=0;i<numSeqs;i++){
1856 if (m->control_pressed) { break; }
1857 groupFile << seqNameVector[i] << '\t' << fileGroup << endl;
1860 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
1863 catch(exception& e) {
1864 m->errorOut(e, "ShhherCommand", "writeGroups");
1869 /**************************************************************************************************/
1871 void ShhherCommand::writeClusters(vector<int> otuCounts){
1873 string thisOutputDir = outputDir;
1874 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1875 string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) +getOutputFileNameTag("counts");
1876 ofstream otuCountsFile;
1877 m->openOutputFile(otuCountsFileName, otuCountsFile);
1879 string bases = flowOrder;
1881 for(int i=0;i<numOTUs;i++){
1883 if (m->control_pressed) {
1886 //output the translated version of the centroid sequence for the otu
1887 if(otuCounts[i] > 0){
1888 int index = centroids[i];
1890 otuCountsFile << "ideal\t";
1891 for(int j=8;j<numFlowCells;j++){
1892 char base = bases[j % 4];
1893 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
1894 otuCountsFile << base;
1897 otuCountsFile << endl;
1899 for(int j=0;j<nSeqsPerOTU[i];j++){
1900 int sequence = aaI[i][j];
1901 otuCountsFile << seqNameVector[sequence] << '\t';
1905 for(int k=0;k<lengths[sequence];k++){
1906 char base = bases[k % 4];
1907 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
1909 for(int s=0;s<freq;s++){
1911 //otuCountsFile << base;
1914 otuCountsFile << newSeq.substr(4) << endl;
1916 otuCountsFile << endl;
1919 otuCountsFile.close();
1920 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
1923 catch(exception& e) {
1924 m->errorOut(e, "ShhherCommand", "writeClusters");
1930 //**********************************************************************************************************************
1932 int ShhherCommand::execute(){
1934 if (abort == true) { if (calledHelp) { return 0; } return 2; }
1936 getSingleLookUp(); if (m->control_pressed) { return 0; }
1937 getJointLookUp(); if (m->control_pressed) { return 0; }
1939 int numFiles = flowFileVector.size();
1941 if (numFiles < processors) { processors = numFiles; }
1943 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1944 if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size()); }
1945 else { createProcesses(flowFileVector); } //each processor processes one file
1947 driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName, 0, flowFileVector.size());
1950 if(compositeFASTAFileName != ""){
1951 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
1952 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
1955 m->mothurOutEndLine();
1956 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
1957 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
1958 m->mothurOutEndLine();
1962 catch(exception& e) {
1963 m->errorOut(e, "ShhherCommand", "execute");
1968 /**************************************************************************************************/
1970 int ShhherCommand::createProcesses(vector<string> filenames){
1972 vector<int> processIDS;
1977 if (filenames.size() < processors) { processors = filenames.size(); }
1979 //divide the groups between the processors
1980 vector<linePair> lines;
1981 vector<int> numFilesToComplete;
1982 int numFilesPerProcessor = filenames.size() / processors;
1983 for (int i = 0; i < processors; i++) {
1984 int startIndex = i * numFilesPerProcessor;
1985 int endIndex = (i+1) * numFilesPerProcessor;
1986 if(i == (processors - 1)){ endIndex = filenames.size(); }
1987 lines.push_back(linePair(startIndex, endIndex));
1988 numFilesToComplete.push_back((endIndex-startIndex));
1991 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1993 //loop through and create all the processes you want
1994 while (process != processors) {
1998 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
2000 }else if (pid == 0){
2001 num = driver(filenames, compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp", lines[process].start, lines[process].end);
2003 //pass numSeqs to parent
2005 string tempFile = compositeFASTAFileName + toString(getpid()) + ".num.temp";
2006 m->openOutputFile(tempFile, out);
2012 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
2013 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
2019 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[0].start, lines[0].end);
2021 //force parent to wait until all the processes are done
2022 for (int i=0;i<processIDS.size();i++) {
2023 int temp = processIDS[i];
2029 //////////////////////////////////////////////////////////////////////////////////////////////////////
2031 /////////////////////// NOT WORKING, ACCESS VIOLATION ON READ OF FLOWGRAMS IN THREAD /////////////////
2033 //////////////////////////////////////////////////////////////////////////////////////////////////////
2034 //Windows version shared memory, so be careful when passing variables through the shhhFlowsData struct.
2035 //Above fork() will clone, so memory is separate, but that's not the case with windows,
2036 //////////////////////////////////////////////////////////////////////////////////////////////////////
2038 vector<shhhFlowsData*> pDataArray;
2039 DWORD dwThreadIdArray[processors-1];
2040 HANDLE hThreadArray[processors-1];
2042 //Create processor worker threads.
2043 for( int i=0; i<processors-1; i++ ){
2044 // Allocate memory for thread data.
2045 string extension = "";
2046 if (i != 0) { extension = toString(i) + ".temp"; }
2048 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);
2049 pDataArray.push_back(tempFlow);
2050 processIDS.push_back(i);
2052 hThreadArray[i] = CreateThread(NULL, 0, ShhhFlowsThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
2055 //using the main process as a worker saves time and memory
2057 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[processors-1].start, lines[processors-1].end);
2059 //Wait until all threads have terminated.
2060 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
2062 //Close all thread handles and free memory allocations.
2063 for(int i=0; i < pDataArray.size(); i++){
2064 for(int j=0; j < pDataArray[i]->outputNames.size(); j++){ outputNames.push_back(pDataArray[i]->outputNames[j]); }
2065 CloseHandle(hThreadArray[i]);
2066 delete pDataArray[i];
2071 for (int i=0;i<processIDS.size();i++) {
2073 string tempFile = compositeFASTAFileName + toString(processIDS[i]) + ".num.temp";
2074 m->openInputFile(tempFile, in);
2078 if (tempNum != numFilesToComplete[i+1]) {
2079 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");
2082 in.close(); m->mothurRemove(tempFile);
2084 if (compositeFASTAFileName != "") {
2085 m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName);
2086 m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName);
2087 m->mothurRemove((compositeFASTAFileName + toString(processIDS[i]) + ".temp"));
2088 m->mothurRemove((compositeNamesFileName + toString(processIDS[i]) + ".temp"));
2095 catch(exception& e) {
2096 m->errorOut(e, "ShhherCommand", "createProcesses");
2100 /**************************************************************************************************/
2102 vector<string> ShhherCommand::parseFlowFiles(string filename){
2104 vector<string> files;
2108 m->openInputFile(filename, in);
2110 int thisNumFLows = 0;
2111 in >> thisNumFLows; m->gobble(in);
2114 if (m->control_pressed) { break; }
2117 string outputFileName = filename + toString(count) + ".temp";
2118 m->openOutputFile(outputFileName, out);
2119 out << thisNumFLows << endl;
2120 files.push_back(outputFileName);
2122 int numLinesWrote = 0;
2123 for (int i = 0; i < largeSize; i++) {
2124 if (in.eof()) { break; }
2125 string line = m->getline(in); m->gobble(in);
2126 out << line << endl;
2131 if (numLinesWrote == 0) { m->mothurRemove(outputFileName); files.pop_back(); }
2136 if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } files.clear(); }
2138 m->mothurOut("\nDivided " + filename + " into " + toString(files.size()) + " files.\n\n");
2142 catch(exception& e) {
2143 m->errorOut(e, "ShhherCommand", "parseFlowFiles");
2147 /**************************************************************************************************/
2149 int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){
2152 int numCompleted = 0;
2154 for(int i=start;i<end;i++){
2156 if (m->control_pressed) { break; }
2158 vector<string> theseFlowFileNames; theseFlowFileNames.push_back(filenames[i]);
2159 if (large) { theseFlowFileNames = parseFlowFiles(filenames[i]); }
2161 if (m->control_pressed) { break; }
2163 double begClock = clock();
2164 unsigned long long begTime;
2166 for (int g = 0; g < theseFlowFileNames.size(); g++) {
2168 string flowFileName = theseFlowFileNames[g];
2169 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n");
2170 m->mothurOut("Reading flowgrams...\n");
2172 vector<string> seqNameVector;
2173 vector<int> lengths;
2174 vector<short> flowDataIntI;
2175 vector<double> flowDataPrI;
2176 map<string, int> nameMap;
2177 vector<short> uniqueFlowgrams;
2178 vector<int> uniqueCount;
2179 vector<int> mapSeqToUnique;
2180 vector<int> mapUniqueToSeq;
2181 vector<int> uniqueLengths;
2184 if (m->debug) { m->mothurOut("[DEBUG]: About to read flowgrams.\n"); }
2185 int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells);
2187 if (m->control_pressed) { break; }
2189 m->mothurOut("Identifying unique flowgrams...\n");
2190 int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI);
2192 if (m->control_pressed) { break; }
2194 m->mothurOut("Calculating distances between flowgrams...\n");
2195 string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
2196 begTime = time(NULL);
2199 flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2201 m->mothurOutEndLine();
2202 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
2205 string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
2206 createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq);
2208 if (m->control_pressed) { break; }
2210 m->mothurOut("\nClustering flowgrams...\n");
2211 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
2212 cluster(listFileName, distFileName, namesFileName);
2214 if (m->control_pressed) { break; }
2216 vector<int> otuData;
2217 vector<int> cumNumSeqs;
2218 vector<int> nSeqsPerOTU;
2219 vector<vector<int> > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2220 vector<vector<int> > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2221 vector<int> seqNumber; //tMaster->anP: the sequence id number sorted by OTU
2222 vector<int> seqIndex; //tMaster->anI; the index that corresponds to seqNumber
2225 int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap);
2227 if (m->control_pressed) { break; }
2229 m->mothurRemove(distFileName);
2230 m->mothurRemove(namesFileName);
2231 m->mothurRemove(listFileName);
2233 vector<double> dist; //adDist - distance of sequences to centroids
2234 vector<short> change; //did the centroid sequence change? 0 = no; 1 = yes
2235 vector<int> centroids; //the representative flowgram for each cluster m
2236 vector<double> weight;
2237 vector<double> singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2238 vector<int> nSeqsBreaks;
2239 vector<int> nOTUsBreaks;
2241 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
2243 dist.assign(numSeqs * numOTUs, 0);
2244 change.assign(numOTUs, 1);
2245 centroids.assign(numOTUs, -1);
2246 weight.assign(numOTUs, 0);
2247 singleTau.assign(numSeqs, 1.0);
2249 nSeqsBreaks.assign(2, 0);
2250 nOTUsBreaks.assign(2, 0);
2253 nSeqsBreaks[1] = numSeqs;
2254 nOTUsBreaks[1] = numOTUs;
2256 if (m->debug) { m->mothurOut("[DEBUG]: done allocating memory, about to denoise.\n"); }
2258 if (m->control_pressed) { break; }
2260 double maxDelta = 0;
2264 begTime = time(NULL);
2266 m->mothurOut("\nDenoising flowgrams...\n");
2267 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
2269 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
2271 if (m->control_pressed) { break; }
2273 double cycClock = clock();
2274 unsigned long long cycTime = time(NULL);
2275 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2277 if (m->control_pressed) { break; }
2279 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2281 if (m->control_pressed) { break; }
2283 maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight);
2285 if (m->control_pressed) { break; }
2287 double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight);
2289 if (m->control_pressed) { break; }
2291 checkCentroids(numOTUs, centroids, weight);
2293 if (m->control_pressed) { break; }
2295 calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths);
2297 if (m->control_pressed) { break; }
2301 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
2305 if (m->control_pressed) { break; }
2307 m->mothurOut("\nFinalizing...\n");
2308 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2310 if (m->control_pressed) { break; }
2312 setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI);
2314 if (m->control_pressed) { break; }
2316 vector<int> otuCounts(numOTUs, 0);
2317 for(int j=0;j<numSeqs;j++) { otuCounts[otuData[j]]++; }
2319 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2321 if (m->control_pressed) { break; }
2323 if ((large) && (g == 0)) { flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; }
2324 string thisOutputDir = outputDir;
2325 if (outputDir == "") { thisOutputDir = m->hasPath(flowFileName); }
2326 string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("qfile");
2327 string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("fasta");
2328 string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("name");
2329 string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + getOutputFileNameTag("counts");
2330 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
2331 int pos = fileRoot.find_first_of('.');
2332 string fileGroup = fileRoot;
2333 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
2334 string groupFileName = thisOutputDir + fileRoot + getOutputFileNameTag("group");
2337 writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; }
2338 writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; }
2339 writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; }
2340 writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; }
2341 writeGroups(groupFileName, fileGroup, numSeqs, seqNameVector); if (m->control_pressed) { break; }
2345 m->appendFiles(qualityFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("qfile")));
2346 m->mothurRemove(qualityFileName);
2347 m->appendFiles(fastaFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("fasta")));
2348 m->mothurRemove(fastaFileName);
2349 m->appendFiles(nameFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("name")));
2350 m->mothurRemove(nameFileName);
2351 m->appendFiles(otuCountsFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("counts")));
2352 m->mothurRemove(otuCountsFileName);
2353 m->appendFiles(groupFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + getOutputFileNameTag("group")));
2354 m->mothurRemove(groupFileName);
2356 m->mothurRemove(theseFlowFileNames[g]);
2361 m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
2364 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
2366 return numCompleted;
2368 }catch(exception& e) {
2369 m->errorOut(e, "ShhherCommand", "driver");
2374 /**************************************************************************************************/
2375 int ShhherCommand::getFlowData(string filename, vector<string>& thisSeqNameVector, vector<int>& thisLengths, vector<short>& thisFlowDataIntI, map<string, int>& thisNameMap, int& numFlowCells){
2380 m->openInputFile(filename, flowFile);
2383 int currentNumFlowCells;
2385 thisSeqNameVector.clear();
2386 thisLengths.clear();
2387 thisFlowDataIntI.clear();
2388 thisNameMap.clear();
2391 flowFile >> numFlowTest;
2393 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); }
2394 else { convert(numFlowTest, numFlowCells); }
2396 if (m->debug) { m->mothurOut("[DEBUG]: numFlowCells = " + toString(numFlowCells) + ".\n"); }
2397 int index = 0;//pcluster
2398 while(!flowFile.eof()){
2400 if (m->control_pressed) { break; }
2402 flowFile >> seqName >> currentNumFlowCells;
2404 thisLengths.push_back(currentNumFlowCells);
2406 thisSeqNameVector.push_back(seqName);
2407 thisNameMap[seqName] = index++;//pcluster
2409 if (m->debug) { m->mothurOut("[DEBUG]: seqName = " + seqName + " length = " + toString(currentNumFlowCells) + " index = " + toString(index) + "\n"); }
2411 for(int i=0;i<numFlowCells;i++){
2412 flowFile >> intensity;
2413 if(intensity > 9.99) { intensity = 9.99; }
2414 int intI = int(100 * intensity + 0.0001);
2415 thisFlowDataIntI.push_back(intI);
2417 m->gobble(flowFile);
2421 int numSeqs = thisSeqNameVector.size();
2423 for(int i=0;i<numSeqs;i++){
2425 if (m->control_pressed) { break; }
2427 int iNumFlowCells = i * numFlowCells;
2428 for(int j=thisLengths[i];j<numFlowCells;j++){
2429 thisFlowDataIntI[iNumFlowCells + j] = 0;
2436 catch(exception& e) {
2437 m->errorOut(e, "ShhherCommand", "getFlowData");
2441 /**************************************************************************************************/
2443 int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int stopSeq, vector<int>& mapUniqueToSeq, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2446 ostringstream outStream;
2447 outStream.setf(ios::fixed, ios::floatfield);
2448 outStream.setf(ios::dec, ios::basefield);
2449 outStream.setf(ios::showpoint);
2450 outStream.precision(6);
2452 int begTime = time(NULL);
2453 double begClock = clock();
2455 for(int i=0;i<stopSeq;i++){
2457 if (m->control_pressed) { break; }
2459 for(int j=0;j<i;j++){
2460 float flowDistance = calcPairwiseDist(numFlowCells, mapUniqueToSeq[i], mapUniqueToSeq[j], mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2462 if(flowDistance < 1e-6){
2463 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
2465 else if(flowDistance <= cutoff){
2466 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
2470 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
2471 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2472 m->mothurOutEndLine();
2476 ofstream distFile(distFileName.c_str());
2477 distFile << outStream.str();
2480 if (m->control_pressed) {}
2482 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
2483 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
2484 m->mothurOutEndLine();
2489 catch(exception& e) {
2490 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
2494 /**************************************************************************************************/
2496 float ShhherCommand::calcPairwiseDist(int numFlowCells, int seqA, int seqB, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2498 int minLength = lengths[mapSeqToUnique[seqA]];
2499 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
2501 int ANumFlowCells = seqA * numFlowCells;
2502 int BNumFlowCells = seqB * numFlowCells;
2506 for(int i=0;i<minLength;i++){
2508 if (m->control_pressed) { break; }
2510 int flowAIntI = flowDataIntI[ANumFlowCells + i];
2511 float flowAPrI = flowDataPrI[ANumFlowCells + i];
2513 int flowBIntI = flowDataIntI[BNumFlowCells + i];
2514 float flowBPrI = flowDataPrI[BNumFlowCells + i];
2515 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
2518 dist /= (float) minLength;
2521 catch(exception& e) {
2522 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
2527 /**************************************************************************************************/
2529 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){
2532 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
2533 uniqueCount.assign(numSeqs, 0); // anWeights
2534 uniqueLengths.assign(numSeqs, 0);
2535 mapSeqToUnique.assign(numSeqs, -1);
2536 mapUniqueToSeq.assign(numSeqs, -1);
2538 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
2540 for(int i=0;i<numSeqs;i++){
2542 if (m->control_pressed) { break; }
2546 vector<short> current(numFlowCells);
2547 for(int j=0;j<numFlowCells;j++){
2548 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
2551 for(int j=0;j<numUniques;j++){
2552 int offset = j * numFlowCells;
2556 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
2557 else { shorterLength = uniqueLengths[j]; }
2559 for(int k=0;k<shorterLength;k++){
2560 if(current[k] != uniqueFlowgrams[offset + k]){
2567 mapSeqToUnique[i] = j;
2570 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
2576 if(index == numUniques){
2577 uniqueLengths[numUniques] = lengths[i];
2578 uniqueCount[numUniques] = 1;
2579 mapSeqToUnique[i] = numUniques;//anMap
2580 mapUniqueToSeq[numUniques] = i;//anF
2582 for(int k=0;k<numFlowCells;k++){
2583 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
2584 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
2590 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
2591 uniqueLengths.resize(numUniques);
2593 flowDataPrI.resize(numSeqs * numFlowCells, 0);
2594 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
2598 catch(exception& e) {
2599 m->errorOut(e, "ShhherCommand", "getUniques");
2603 /**************************************************************************************************/
2604 int ShhherCommand::createNamesFile(int numSeqs, int numUniques, string filename, vector<string>& seqNameVector, vector<int>& mapSeqToUnique, vector<int>& mapUniqueToSeq){
2607 vector<string> duplicateNames(numUniques, "");
2608 for(int i=0;i<numSeqs;i++){
2609 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
2613 m->openOutputFile(filename, nameFile);
2615 for(int i=0;i<numUniques;i++){
2617 if (m->control_pressed) { break; }
2619 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2620 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2627 catch(exception& e) {
2628 m->errorOut(e, "ShhherCommand", "createNamesFile");
2632 //**********************************************************************************************************************
2634 int ShhherCommand::cluster(string filename, string distFileName, string namesFileName){
2637 ReadMatrix* read = new ReadColumnMatrix(distFileName);
2638 read->setCutoff(cutoff);
2640 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
2641 clusterNameMap->readMap();
2642 read->read(clusterNameMap);
2644 ListVector* list = read->getListVector();
2645 SparseDistanceMatrix* matrix = read->getDMatrix();
2648 delete clusterNameMap;
2650 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
2652 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
2653 string tag = cluster->getTag();
2655 double clusterCutoff = cutoff;
2656 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
2658 if (m->control_pressed) { break; }
2660 cluster->update(clusterCutoff);
2663 list->setLabel(toString(cutoff));
2666 m->openOutputFile(filename, listFile);
2667 list->print(listFile);
2670 delete matrix; delete cluster; delete rabund; delete list;
2674 catch(exception& e) {
2675 m->errorOut(e, "ShhherCommand", "cluster");
2679 /**************************************************************************************************/
2681 int ShhherCommand::getOTUData(int numSeqs, string fileName, vector<int>& otuData,
2682 vector<int>& cumNumSeqs,
2683 vector<int>& nSeqsPerOTU,
2684 vector<vector<int> >& aaP, //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2685 vector<vector<int> >& aaI, //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2686 vector<int>& seqNumber, //tMaster->anP: the sequence id number sorted by OTU
2687 vector<int>& seqIndex,
2688 map<string, int>& nameMap){
2692 m->openInputFile(fileName, listFile);
2696 listFile >> label >> numOTUs;
2698 if (m->debug) { m->mothurOut("[DEBUG]: Getting OTU Data...\n"); }
2700 otuData.assign(numSeqs, 0);
2701 cumNumSeqs.assign(numOTUs, 0);
2702 nSeqsPerOTU.assign(numOTUs, 0);
2703 aaP.clear();aaP.resize(numOTUs);
2709 string singleOTU = "";
2711 for(int i=0;i<numOTUs;i++){
2713 if (m->control_pressed) { break; }
2714 if (m->debug) { m->mothurOut("[DEBUG]: processing OTU " + toString(i) + ".\n"); }
2716 listFile >> singleOTU;
2718 istringstream otuString(singleOTU);
2722 string seqName = "";
2724 for(int j=0;j<singleOTU.length();j++){
2725 char letter = otuString.get();
2731 map<string,int>::iterator nmIt = nameMap.find(seqName);
2732 int index = nmIt->second;
2734 nameMap.erase(nmIt);
2738 aaP[i].push_back(index);
2743 map<string,int>::iterator nmIt = nameMap.find(seqName);
2745 int index = nmIt->second;
2746 nameMap.erase(nmIt);
2750 aaP[i].push_back(index);
2755 sort(aaP[i].begin(), aaP[i].end());
2756 for(int j=0;j<nSeqsPerOTU[i];j++){
2757 seqNumber.push_back(aaP[i][j]);
2759 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
2760 aaP[i].push_back(0);
2766 for(int i=1;i<numOTUs;i++){
2767 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
2770 seqIndex = seqNumber;
2777 catch(exception& e) {
2778 m->errorOut(e, "ShhherCommand", "getOTUData");
2782 /**************************************************************************************************/
2784 int ShhherCommand::calcCentroidsDriver(int numOTUs,
2785 vector<int>& cumNumSeqs,
2786 vector<int>& nSeqsPerOTU,
2787 vector<int>& seqIndex,
2788 vector<short>& change, //did the centroid sequence change? 0 = no; 1 = yes
2789 vector<int>& centroids, //the representative flowgram for each cluster m
2790 vector<double>& singleTau, //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2791 vector<int>& mapSeqToUnique,
2792 vector<short>& uniqueFlowgrams,
2793 vector<short>& flowDataIntI,
2794 vector<int>& lengths,
2796 vector<int>& seqNumber){
2798 //this function gets the most likely homopolymer length at a flow position for a group of sequences
2803 for(int i=0;i<numOTUs;i++){
2805 if (m->control_pressed) { break; }
2809 int minFlowGram = 100000000;
2810 double minFlowValue = 1e8;
2811 change[i] = 0; //FALSE
2813 for(int j=0;j<nSeqsPerOTU[i];j++){
2814 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
2817 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
2818 vector<double> adF(nSeqsPerOTU[i]);
2819 vector<int> anL(nSeqsPerOTU[i]);
2821 for(int j=0;j<nSeqsPerOTU[i];j++){
2822 int index = cumNumSeqs[i] + j;
2823 int nI = seqIndex[index];
2824 int nIU = mapSeqToUnique[nI];
2827 for(k=0;k<position;k++){
2833 anL[position] = nIU;
2834 adF[position] = 0.0000;
2839 for(int j=0;j<nSeqsPerOTU[i];j++){
2840 int index = cumNumSeqs[i] + j;
2841 int nI = seqIndex[index];
2843 double tauValue = singleTau[seqNumber[index]];
2845 for(int k=0;k<position;k++){
2846 double dist = getDistToCentroid(anL[k], nI, lengths[nI], uniqueFlowgrams, flowDataIntI, numFlowCells);
2847 adF[k] += dist * tauValue;
2851 for(int j=0;j<position;j++){
2852 if(adF[j] < minFlowValue){
2854 minFlowValue = adF[j];
2858 if(centroids[i] != anL[minFlowGram]){
2860 centroids[i] = anL[minFlowGram];
2863 else if(centroids[i] != -1){
2871 catch(exception& e) {
2872 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
2876 /**************************************************************************************************/
2878 double ShhherCommand::getDistToCentroid(int cent, int flow, int length, vector<short>& uniqueFlowgrams,
2879 vector<short>& flowDataIntI, int numFlowCells){
2882 int flowAValue = cent * numFlowCells;
2883 int flowBValue = flow * numFlowCells;
2887 for(int i=0;i<length;i++){
2888 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
2893 return dist / (double)length;
2895 catch(exception& e) {
2896 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
2900 /**************************************************************************************************/
2902 double ShhherCommand::getNewWeights(int numOTUs, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU, vector<double>& singleTau, vector<int>& seqNumber, vector<double>& weight){
2905 double maxChange = 0;
2907 for(int i=0;i<numOTUs;i++){
2909 if (m->control_pressed) { break; }
2911 double difference = weight[i];
2914 for(int j=0;j<nSeqsPerOTU[i];j++){
2915 int index = cumNumSeqs[i] + j;
2916 double tauValue = singleTau[seqNumber[index]];
2917 weight[i] += tauValue;
2920 difference = fabs(weight[i] - difference);
2921 if(difference > maxChange){ maxChange = difference; }
2925 catch(exception& e) {
2926 m->errorOut(e, "ShhherCommand", "getNewWeights");
2931 /**************************************************************************************************/
2933 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){
2937 vector<long double> P(numSeqs, 0);
2940 for(int i=0;i<numOTUs;i++){
2941 if(weight[i] > MIN_WEIGHT){
2947 for(int i=0;i<numOTUs;i++){
2949 if (m->control_pressed) { break; }
2951 for(int j=0;j<nSeqsPerOTU[i];j++){
2952 int index = cumNumSeqs[i] + j;
2953 int nI = seqIndex[index];
2954 double singleDist = dist[seqNumber[index]];
2956 P[nI] += weight[i] * exp(-singleDist * sigma);
2960 for(int i=0;i<numSeqs;i++){
2961 if(P[i] == 0){ P[i] = DBL_EPSILON; }
2966 nLL = nLL -(double)numSeqs * log(sigma);
2970 catch(exception& e) {
2971 m->errorOut(e, "ShhherCommand", "getNewWeights");
2976 /**************************************************************************************************/
2978 int ShhherCommand::checkCentroids(int numOTUs, vector<int>& centroids, vector<double>& weight){
2980 vector<int> unique(numOTUs, 1);
2982 for(int i=0;i<numOTUs;i++){
2983 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
2988 for(int i=0;i<numOTUs;i++){
2990 if (m->control_pressed) { break; }
2993 for(int j=i+1;j<numOTUs;j++){
2996 if(centroids[j] == centroids[i]){
3000 weight[i] += weight[j];
3010 catch(exception& e) {
3011 m->errorOut(e, "ShhherCommand", "checkCentroids");
3015 /**************************************************************************************************/
3017 void ShhherCommand::calcNewDistances(int numSeqs, int numOTUs, vector<int>& nSeqsPerOTU, vector<double>& dist,
3018 vector<double>& weight, vector<short>& change, vector<int>& centroids,
3019 vector<vector<int> >& aaP, vector<double>& singleTau, vector<vector<int> >& aaI,
3020 vector<int>& seqNumber, vector<int>& seqIndex,
3021 vector<short>& uniqueFlowgrams,
3022 vector<short>& flowDataIntI, int numFlowCells, vector<int>& lengths){
3027 vector<double> newTau(numOTUs,0);
3028 vector<double> norms(numSeqs, 0);
3029 nSeqsPerOTU.assign(numOTUs, 0);
3031 for(int i=0;i<numSeqs;i++){
3033 if (m->control_pressed) { break; }
3035 int indexOffset = i * numOTUs;
3037 double offset = 1e8;
3039 for(int j=0;j<numOTUs;j++){
3041 if(weight[j] > MIN_WEIGHT && change[j] == 1){
3042 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i], uniqueFlowgrams, flowDataIntI, numFlowCells);
3045 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
3046 offset = dist[indexOffset + j];
3050 for(int j=0;j<numOTUs;j++){
3051 if(weight[j] > MIN_WEIGHT){
3052 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
3053 norms[i] += newTau[j];
3060 for(int j=0;j<numOTUs;j++){
3061 newTau[j] /= norms[i];
3064 for(int j=0;j<numOTUs;j++){
3065 if(newTau[j] > MIN_TAU){
3067 int oldTotal = total;
3071 singleTau.resize(total, 0);
3072 seqNumber.resize(total, 0);
3073 seqIndex.resize(total, 0);
3075 singleTau[oldTotal] = newTau[j];
3077 aaP[j][nSeqsPerOTU[j]] = oldTotal;
3078 aaI[j][nSeqsPerOTU[j]] = i;
3086 catch(exception& e) {
3087 m->errorOut(e, "ShhherCommand", "calcNewDistances");
3091 /**************************************************************************************************/
3093 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){
3096 for(int i=0;i<numOTUs;i++){
3098 if (m->control_pressed) { return 0; }
3100 cumNumSeqs[i] = index;
3101 for(int j=0;j<nSeqsPerOTU[i];j++){
3102 seqNumber[index] = aaP[i][j];
3103 seqIndex[index] = aaI[i][j];
3111 catch(exception& e) {
3112 m->errorOut(e, "ShhherCommand", "fill");
3116 /**************************************************************************************************/
3118 void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector<int>& seqNumber, vector<int>& seqIndex, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU,
3119 vector<int>& otuData, vector<double>& singleTau, vector<double>& dist, vector<vector<int> >& aaP, vector<vector<int> >& aaI){
3122 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
3124 for(int i=0;i<numOTUs;i++){
3126 if (m->control_pressed) { break; }
3128 for(int j=0;j<nSeqsPerOTU[i];j++){
3129 int index = cumNumSeqs[i] + j;
3130 double tauValue = singleTau[seqNumber[index]];
3131 int sIndex = seqIndex[index];
3132 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
3136 for(int i=0;i<numSeqs;i++){
3137 double maxTau = -1.0000;
3139 for(int j=0;j<numOTUs;j++){
3140 if(bigTauMatrix[i * numOTUs + j] > maxTau){
3141 maxTau = bigTauMatrix[i * numOTUs + j];
3146 otuData[i] = maxOTU;
3149 nSeqsPerOTU.assign(numOTUs, 0);
3151 for(int i=0;i<numSeqs;i++){
3152 int index = otuData[i];
3154 singleTau[i] = 1.0000;
3157 aaP[index][nSeqsPerOTU[index]] = i;
3158 aaI[index][nSeqsPerOTU[index]] = i;
3160 nSeqsPerOTU[index]++;
3163 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
3165 catch(exception& e) {
3166 m->errorOut(e, "ShhherCommand", "setOTUs");
3170 /**************************************************************************************************/
3172 void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string qualityFileName, vector<int> otuCounts, vector<int>& nSeqsPerOTU, vector<int>& seqNumber,
3173 vector<double>& singleTau, vector<short>& flowDataIntI, vector<short>& uniqueFlowgrams, vector<int>& cumNumSeqs,
3174 vector<int>& mapUniqueToSeq, vector<string>& seqNameVector, vector<int>& centroids, vector<vector<int> >& aaI){
3178 ofstream qualityFile;
3179 m->openOutputFile(qualityFileName, qualityFile);
3181 qualityFile.setf(ios::fixed, ios::floatfield);
3182 qualityFile.setf(ios::showpoint);
3183 qualityFile << setprecision(6);
3185 vector<vector<int> > qualities(numOTUs);
3186 vector<double> pr(HOMOPS, 0);
3189 for(int i=0;i<numOTUs;i++){
3191 if (m->control_pressed) { break; }
3196 if(nSeqsPerOTU[i] > 0){
3197 qualities[i].assign(1024, -1);
3199 while(index < numFlowCells){
3200 double maxPrValue = 1e8;
3201 short maxPrIndex = -1;
3202 double count = 0.0000;
3204 pr.assign(HOMOPS, 0);
3206 for(int j=0;j<nSeqsPerOTU[i];j++){
3207 int lIndex = cumNumSeqs[i] + j;
3208 double tauValue = singleTau[seqNumber[lIndex]];
3209 int sequenceIndex = aaI[i][j];
3210 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
3214 for(int s=0;s<HOMOPS;s++){
3215 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
3219 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
3220 maxPrValue = pr[maxPrIndex];
3222 if(count > MIN_COUNT){
3224 double norm = 0.0000;
3226 for(int s=0;s<HOMOPS;s++){
3227 norm += exp(-(pr[s] - maxPrValue));
3230 for(int s=1;s<=maxPrIndex;s++){
3232 double temp = 0.0000;
3234 U += exp(-(pr[s-1]-maxPrValue))/norm;
3242 temp = floor(-10 * temp);
3243 value = (int)floor(temp);
3244 if(value > 100){ value = 100; }
3246 qualities[i][base] = (int)value;
3256 if(otuCounts[i] > 0){
3257 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
3259 int j=4; //need to get past the first four bases
3260 while(qualities[i][j] != -1){
3261 qualityFile << qualities[i][j] << ' ';
3262 if (j > qualities[i].size()) { break; }
3265 qualityFile << endl;
3268 qualityFile.close();
3269 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
3272 catch(exception& e) {
3273 m->errorOut(e, "ShhherCommand", "writeQualities");
3278 /**************************************************************************************************/
3280 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){
3284 m->openOutputFile(fastaFileName, fastaFile);
3286 vector<string> names(numOTUs, "");
3288 for(int i=0;i<numOTUs;i++){
3290 if (m->control_pressed) { break; }
3292 int index = centroids[i];
3294 if(otuCounts[i] > 0){
3295 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
3299 for(int j=0;j<numFlowCells;j++){
3301 char base = flowOrder[j % 4];
3302 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
3307 if (newSeq.length() >= 4) { fastaFile << newSeq.substr(4) << endl; }
3308 else { fastaFile << "NNNN" << endl; }
3313 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
3315 if(thisCompositeFASTAFileName != ""){
3316 m->appendFiles(fastaFileName, thisCompositeFASTAFileName);
3319 catch(exception& e) {
3320 m->errorOut(e, "ShhherCommand", "writeSequences");
3325 /**************************************************************************************************/
3327 void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string nameFileName, vector<int> otuCounts, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& nSeqsPerOTU){
3331 m->openOutputFile(nameFileName, nameFile);
3333 for(int i=0;i<numOTUs;i++){
3335 if (m->control_pressed) { break; }
3337 if(otuCounts[i] > 0){
3338 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
3340 for(int j=1;j<nSeqsPerOTU[i];j++){
3341 nameFile << ',' << seqNameVector[aaI[i][j]];
3348 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
3351 if(thisCompositeNamesFileName != ""){
3352 m->appendFiles(nameFileName, thisCompositeNamesFileName);
3355 catch(exception& e) {
3356 m->errorOut(e, "ShhherCommand", "writeNames");
3361 /**************************************************************************************************/
3363 void ShhherCommand::writeGroups(string groupFileName, string fileRoot, int numSeqs, vector<string>& seqNameVector){
3366 m->openOutputFile(groupFileName, groupFile);
3368 for(int i=0;i<numSeqs;i++){
3369 if (m->control_pressed) { break; }
3370 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
3373 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
3376 catch(exception& e) {
3377 m->errorOut(e, "ShhherCommand", "writeGroups");
3382 /**************************************************************************************************/
3384 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){
3386 ofstream otuCountsFile;
3387 m->openOutputFile(otuCountsFileName, otuCountsFile);
3389 string bases = flowOrder;
3391 for(int i=0;i<numOTUs;i++){
3393 if (m->control_pressed) {
3396 //output the translated version of the centroid sequence for the otu
3397 if(otuCounts[i] > 0){
3398 int index = centroids[i];
3400 otuCountsFile << "ideal\t";
3401 for(int j=8;j<numFlowCells;j++){
3402 char base = bases[j % 4];
3403 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
3404 otuCountsFile << base;
3407 otuCountsFile << endl;
3409 for(int j=0;j<nSeqsPerOTU[i];j++){
3410 int sequence = aaI[i][j];
3411 otuCountsFile << seqNameVector[sequence] << '\t';
3415 for(int k=0;k<lengths[sequence];k++){
3416 char base = bases[k % 4];
3417 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
3419 for(int s=0;s<freq;s++){
3421 //otuCountsFile << base;
3425 if (newSeq.length() >= 4) { otuCountsFile << newSeq.substr(4) << endl; }
3426 else { otuCountsFile << "NNNN" << endl; }
3428 otuCountsFile << endl;
3431 otuCountsFile.close();
3432 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
3435 catch(exception& e) {
3436 m->errorOut(e, "ShhherCommand", "writeClusters");
3441 /**************************************************************************************************/
3443 void ShhherCommand::getSingleLookUp(){
3445 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
3446 singleLookUp.assign(HOMOPS * NUMBINS, 0);
3449 ifstream lookUpFile;
3450 m->openInputFile(lookupFileName, lookUpFile);
3452 for(int i=0;i<HOMOPS;i++){
3454 if (m->control_pressed) { break; }
3457 lookUpFile >> logFracFreq;
3459 for(int j=0;j<NUMBINS;j++) {
3460 lookUpFile >> singleLookUp[index];
3466 catch(exception& e) {
3467 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
3472 /**************************************************************************************************/
3474 void ShhherCommand::getJointLookUp(){
3477 // the most likely joint probability (-log) that two intenities have the same polymer length
3478 jointLookUp.resize(NUMBINS * NUMBINS, 0);
3480 for(int i=0;i<NUMBINS;i++){
3482 if (m->control_pressed) { break; }
3484 for(int j=0;j<NUMBINS;j++){
3486 double minSum = 100000000;
3488 for(int k=0;k<HOMOPS;k++){
3489 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
3491 if(sum < minSum) { minSum = sum; }
3493 jointLookUp[i * NUMBINS + j] = minSum;
3497 catch(exception& e) {
3498 m->errorOut(e, "ShhherCommand", "getJointLookUp");
3503 /**************************************************************************************************/
3505 double ShhherCommand::getProbIntensity(int intIntensity){
3508 double minNegLogProb = 100000000;
3511 for(int i=0;i<HOMOPS;i++){//loop signal strength
3513 if (m->control_pressed) { break; }
3515 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
3516 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
3519 return minNegLogProb;
3521 catch(exception& e) {
3522 m->errorOut(e, "ShhherCommand", "getProbIntensity");