5 * Created by Pat Schloss on 12/27/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "shhhercommand.h"
12 //**********************************************************************************************************************
13 vector<string> ShhherCommand::setParameters(){
15 CommandParameter pflow("flow", "InputTypes", "", "", "none", "fileflow", "none","fasta-name-group-counts-qfile",false,false,true); parameters.push_back(pflow);
16 CommandParameter pfile("file", "InputTypes", "", "", "none", "fileflow", "none","fasta-name-group-counts-qfile",false,false,true); parameters.push_back(pfile);
17 CommandParameter plookup("lookup", "InputTypes", "", "", "none", "none", "none","",false,false,true); parameters.push_back(plookup);
18 CommandParameter pcutoff("cutoff", "Number", "", "0.01", "", "", "","",false,false); parameters.push_back(pcutoff);
19 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
20 CommandParameter pmaxiter("maxiter", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(pmaxiter);
21 CommandParameter plarge("large", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(plarge);
22 CommandParameter psigma("sigma", "Number", "", "60", "", "", "","",false,false); parameters.push_back(psigma);
23 CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "","",false,false); parameters.push_back(pmindelta);
24 CommandParameter porder("order", "Multiple", "A-B-I", "A", "", "", "","",false,false, true); parameters.push_back(porder); CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
25 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
27 vector<string> myArray;
28 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
32 m->errorOut(e, "ShhherCommand", "setParameters");
36 //**********************************************************************************************************************
37 string ShhherCommand::getHelpString(){
39 string helpString = "";
40 helpString += "The shhh.flows command reads a file containing flowgrams and creates a file of corrected sequences.\n";
41 helpString += "The shhh.flows command parameters are flow, file, lookup, cutoff, processors, large, maxiter, sigma, mindelta and order.\n";
42 helpString += "The flow parameter is used to input your flow file.\n";
43 helpString += "The file parameter is used to input the *flow.files file created by trim.flows.\n";
44 helpString += "The lookup parameter is used specify the lookup file you would like to use. http://www.mothur.org/wiki/Lookup_files.\n";
45 helpString += "The order parameter options are A, B or I. Default=A. A = TACG and B = TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC and I = TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC.\n";
49 m->errorOut(e, "ShhherCommand", "getHelpString");
53 //**********************************************************************************************************************
54 string ShhherCommand::getOutputPattern(string type) {
58 if (type == "fasta") { pattern = "[filename],shhh.fasta"; }
59 else if (type == "name") { pattern = "[filename],shhh.names"; }
60 else if (type == "group") { pattern = "[filename],shhh.groups"; }
61 else if (type == "counts") { pattern = "[filename],shhh.counts"; }
62 else if (type == "qfile") { pattern = "[filename],shhh.qual"; }
63 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
68 m->errorOut(e, "ShhherCommand", "getOutputPattern");
72 //**********************************************************************************************************************
74 ShhherCommand::ShhherCommand(){
76 abort = true; calledHelp = true;
79 //initialize outputTypes
80 vector<string> tempOutNames;
81 outputTypes["fasta"] = tempOutNames;
82 outputTypes["name"] = tempOutNames;
83 outputTypes["group"] = tempOutNames;
84 outputTypes["counts"] = tempOutNames;
85 outputTypes["qfile"] = tempOutNames;
89 m->errorOut(e, "ShhherCommand", "ShhherCommand");
94 //**********************************************************************************************************************
96 ShhherCommand::ShhherCommand(string option) {
100 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
101 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
105 abort = false; calledHelp = false;
107 //allow user to run help
108 if(option == "help") { help(); abort = true; calledHelp = true; }
109 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
112 vector<string> myArray = setParameters();
114 OptionParser parser(option);
115 map<string,string> parameters = parser.getParameters();
117 ValidParameters validParameter;
118 map<string,string>::iterator it;
120 //check to make sure all parameters are valid for command
121 for (it = parameters.begin(); it != parameters.end(); it++) {
122 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
125 //initialize outputTypes
126 vector<string> tempOutNames;
127 outputTypes["fasta"] = tempOutNames;
128 outputTypes["name"] = tempOutNames;
129 outputTypes["group"] = tempOutNames;
130 outputTypes["counts"] = tempOutNames;
131 outputTypes["qfile"] = tempOutNames;
134 //if the user changes the input directory command factory will send this info to us in the output parameter
135 string inputDir = validParameter.validFile(parameters, "inputdir", false);
136 if (inputDir == "not found"){ inputDir = ""; }
139 it = parameters.find("flow");
140 //user has given a template file
141 if(it != parameters.end()){
142 path = m->hasPath(it->second);
143 //if the user has not given a path then, add inputdir. else leave path alone.
144 if (path == "") { parameters["flow"] = inputDir + it->second; }
147 it = parameters.find("lookup");
148 //user has given a template file
149 if(it != parameters.end()){
150 path = m->hasPath(it->second);
151 //if the user has not given a path then, add inputdir. else leave path alone.
152 if (path == "") { parameters["lookup"] = inputDir + it->second; }
155 it = parameters.find("file");
156 //user has given a template file
157 if(it != parameters.end()){
158 path = m->hasPath(it->second);
159 //if the user has not given a path then, add inputdir. else leave path alone.
160 if (path == "") { parameters["file"] = inputDir + it->second; }
164 //if the user changes the output directory command factory will send this info to us in the output parameter
165 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
167 //check for required parameters
168 flowFileName = validParameter.validFile(parameters, "flow", true);
169 flowFilesFileName = validParameter.validFile(parameters, "file", true);
170 if (flowFileName == "not found" && flowFilesFileName == "not found") {
171 m->mothurOut("values for either flow or file must be provided for the shhh.flows command.");
172 m->mothurOutEndLine();
175 else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
177 if(flowFileName != "not found"){
178 compositeFASTAFileName = "";
179 compositeNamesFileName = "";
184 string thisoutputDir = outputDir;
185 if (outputDir == "") { thisoutputDir = m->hasPath(flowFilesFileName); } //if user entered a file with a path then preserve it
187 //we want to rip off .files, and also .flow if its there
188 string fileroot = m->getRootName(m->getSimpleName(flowFilesFileName));
189 if (fileroot[fileroot.length()-1] == '.') { fileroot = fileroot.substr(0, fileroot.length()-1); } //rip off dot
190 string extension = m->getExtension(fileroot);
191 if (extension == ".flow") { fileroot = m->getRootName(fileroot); }
192 else { fileroot += "."; } //add back if needed
194 compositeFASTAFileName = thisoutputDir + fileroot + "shhh.fasta";
195 m->openOutputFile(compositeFASTAFileName, temp);
198 compositeNamesFileName = thisoutputDir + fileroot + "shhh.names";
199 m->openOutputFile(compositeNamesFileName, temp);
203 if(flowFilesFileName != "not found"){
206 ifstream flowFilesFile;
207 m->openInputFile(flowFilesFileName, flowFilesFile);
208 while(flowFilesFile){
209 fName = m->getline(flowFilesFile);
211 //test if file is valid
213 int ableToOpen = m->openInputFile(fName, in, "noerror");
215 if (ableToOpen == 1) {
216 if (inputDir != "") { //default path is set
217 string tryPath = inputDir + fName;
218 m->mothurOut("Unable to open " + fName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
220 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
226 if (ableToOpen == 1) {
227 if (m->getDefaultPath() != "") { //default path is set
228 string tryPath = m->getDefaultPath() + m->getSimpleName(fName);
229 m->mothurOut("Unable to open " + fName + ". Trying default " + tryPath); m->mothurOutEndLine();
231 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
237 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
238 if (ableToOpen == 1) {
239 string exepath = m->argv;
240 string tempPath = exepath;
241 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
242 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
244 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(fName);
245 m->mothurOut("Unable to open " + fName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
247 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
252 if (ableToOpen == 1) { m->mothurOut("Unable to open " + fName + ". Disregarding. "); m->mothurOutEndLine(); }
253 else { flowFileVector.push_back(fName); }
254 m->gobble(flowFilesFile);
256 flowFilesFile.close();
257 if (flowFileVector.size() == 0) { m->mothurOut("[ERROR]: no valid files."); m->mothurOutEndLine(); abort = true; }
260 if (outputDir == "") { outputDir = m->hasPath(flowFileName); }
261 flowFileVector.push_back(flowFileName);
264 //check for optional parameter and set defaults
265 // ...at some point should added some additional type checking...
267 temp = validParameter.validFile(parameters, "lookup", true);
268 if (temp == "not found") {
269 string path = m->argv;
270 string tempPath = path;
271 for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
272 path = path.substr(0, (tempPath.find_last_of('m')));
274 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
275 path += "lookupFiles/";
277 path += "lookupFiles\\";
279 lookupFileName = m->getFullPathName(path) + "LookUp_Titanium.pat";
283 ableToOpen = m->openInputFile(lookupFileName, in, "noerror");
286 //if you can't open it, try input location
287 if (ableToOpen == 1) {
288 if (inputDir != "") { //default path is set
289 string tryPath = inputDir + m->getSimpleName(lookupFileName);
290 m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
292 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
294 lookupFileName = tryPath;
298 //if you can't open it, try default location
299 if (ableToOpen == 1) {
300 if (m->getDefaultPath() != "") { //default path is set
301 string tryPath = m->getDefaultPath() + m->getSimpleName(lookupFileName);
302 m->mothurOut("Unable to open " + lookupFileName + ". Trying default " + tryPath); m->mothurOutEndLine();
304 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
306 lookupFileName = tryPath;
310 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
311 if (ableToOpen == 1) {
312 string exepath = m->argv;
313 string tempPath = exepath;
314 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
315 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
317 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
318 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
320 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
322 lookupFileName = tryPath;
325 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
327 else if(temp == "not open") {
329 lookupFileName = validParameter.validFile(parameters, "lookup", false);
331 //if you can't open it its not inputDir, try mothur excutable location
332 string exepath = m->argv;
333 string tempPath = exepath;
334 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
335 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
337 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
338 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
340 int ableToOpen = m->openInputFile(tryPath, in2, "noerror");
342 lookupFileName = tryPath;
344 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
345 }else { lookupFileName = temp; }
347 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
348 m->setProcessors(temp);
349 m->mothurConvert(temp, processors);
351 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
352 m->mothurConvert(temp, cutoff);
354 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
355 m->mothurConvert(temp, minDelta);
357 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
358 m->mothurConvert(temp, maxIters);
360 temp = validParameter.validFile(parameters, "large", false); if (temp == "not found"){ temp = "0"; }
361 m->mothurConvert(temp, largeSize);
362 if (largeSize != 0) { large = true; }
363 else { large = false; }
364 if (largeSize < 0) { m->mothurOut("The value of the large cannot be negative.\n"); }
367 if (large) { m->mothurOut("The large parameter is not available with the MPI-Enabled version.\n"); large=false; }
371 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
372 m->mothurConvert(temp, sigma);
374 temp = validParameter.validFile(parameters, "order", false); if (temp == "not found"){ temp = "A"; }
375 if (temp.length() > 1) { m->mothurOut("[ERROR]: " + temp + " is not a valid option for order. order options are A, B, or I. A = TACG, B = TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC, and I = TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC.\n"); abort=true;
378 if (toupper(temp[0]) == 'A') { flowOrder = "TACG"; }
379 else if(toupper(temp[0]) == 'B'){
380 flowOrder = "TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC"; }
381 else if(toupper(temp[0]) == 'I'){
382 flowOrder = "TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC"; }
384 m->mothurOut("[ERROR]: " + temp + " is not a valid option for order. order options are A, B, or I. A = TACG, B = TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC, and I = TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC.\n"); abort=true;
394 catch(exception& e) {
395 m->errorOut(e, "ShhherCommand", "ShhherCommand");
399 //**********************************************************************************************************************
401 int ShhherCommand::execute(){
403 if (abort == true) { if (calledHelp) { return 0; } return 2; }
410 for(int i=1;i<ncpus;i++){
411 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
413 if(abort == 1){ return 0; }
417 m->mothurOut("\nGetting preliminary data...\n");
418 getSingleLookUp(); if (m->control_pressed) { return 0; }
419 getJointLookUp(); if (m->control_pressed) { return 0; }
421 vector<string> flowFileVector;
422 if(flowFilesFileName != "not found"){
425 ifstream flowFilesFile;
426 m->openInputFile(flowFilesFileName, flowFilesFile);
427 while(flowFilesFile){
428 fName = m->getline(flowFilesFile);
429 flowFileVector.push_back(fName);
430 m->gobble(flowFilesFile);
434 flowFileVector.push_back(flowFileName);
437 int numFiles = flowFileVector.size();
439 for(int i=1;i<ncpus;i++){
440 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
443 for(int i=0;i<numFiles;i++){
445 if (m->control_pressed) { break; }
447 double begClock = clock();
448 unsigned long long begTime = time(NULL);
450 flowFileName = flowFileVector[i];
452 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
453 m->mothurOut("Reading flowgrams...\n");
456 if (m->control_pressed) { break; }
458 m->mothurOut("Identifying unique flowgrams...\n");
461 if (m->control_pressed) { break; }
463 m->mothurOut("Calculating distances between flowgrams...\n");
465 strcpy(fileName, flowFileName.c_str());
467 for(int i=1;i<ncpus;i++){
468 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
470 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
471 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
472 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
473 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
474 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
475 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
476 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
477 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
478 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
479 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
482 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
484 if (m->control_pressed) { break; }
487 for(int i=1;i<ncpus;i++){
488 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
490 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
491 m->mothurRemove((distFileName + ".temp." + toString(i)));
494 string namesFileName = createNamesFile();
496 if (m->control_pressed) { break; }
498 m->mothurOut("\nClustering flowgrams...\n");
499 string listFileName = cluster(distFileName, namesFileName);
501 if (m->control_pressed) { break; }
505 getOTUData(listFileName);
507 m->mothurRemove(distFileName);
508 m->mothurRemove(namesFileName);
509 m->mothurRemove(listFileName);
511 if (m->control_pressed) { break; }
515 if (m->control_pressed) { break; }
518 for(int i=1;i<ncpus;i++){
519 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
520 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
521 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
522 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
525 if (m->control_pressed) { break; }
530 int numOTUsOnCPU = numOTUs / ncpus;
531 int numSeqsOnCPU = numSeqs / ncpus;
532 m->mothurOut("\nDenoising flowgrams...\n");
533 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
535 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
537 double cycClock = clock();
538 unsigned long long cycTime = time(NULL);
541 if (m->control_pressed) { break; }
543 int total = singleTau.size();
544 for(int i=1;i<ncpus;i++){
545 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
546 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
547 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
549 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
550 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
551 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
552 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
553 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
556 calcCentroidsDriver(0, numOTUsOnCPU);
558 for(int i=1;i<ncpus;i++){
559 int otuStart = i * numOTUs / ncpus;
560 int otuStop = (i + 1) * numOTUs / ncpus;
562 vector<int> tempCentroids(numOTUs, 0);
563 vector<short> tempChange(numOTUs, 0);
565 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
566 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
568 for(int j=otuStart;j<otuStop;j++){
569 centroids[j] = tempCentroids[j];
570 change[j] = tempChange[j];
574 maxDelta = getNewWeights(); if (m->control_pressed) { break; }
575 double nLL = getLikelihood(); if (m->control_pressed) { break; }
576 checkCentroids(); if (m->control_pressed) { break; }
578 for(int i=1;i<ncpus;i++){
579 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
580 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
581 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
584 calcNewDistancesParent(0, numSeqsOnCPU);
586 total = singleTau.size();
588 for(int i=1;i<ncpus;i++){
590 int seqStart = i * numSeqs / ncpus;
591 int seqStop = (i + 1) * numSeqs / ncpus;
593 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
595 vector<int> childSeqIndex(childTotal, 0);
596 vector<double> childSingleTau(childTotal, 0);
597 vector<double> childDist(numSeqs * numOTUs, 0);
598 vector<int> otuIndex(childTotal, 0);
600 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
601 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
602 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
603 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
605 int oldTotal = total;
607 singleTau.resize(total, 0);
608 seqIndex.resize(total, 0);
609 seqNumber.resize(total, 0);
613 for(int j=oldTotal;j<total;j++){
614 int otuI = otuIndex[childIndex];
615 int seqI = childSeqIndex[childIndex];
617 singleTau[j] = childSingleTau[childIndex];
619 aaP[otuI][nSeqsPerOTU[otuI]] = j;
620 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
625 int index = seqStart * numOTUs;
626 for(int j=seqStart;j<seqStop;j++){
627 for(int k=0;k<numOTUs;k++){
628 dist[index] = childDist[index];
636 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
638 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
640 for(int i=1;i<ncpus;i++){
641 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
646 for(int i=1;i<ncpus;i++){
647 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
653 if (m->control_pressed) { break; }
655 m->mothurOut("\nFinalizing...\n");
658 if (m->control_pressed) { break; }
662 vector<int> otuCounts(numOTUs, 0);
663 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
664 calcCentroidsDriver(0, numOTUs);
666 if (m->control_pressed) { break; }
668 writeQualities(otuCounts); if (m->control_pressed) { break; }
669 writeSequences(otuCounts); if (m->control_pressed) { break; }
670 writeNames(otuCounts); if (m->control_pressed) { break; }
671 writeClusters(otuCounts); if (m->control_pressed) { break; }
672 writeGroups(); if (m->control_pressed) { break; }
675 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
681 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
682 if(abort){ return 0; }
685 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
687 for(int i=0;i<numFiles;i++){
689 if (m->control_pressed) { break; }
691 //Now into the pyrodist part
695 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
696 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
697 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
698 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
700 flowDataIntI.resize(numSeqs * numFlowCells);
701 flowDataPrI.resize(numSeqs * numFlowCells);
702 mapUniqueToSeq.resize(numSeqs);
703 mapSeqToUnique.resize(numSeqs);
704 lengths.resize(numSeqs);
705 jointLookUp.resize(NUMBINS * NUMBINS);
707 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
708 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
709 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
710 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
711 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
712 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
713 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
715 flowFileName = string(fileName);
716 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
717 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
719 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
721 if (m->control_pressed) { break; }
724 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
726 //Now into the pyronoise part
727 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
729 singleLookUp.resize(HOMOPS * NUMBINS);
730 uniqueFlowgrams.resize(numUniques * numFlowCells);
731 weight.resize(numOTUs);
732 centroids.resize(numOTUs);
733 change.resize(numOTUs);
734 dist.assign(numOTUs * numSeqs, 0);
735 nSeqsPerOTU.resize(numOTUs);
736 cumNumSeqs.resize(numOTUs);
738 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
739 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
740 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
742 int startOTU = pid * numOTUs / ncpus;
743 int endOTU = (pid + 1) * numOTUs / ncpus;
745 int startSeq = pid * numSeqs / ncpus;
746 int endSeq = (pid + 1) * numSeqs /ncpus;
752 if (m->control_pressed) { break; }
754 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
755 singleTau.assign(total, 0.0000);
756 seqNumber.assign(total, 0);
757 seqIndex.assign(total, 0);
759 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
760 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
761 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
762 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
763 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
764 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
765 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
767 calcCentroidsDriver(startOTU, endOTU);
769 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
770 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
772 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
773 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
774 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
776 vector<int> otuIndex(total, 0);
777 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
778 total = otuIndex.size();
780 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
781 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
782 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
783 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
784 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
786 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
791 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
793 MPI_Barrier(MPI_COMM_WORLD);
796 if(compositeFASTAFileName != ""){
797 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
798 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
801 m->mothurOutEndLine();
802 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
803 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
804 m->mothurOutEndLine();
809 catch(exception& e) {
810 m->errorOut(e, "ShhherCommand", "execute");
814 /**************************************************************************************************/
815 string ShhherCommand::createNamesFile(){
818 vector<string> duplicateNames(numUniques, "");
819 for(int i=0;i<numSeqs;i++){
820 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
822 map<string, string> variables;
823 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName));
824 string nameFileName = getOutputFileName("name",variables);
827 m->openOutputFile(nameFileName, nameFile);
829 for(int i=0;i<numUniques;i++){
831 if (m->control_pressed) { break; }
833 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
834 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
840 catch(exception& e) {
841 m->errorOut(e, "ShhherCommand", "createNamesFile");
845 /**************************************************************************************************/
847 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
849 ostringstream outStream;
850 outStream.setf(ios::fixed, ios::floatfield);
851 outStream.setf(ios::dec, ios::basefield);
852 outStream.setf(ios::showpoint);
853 outStream.precision(6);
855 int begTime = time(NULL);
856 double begClock = clock();
858 for(int i=startSeq;i<stopSeq;i++){
860 if (m->control_pressed) { break; }
862 for(int j=0;j<i;j++){
863 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
865 if(flowDistance < 1e-6){
866 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
868 else if(flowDistance <= cutoff){
869 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
873 m->mothurOutJustToScreen(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
877 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
878 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
880 if (m->control_pressed) { return fDistFileName; }
882 m->mothurOutJustToScreen(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
884 ofstream distFile(fDistFileName.c_str());
885 distFile << outStream.str();
888 return fDistFileName;
890 catch(exception& e) {
891 m->errorOut(e, "ShhherCommand", "flowDistMPI");
895 /**************************************************************************************************/
897 void ShhherCommand::getOTUData(string listFileName){
901 m->openInputFile(listFileName, listFile);
904 listFile >> label >> numOTUs;
906 otuData.assign(numSeqs, 0);
907 cumNumSeqs.assign(numOTUs, 0);
908 nSeqsPerOTU.assign(numOTUs, 0);
909 aaP.clear();aaP.resize(numOTUs);
915 string singleOTU = "";
917 for(int i=0;i<numOTUs;i++){
919 if (m->control_pressed) { break; }
921 listFile >> singleOTU;
923 istringstream otuString(singleOTU);
929 for(int j=0;j<singleOTU.length();j++){
930 char letter = otuString.get();
936 map<string,int>::iterator nmIt = nameMap.find(seqName);
937 int index = nmIt->second;
943 aaP[i].push_back(index);
948 map<string,int>::iterator nmIt = nameMap.find(seqName);
950 int index = nmIt->second;
955 aaP[i].push_back(index);
960 sort(aaP[i].begin(), aaP[i].end());
961 for(int j=0;j<nSeqsPerOTU[i];j++){
962 seqNumber.push_back(aaP[i][j]);
964 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
971 for(int i=1;i<numOTUs;i++){
972 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
975 seqIndex = seqNumber;
980 catch(exception& e) {
981 m->errorOut(e, "ShhherCommand", "getOTUData");
986 /**************************************************************************************************/
988 void ShhherCommand::initPyroCluster(){
990 if (numOTUs < processors) { processors = 1; }
992 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
994 dist.assign(numSeqs * numOTUs, 0);
995 change.assign(numOTUs, 1);
996 centroids.assign(numOTUs, -1);
997 weight.assign(numOTUs, 0);
998 singleTau.assign(numSeqs, 1.0);
1000 nSeqsBreaks.assign(processors+1, 0);
1001 nOTUsBreaks.assign(processors+1, 0);
1003 if (m->debug) { m->mothurOut("[DEBUG]: made it through the memory allocation.\n"); }
1006 for(int i=0;i<processors;i++){
1007 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
1008 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
1010 nSeqsBreaks[processors] = numSeqs;
1011 nOTUsBreaks[processors] = numOTUs;
1013 catch(exception& e) {
1014 m->errorOut(e, "ShhherCommand", "initPyroCluster");
1019 /**************************************************************************************************/
1021 void ShhherCommand::fill(){
1024 for(int i=0;i<numOTUs;i++){
1026 if (m->control_pressed) { break; }
1028 cumNumSeqs[i] = index;
1029 for(int j=0;j<nSeqsPerOTU[i];j++){
1030 seqNumber[index] = aaP[i][j];
1031 seqIndex[index] = aaI[i][j];
1037 catch(exception& e) {
1038 m->errorOut(e, "ShhherCommand", "fill");
1043 /**************************************************************************************************/
1045 void ShhherCommand::getFlowData(){
1048 m->openInputFile(flowFileName, flowFile);
1051 seqNameVector.clear();
1053 flowDataIntI.clear();
1057 int currentNumFlowCells;
1062 flowFile >> numFlowTest;
1064 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); }
1065 else { convert(numFlowTest, numFlowCells); }
1067 int index = 0;//pcluster
1068 while(!flowFile.eof()){
1070 if (m->control_pressed) { break; }
1072 flowFile >> seqName >> currentNumFlowCells;
1073 lengths.push_back(currentNumFlowCells);
1075 seqNameVector.push_back(seqName);
1076 nameMap[seqName] = index++;//pcluster
1078 for(int i=0;i<numFlowCells;i++){
1079 flowFile >> intensity;
1080 if(intensity > 9.99) { intensity = 9.99; }
1081 int intI = int(100 * intensity + 0.0001);
1082 flowDataIntI.push_back(intI);
1084 m->gobble(flowFile);
1088 numSeqs = seqNameVector.size();
1090 for(int i=0;i<numSeqs;i++){
1092 if (m->control_pressed) { break; }
1094 int iNumFlowCells = i * numFlowCells;
1095 for(int j=lengths[i];j<numFlowCells;j++){
1096 flowDataIntI[iNumFlowCells + j] = 0;
1101 catch(exception& e) {
1102 m->errorOut(e, "ShhherCommand", "getFlowData");
1106 /**************************************************************************************************/
1107 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1110 vector<double> newTau(numOTUs,0);
1111 vector<double> norms(numSeqs, 0);
1116 for(int i=startSeq;i<stopSeq;i++){
1118 if (m->control_pressed) { break; }
1120 double offset = 1e8;
1121 int indexOffset = i * numOTUs;
1123 for(int j=0;j<numOTUs;j++){
1125 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1126 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1128 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1129 offset = dist[indexOffset + j];
1133 for(int j=0;j<numOTUs;j++){
1134 if(weight[j] > MIN_WEIGHT){
1135 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1136 norms[i] += newTau[j];
1143 for(int j=0;j<numOTUs;j++){
1145 newTau[j] /= norms[i];
1147 if(newTau[j] > MIN_TAU){
1148 otuIndex.push_back(j);
1149 seqIndex.push_back(i);
1150 singleTau.push_back(newTau[j]);
1156 catch(exception& e) {
1157 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1162 /**************************************************************************************************/
1164 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1169 vector<double> newTau(numOTUs,0);
1170 vector<double> norms(numSeqs, 0);
1171 nSeqsPerOTU.assign(numOTUs, 0);
1173 for(int i=startSeq;i<stopSeq;i++){
1175 if (m->control_pressed) { break; }
1177 int indexOffset = i * numOTUs;
1179 double offset = 1e8;
1181 for(int j=0;j<numOTUs;j++){
1183 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1184 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1187 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1188 offset = dist[indexOffset + j];
1192 for(int j=0;j<numOTUs;j++){
1193 if(weight[j] > MIN_WEIGHT){
1194 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1195 norms[i] += newTau[j];
1202 for(int j=0;j<numOTUs;j++){
1203 newTau[j] /= norms[i];
1206 for(int j=0;j<numOTUs;j++){
1207 if(newTau[j] > MIN_TAU){
1209 int oldTotal = total;
1213 singleTau.resize(total, 0);
1214 seqNumber.resize(total, 0);
1215 seqIndex.resize(total, 0);
1217 singleTau[oldTotal] = newTau[j];
1219 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1220 aaI[j][nSeqsPerOTU[j]] = i;
1228 catch(exception& e) {
1229 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1234 /**************************************************************************************************/
1236 void ShhherCommand::setOTUs(){
1239 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1241 for(int i=0;i<numOTUs;i++){
1243 if (m->control_pressed) { break; }
1245 for(int j=0;j<nSeqsPerOTU[i];j++){
1246 int index = cumNumSeqs[i] + j;
1247 double tauValue = singleTau[seqNumber[index]];
1248 int sIndex = seqIndex[index];
1249 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1253 for(int i=0;i<numSeqs;i++){
1254 double maxTau = -1.0000;
1256 for(int j=0;j<numOTUs;j++){
1257 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1258 maxTau = bigTauMatrix[i * numOTUs + j];
1263 otuData[i] = maxOTU;
1266 nSeqsPerOTU.assign(numOTUs, 0);
1268 for(int i=0;i<numSeqs;i++){
1269 int index = otuData[i];
1271 singleTau[i] = 1.0000;
1274 aaP[index][nSeqsPerOTU[index]] = i;
1275 aaI[index][nSeqsPerOTU[index]] = i;
1277 nSeqsPerOTU[index]++;
1281 catch(exception& e) {
1282 m->errorOut(e, "ShhherCommand", "setOTUs");
1287 /**************************************************************************************************/
1289 void ShhherCommand::getUniques(){
1294 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
1295 uniqueCount.assign(numSeqs, 0); // anWeights
1296 uniqueLengths.assign(numSeqs, 0);
1297 mapSeqToUnique.assign(numSeqs, -1);
1298 mapUniqueToSeq.assign(numSeqs, -1);
1300 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
1302 for(int i=0;i<numSeqs;i++){
1304 if (m->control_pressed) { break; }
1308 vector<short> current(numFlowCells);
1309 for(int j=0;j<numFlowCells;j++){
1310 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
1313 for(int j=0;j<numUniques;j++){
1314 int offset = j * numFlowCells;
1318 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
1319 else { shorterLength = uniqueLengths[j]; }
1321 for(int k=0;k<shorterLength;k++){
1322 if(current[k] != uniqueFlowgrams[offset + k]){
1329 mapSeqToUnique[i] = j;
1332 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
1338 if(index == numUniques){
1339 uniqueLengths[numUniques] = lengths[i];
1340 uniqueCount[numUniques] = 1;
1341 mapSeqToUnique[i] = numUniques;//anMap
1342 mapUniqueToSeq[numUniques] = i;//anF
1344 for(int k=0;k<numFlowCells;k++){
1345 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
1346 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
1352 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
1353 uniqueLengths.resize(numUniques);
1355 flowDataPrI.resize(numSeqs * numFlowCells, 0);
1356 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
1358 catch(exception& e) {
1359 m->errorOut(e, "ShhherCommand", "getUniques");
1364 /**************************************************************************************************/
1366 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
1368 int minLength = lengths[mapSeqToUnique[seqA]];
1369 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
1371 int ANumFlowCells = seqA * numFlowCells;
1372 int BNumFlowCells = seqB * numFlowCells;
1376 for(int i=0;i<minLength;i++){
1378 if (m->control_pressed) { break; }
1380 int flowAIntI = flowDataIntI[ANumFlowCells + i];
1381 float flowAPrI = flowDataPrI[ANumFlowCells + i];
1383 int flowBIntI = flowDataIntI[BNumFlowCells + i];
1384 float flowBPrI = flowDataPrI[BNumFlowCells + i];
1385 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
1388 dist /= (float) minLength;
1391 catch(exception& e) {
1392 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
1397 //**********************************************************************************************************************/
1399 string ShhherCommand::cluster(string distFileName, string namesFileName){
1402 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1403 read->setCutoff(cutoff);
1405 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1406 clusterNameMap->readMap();
1407 read->read(clusterNameMap);
1409 ListVector* list = read->getListVector();
1410 SparseDistanceMatrix* matrix = read->getDMatrix();
1413 delete clusterNameMap;
1415 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1417 float adjust = -1.0;
1418 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest", adjust);
1419 string tag = cluster->getTag();
1421 double clusterCutoff = cutoff;
1422 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1424 if (m->control_pressed) { break; }
1426 cluster->update(clusterCutoff);
1429 list->setLabel(toString(cutoff));
1431 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
1433 m->openOutputFile(listFileName, listFile);
1434 list->print(listFile);
1437 delete matrix; delete cluster; delete rabund; delete list;
1439 return listFileName;
1441 catch(exception& e) {
1442 m->errorOut(e, "ShhherCommand", "cluster");
1447 /**************************************************************************************************/
1449 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1451 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1456 for(int i=start;i<finish;i++){
1458 if (m->control_pressed) { break; }
1462 int minFlowGram = 100000000;
1463 double minFlowValue = 1e8;
1464 change[i] = 0; //FALSE
1466 for(int j=0;j<nSeqsPerOTU[i];j++){
1467 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1470 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1471 vector<double> adF(nSeqsPerOTU[i]);
1472 vector<int> anL(nSeqsPerOTU[i]);
1474 for(int j=0;j<nSeqsPerOTU[i];j++){
1475 int index = cumNumSeqs[i] + j;
1476 int nI = seqIndex[index];
1477 int nIU = mapSeqToUnique[nI];
1480 for(k=0;k<position;k++){
1486 anL[position] = nIU;
1487 adF[position] = 0.0000;
1492 for(int j=0;j<nSeqsPerOTU[i];j++){
1493 int index = cumNumSeqs[i] + j;
1494 int nI = seqIndex[index];
1496 double tauValue = singleTau[seqNumber[index]];
1498 for(int k=0;k<position;k++){
1499 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1500 adF[k] += dist * tauValue;
1504 for(int j=0;j<position;j++){
1505 if(adF[j] < minFlowValue){
1507 minFlowValue = adF[j];
1511 if(centroids[i] != anL[minFlowGram]){
1513 centroids[i] = anL[minFlowGram];
1516 else if(centroids[i] != -1){
1522 catch(exception& e) {
1523 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1528 /**************************************************************************************************/
1530 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1533 int flowAValue = cent * numFlowCells;
1534 int flowBValue = flow * numFlowCells;
1538 for(int i=0;i<length;i++){
1539 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1544 return dist / (double)length;
1546 catch(exception& e) {
1547 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1552 /**************************************************************************************************/
1554 double ShhherCommand::getNewWeights(){
1557 double maxChange = 0;
1559 for(int i=0;i<numOTUs;i++){
1561 if (m->control_pressed) { break; }
1563 double difference = weight[i];
1566 for(int j=0;j<nSeqsPerOTU[i];j++){
1567 int index = cumNumSeqs[i] + j;
1568 double tauValue = singleTau[seqNumber[index]];
1569 weight[i] += tauValue;
1572 difference = fabs(weight[i] - difference);
1573 if(difference > maxChange){ maxChange = difference; }
1577 catch(exception& e) {
1578 m->errorOut(e, "ShhherCommand", "getNewWeights");
1583 /**************************************************************************************************/
1585 double ShhherCommand::getLikelihood(){
1589 vector<long double> P(numSeqs, 0);
1592 for(int i=0;i<numOTUs;i++){
1593 if(weight[i] > MIN_WEIGHT){
1599 for(int i=0;i<numOTUs;i++){
1601 if (m->control_pressed) { break; }
1603 for(int j=0;j<nSeqsPerOTU[i];j++){
1604 int index = cumNumSeqs[i] + j;
1605 int nI = seqIndex[index];
1606 double singleDist = dist[seqNumber[index]];
1608 P[nI] += weight[i] * exp(-singleDist * sigma);
1612 for(int i=0;i<numSeqs;i++){
1613 if(P[i] == 0){ P[i] = DBL_EPSILON; }
1618 nLL = nLL -(double)numSeqs * log(sigma);
1622 catch(exception& e) {
1623 m->errorOut(e, "ShhherCommand", "getNewWeights");
1628 /**************************************************************************************************/
1630 void ShhherCommand::checkCentroids(){
1632 vector<int> unique(numOTUs, 1);
1634 for(int i=0;i<numOTUs;i++){
1635 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1640 for(int i=0;i<numOTUs;i++){
1642 if (m->control_pressed) { break; }
1645 for(int j=i+1;j<numOTUs;j++){
1648 if(centroids[j] == centroids[i]){
1652 weight[i] += weight[j];
1660 catch(exception& e) {
1661 m->errorOut(e, "ShhherCommand", "checkCentroids");
1665 /**************************************************************************************************/
1669 void ShhherCommand::writeQualities(vector<int> otuCounts){
1672 string thisOutputDir = outputDir;
1673 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1674 map<string, string> variables;
1675 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName));
1676 string qualityFileName = getOutputFileName("qfile",variables);
1678 ofstream qualityFile;
1679 m->openOutputFile(qualityFileName, qualityFile);
1681 qualityFile.setf(ios::fixed, ios::floatfield);
1682 qualityFile.setf(ios::showpoint);
1683 qualityFile << setprecision(6);
1685 vector<vector<int> > qualities(numOTUs);
1686 vector<double> pr(HOMOPS, 0);
1689 for(int i=0;i<numOTUs;i++){
1691 if (m->control_pressed) { break; }
1695 if(nSeqsPerOTU[i] > 0){
1697 while(index < numFlowCells){
1698 double maxPrValue = 1e8;
1699 short maxPrIndex = -1;
1700 double count = 0.0000;
1702 pr.assign(HOMOPS, 0);
1704 for(int j=0;j<nSeqsPerOTU[i];j++){
1705 int lIndex = cumNumSeqs[i] + j;
1706 double tauValue = singleTau[seqNumber[lIndex]];
1707 int sequenceIndex = aaI[i][j];
1708 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1712 for(int s=0;s<HOMOPS;s++){
1713 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1717 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1718 maxPrValue = pr[maxPrIndex];
1720 if(count > MIN_COUNT){
1722 double norm = 0.0000;
1724 for(int s=0;s<HOMOPS;s++){
1725 norm += exp(-(pr[s] - maxPrValue));
1728 for(int s=1;s<=maxPrIndex;s++){
1730 double temp = 0.0000;
1732 U += exp(-(pr[s-1]-maxPrValue))/norm;
1740 temp = floor(-10 * temp);
1741 value = (int)floor(temp);
1742 if(value > 100){ value = 100; }
1744 qualities[i].push_back((int)value);
1753 if(otuCounts[i] > 0){
1754 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1755 for (int j = 4; j < qualities[i].size(); j++) { qualityFile << qualities[i][j] << ' '; }
1756 qualityFile << endl;
1759 qualityFile.close();
1760 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
1763 catch(exception& e) {
1764 m->errorOut(e, "ShhherCommand", "writeQualities");
1769 /**************************************************************************************************/
1771 void ShhherCommand::writeSequences(vector<int> otuCounts){
1773 string thisOutputDir = outputDir;
1774 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1775 map<string, string> variables;
1776 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1777 string fastaFileName = getOutputFileName("fasta",variables);
1779 m->openOutputFile(fastaFileName, fastaFile);
1781 vector<string> names(numOTUs, "");
1783 for(int i=0;i<numOTUs;i++){
1785 if (m->control_pressed) { break; }
1787 int index = centroids[i];
1789 if(otuCounts[i] > 0){
1790 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
1794 for(int j=0;j<numFlowCells;j++){
1796 char base = flowOrder[j % flowOrder.length()];
1797 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
1802 fastaFile << newSeq.substr(4) << endl;
1807 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
1809 if(compositeFASTAFileName != ""){
1810 m->appendFiles(fastaFileName, compositeFASTAFileName);
1813 catch(exception& e) {
1814 m->errorOut(e, "ShhherCommand", "writeSequences");
1819 /**************************************************************************************************/
1821 void ShhherCommand::writeNames(vector<int> otuCounts){
1823 string thisOutputDir = outputDir;
1824 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1825 map<string, string> variables;
1826 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1827 string nameFileName = getOutputFileName("name",variables);
1829 m->openOutputFile(nameFileName, nameFile);
1831 for(int i=0;i<numOTUs;i++){
1833 if (m->control_pressed) { break; }
1835 if(otuCounts[i] > 0){
1836 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
1838 for(int j=1;j<nSeqsPerOTU[i];j++){
1839 nameFile << ',' << seqNameVector[aaI[i][j]];
1846 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
1849 if(compositeNamesFileName != ""){
1850 m->appendFiles(nameFileName, compositeNamesFileName);
1853 catch(exception& e) {
1854 m->errorOut(e, "ShhherCommand", "writeNames");
1859 /**************************************************************************************************/
1861 void ShhherCommand::writeGroups(){
1863 string thisOutputDir = outputDir;
1864 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1865 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
1866 int pos = fileRoot.find_first_of('.');
1867 string fileGroup = fileRoot;
1868 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
1869 map<string, string> variables;
1870 variables["[filename]"] = thisOutputDir + fileRoot;
1871 string groupFileName = getOutputFileName("group",variables);
1873 m->openOutputFile(groupFileName, groupFile);
1875 for(int i=0;i<numSeqs;i++){
1876 if (m->control_pressed) { break; }
1877 groupFile << seqNameVector[i] << '\t' << fileGroup << endl;
1880 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
1883 catch(exception& e) {
1884 m->errorOut(e, "ShhherCommand", "writeGroups");
1889 /**************************************************************************************************/
1891 void ShhherCommand::writeClusters(vector<int> otuCounts){
1893 string thisOutputDir = outputDir;
1894 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1895 map<string, string> variables;
1896 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1897 string otuCountsFileName = getOutputFileName("counts",variables);
1898 ofstream otuCountsFile;
1899 m->openOutputFile(otuCountsFileName, otuCountsFile);
1901 string bases = flowOrder;
1903 for(int i=0;i<numOTUs;i++){
1905 if (m->control_pressed) {
1908 //output the translated version of the centroid sequence for the otu
1909 if(otuCounts[i] > 0){
1910 int index = centroids[i];
1912 otuCountsFile << "ideal\t";
1913 for(int j=8;j<numFlowCells;j++){
1914 char base = bases[j % bases.length()];
1915 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
1916 otuCountsFile << base;
1919 otuCountsFile << endl;
1921 for(int j=0;j<nSeqsPerOTU[i];j++){
1922 int sequence = aaI[i][j];
1923 otuCountsFile << seqNameVector[sequence] << '\t';
1927 for(int k=0;k<lengths[sequence];k++){
1928 char base = bases[k % bases.length()];
1929 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
1931 for(int s=0;s<freq;s++){
1933 //otuCountsFile << base;
1936 otuCountsFile << newSeq.substr(4) << endl;
1938 otuCountsFile << endl;
1941 otuCountsFile.close();
1942 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
1945 catch(exception& e) {
1946 m->errorOut(e, "ShhherCommand", "writeClusters");
1952 //**********************************************************************************************************************
1954 int ShhherCommand::execute(){
1956 if (abort == true) { if (calledHelp) { return 0; } return 2; }
1958 getSingleLookUp(); if (m->control_pressed) { return 0; }
1959 getJointLookUp(); if (m->control_pressed) { return 0; }
1961 int numFiles = flowFileVector.size();
1963 if (numFiles < processors) { processors = numFiles; }
1965 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1966 if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName); }
1967 else { createProcesses(flowFileVector); } //each processor processes one file
1969 driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName);
1972 if(compositeFASTAFileName != ""){
1973 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
1974 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
1977 m->mothurOutEndLine();
1978 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
1979 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
1980 m->mothurOutEndLine();
1984 catch(exception& e) {
1985 m->errorOut(e, "ShhherCommand", "execute");
1990 //********************************************************************************************************************
1991 //sorts biggest to smallest
1992 inline bool compareFileSizes(string left, string right){
1997 //get num bytes in file
1998 string filename = left;
1999 pFile = fopen (filename.c_str(),"rb");
2000 string error = "Error opening " + filename;
2001 if (pFile==NULL) perror (error.c_str());
2003 fseek (pFile, 0, SEEK_END);
2004 leftsize=ftell (pFile);
2011 //get num bytes in file
2013 pFile2 = fopen (filename.c_str(),"rb");
2014 error = "Error opening " + filename;
2015 if (pFile2==NULL) perror (error.c_str());
2017 fseek (pFile2, 0, SEEK_END);
2018 rightsize=ftell (pFile2);
2022 return (leftsize > rightsize);
2024 /**************************************************************************************************/
2026 int ShhherCommand::createProcesses(vector<string> filenames){
2028 vector<int> processIDS;
2033 if (filenames.size() < processors) { processors = filenames.size(); }
2035 //sort file names by size to divide load better
2036 sort(filenames.begin(), filenames.end(), compareFileSizes);
2038 vector < vector <string> > dividedFiles; //dividedFiles[1] = vector of filenames for process 1...
2039 dividedFiles.resize(processors);
2041 //for each file, figure out which process will complete it
2042 //want to divide the load intelligently so the big files are spread between processes
2043 for (int i = 0; i < filenames.size(); i++) {
2044 int processToAssign = (i+1) % processors;
2045 if (processToAssign == 0) { processToAssign = processors; }
2047 dividedFiles[(processToAssign-1)].push_back(filenames[i]);
2050 //now lets reverse the order of ever other process, so we balance big files running with little ones
2051 for (int i = 0; i < processors; i++) {
2052 int remainder = ((i+1) % processors);
2053 if (remainder) { reverse(dividedFiles[i].begin(), dividedFiles[i].end()); }
2057 //divide the groups between the processors
2058 /*vector<linePair> lines;
2059 vector<int> numFilesToComplete;
2060 int numFilesPerProcessor = filenames.size() / processors;
2061 for (int i = 0; i < processors; i++) {
2062 int startIndex = i * numFilesPerProcessor;
2063 int endIndex = (i+1) * numFilesPerProcessor;
2064 if(i == (processors - 1)){ endIndex = filenames.size(); }
2065 lines.push_back(linePair(startIndex, endIndex));
2066 numFilesToComplete.push_back((endIndex-startIndex));
2069 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
2071 //loop through and create all the processes you want
2072 while (process != processors) {
2076 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
2078 }else if (pid == 0){
2079 num = driver(dividedFiles[process], compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp");
2081 //pass numSeqs to parent
2083 string tempFile = compositeFASTAFileName + toString(getpid()) + ".num.temp";
2084 m->openOutputFile(tempFile, out);
2090 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
2091 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
2097 driver(dividedFiles[0], compositeFASTAFileName, compositeNamesFileName);
2099 //force parent to wait until all the processes are done
2100 for (int i=0;i<processIDS.size();i++) {
2101 int temp = processIDS[i];
2107 //////////////////////////////////////////////////////////////////////////////////////////////////////
2109 /////////////////////// NOT WORKING, ACCESS VIOLATION ON READ OF FLOWGRAMS IN THREAD /////////////////
2111 //////////////////////////////////////////////////////////////////////////////////////////////////////
2112 //Windows version shared memory, so be careful when passing variables through the shhhFlowsData struct.
2113 //Above fork() will clone, so memory is separate, but that's not the case with windows,
2114 //////////////////////////////////////////////////////////////////////////////////////////////////////
2116 vector<shhhFlowsData*> pDataArray;
2117 DWORD dwThreadIdArray[processors-1];
2118 HANDLE hThreadArray[processors-1];
2120 //Create processor worker threads.
2121 for( int i=0; i<processors-1; i++ ){
2122 // Allocate memory for thread data.
2123 string extension = "";
2124 if (i != 0) { extension = toString(i) + ".temp"; }
2126 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);
2127 pDataArray.push_back(tempFlow);
2128 processIDS.push_back(i);
2130 hThreadArray[i] = CreateThread(NULL, 0, ShhhFlowsThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
2133 //using the main process as a worker saves time and memory
2135 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[processors-1].start, lines[processors-1].end);
2137 //Wait until all threads have terminated.
2138 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
2140 //Close all thread handles and free memory allocations.
2141 for(int i=0; i < pDataArray.size(); i++){
2142 for(int j=0; j < pDataArray[i]->outputNames.size(); j++){ outputNames.push_back(pDataArray[i]->outputNames[j]); }
2143 CloseHandle(hThreadArray[i]);
2144 delete pDataArray[i];
2149 for (int i=0;i<processIDS.size();i++) {
2151 string tempFile = compositeFASTAFileName + toString(processIDS[i]) + ".num.temp";
2152 m->openInputFile(tempFile, in);
2156 if (tempNum != dividedFiles[i+1].size()) {
2157 m->mothurOut("[ERROR]: main process expected " + toString(processIDS[i]) + " to complete " + toString(dividedFiles[i+1].size()) + " 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");
2160 in.close(); m->mothurRemove(tempFile);
2162 if (compositeFASTAFileName != "") {
2163 m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName);
2164 m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName);
2165 m->mothurRemove((compositeFASTAFileName + toString(processIDS[i]) + ".temp"));
2166 m->mothurRemove((compositeNamesFileName + toString(processIDS[i]) + ".temp"));
2173 catch(exception& e) {
2174 m->errorOut(e, "ShhherCommand", "createProcesses");
2178 /**************************************************************************************************/
2180 vector<string> ShhherCommand::parseFlowFiles(string filename){
2182 vector<string> files;
2186 m->openInputFile(filename, in);
2188 int thisNumFLows = 0;
2189 in >> thisNumFLows; m->gobble(in);
2192 if (m->control_pressed) { break; }
2195 string outputFileName = filename + toString(count) + ".temp";
2196 m->openOutputFile(outputFileName, out);
2197 out << thisNumFLows << endl;
2198 files.push_back(outputFileName);
2200 int numLinesWrote = 0;
2201 for (int i = 0; i < largeSize; i++) {
2202 if (in.eof()) { break; }
2203 string line = m->getline(in); m->gobble(in);
2204 out << line << endl;
2209 if (numLinesWrote == 0) { m->mothurRemove(outputFileName); files.pop_back(); }
2214 if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } files.clear(); }
2216 m->mothurOut("\nDivided " + filename + " into " + toString(files.size()) + " files.\n\n");
2220 catch(exception& e) {
2221 m->errorOut(e, "ShhherCommand", "parseFlowFiles");
2225 /**************************************************************************************************/
2227 int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName){
2230 int numCompleted = 0;
2232 for(int i=0;i<filenames.size();i++){
2234 if (m->control_pressed) { break; }
2236 vector<string> theseFlowFileNames; theseFlowFileNames.push_back(filenames[i]);
2237 if (large) { theseFlowFileNames = parseFlowFiles(filenames[i]); }
2239 if (m->control_pressed) { break; }
2241 double begClock = clock();
2242 unsigned long long begTime;
2244 string fileNameForOutput = filenames[i];
2246 for (int g = 0; g < theseFlowFileNames.size(); g++) {
2248 string flowFileName = theseFlowFileNames[g];
2249 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n");
2250 m->mothurOut("Reading flowgrams...\n");
2252 vector<string> seqNameVector;
2253 vector<int> lengths;
2254 vector<short> flowDataIntI;
2255 vector<double> flowDataPrI;
2256 map<string, int> nameMap;
2257 vector<short> uniqueFlowgrams;
2258 vector<int> uniqueCount;
2259 vector<int> mapSeqToUnique;
2260 vector<int> mapUniqueToSeq;
2261 vector<int> uniqueLengths;
2264 if (m->debug) { m->mothurOut("[DEBUG]: About to read flowgrams.\n"); }
2265 int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells);
2267 if (m->control_pressed) { break; }
2269 m->mothurOut("Identifying unique flowgrams...\n");
2270 int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI);
2272 if (m->control_pressed) { break; }
2274 m->mothurOut("Calculating distances between flowgrams...\n");
2275 string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
2276 begTime = time(NULL);
2279 flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2281 m->mothurOutEndLine();
2282 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
2285 string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
2286 createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq);
2288 if (m->control_pressed) { break; }
2290 m->mothurOut("\nClustering flowgrams...\n");
2291 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
2292 cluster(listFileName, distFileName, namesFileName);
2294 if (m->control_pressed) { break; }
2296 vector<int> otuData;
2297 vector<int> cumNumSeqs;
2298 vector<int> nSeqsPerOTU;
2299 vector<vector<int> > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2300 vector<vector<int> > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2301 vector<int> seqNumber; //tMaster->anP: the sequence id number sorted by OTU
2302 vector<int> seqIndex; //tMaster->anI; the index that corresponds to seqNumber
2305 int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap);
2307 if (m->control_pressed) { break; }
2309 m->mothurRemove(distFileName);
2310 m->mothurRemove(namesFileName);
2311 m->mothurRemove(listFileName);
2313 vector<double> dist; //adDist - distance of sequences to centroids
2314 vector<short> change; //did the centroid sequence change? 0 = no; 1 = yes
2315 vector<int> centroids; //the representative flowgram for each cluster m
2316 vector<double> weight;
2317 vector<double> singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2318 vector<int> nSeqsBreaks;
2319 vector<int> nOTUsBreaks;
2321 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
2323 dist.assign(numSeqs * numOTUs, 0);
2324 change.assign(numOTUs, 1);
2325 centroids.assign(numOTUs, -1);
2326 weight.assign(numOTUs, 0);
2327 singleTau.assign(numSeqs, 1.0);
2329 nSeqsBreaks.assign(2, 0);
2330 nOTUsBreaks.assign(2, 0);
2333 nSeqsBreaks[1] = numSeqs;
2334 nOTUsBreaks[1] = numOTUs;
2336 if (m->debug) { m->mothurOut("[DEBUG]: done allocating memory, about to denoise.\n"); }
2338 if (m->control_pressed) { break; }
2340 double maxDelta = 0;
2344 begTime = time(NULL);
2346 m->mothurOut("\nDenoising flowgrams...\n");
2347 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
2349 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
2351 if (m->control_pressed) { break; }
2353 double cycClock = clock();
2354 unsigned long long cycTime = time(NULL);
2355 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2357 if (m->control_pressed) { break; }
2359 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2361 if (m->control_pressed) { break; }
2363 maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight);
2365 if (m->control_pressed) { break; }
2367 double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight);
2369 if (m->control_pressed) { break; }
2371 checkCentroids(numOTUs, centroids, weight);
2373 if (m->control_pressed) { break; }
2375 calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths);
2377 if (m->control_pressed) { break; }
2381 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
2385 if (m->control_pressed) { break; }
2387 m->mothurOut("\nFinalizing...\n");
2388 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2390 if (m->debug) { m->mothurOut("[DEBUG]: done fill().\n"); }
2392 if (m->control_pressed) { break; }
2394 setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI);
2396 if (m->debug) { m->mothurOut("[DEBUG]: done setOTUs().\n"); }
2398 if (m->control_pressed) { break; }
2400 vector<int> otuCounts(numOTUs, 0);
2401 for(int j=0;j<numSeqs;j++) { otuCounts[otuData[j]]++; }
2403 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2405 if (m->debug) { m->mothurOut("[DEBUG]: done calcCentroidsDriver().\n"); }
2407 if (m->control_pressed) { break; }
2409 if ((large) && (g == 0)) { flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; }
2410 string thisOutputDir = outputDir;
2411 if (outputDir == "") { thisOutputDir = m->hasPath(flowFileName); }
2412 map<string, string> variables;
2413 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
2414 string qualityFileName = getOutputFileName("qfile",variables);
2415 string fastaFileName = getOutputFileName("fasta",variables);
2416 string nameFileName = getOutputFileName("name",variables);
2417 string otuCountsFileName = getOutputFileName("counts",variables);
2418 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
2419 int pos = fileRoot.find_first_of('.');
2420 string fileGroup = fileRoot;
2421 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
2422 string groupFileName = getOutputFileName("group",variables);
2425 writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; }
2426 writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; }
2427 writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; }
2428 writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; }
2429 writeGroups(groupFileName, fileGroup, numSeqs, seqNameVector); if (m->control_pressed) { break; }
2433 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0]));
2434 m->appendFiles(qualityFileName, getOutputFileName("qfile",variables));
2435 m->mothurRemove(qualityFileName);
2436 m->appendFiles(fastaFileName, getOutputFileName("fasta",variables));
2437 m->mothurRemove(fastaFileName);
2438 m->appendFiles(nameFileName, getOutputFileName("name",variables));
2439 m->mothurRemove(nameFileName);
2440 m->appendFiles(otuCountsFileName, getOutputFileName("counts",variables));
2441 m->mothurRemove(otuCountsFileName);
2442 m->appendFiles(groupFileName, getOutputFileName("group",variables));
2443 m->mothurRemove(groupFileName);
2445 m->mothurRemove(theseFlowFileNames[g]);
2450 m->mothurOut("Total time to process " + fileNameForOutput + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
2453 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
2455 return numCompleted;
2457 }catch(exception& e) {
2458 m->errorOut(e, "ShhherCommand", "driver");
2463 /**************************************************************************************************/
2464 int ShhherCommand::getFlowData(string filename, vector<string>& thisSeqNameVector, vector<int>& thisLengths, vector<short>& thisFlowDataIntI, map<string, int>& thisNameMap, int& numFlowCells){
2469 m->openInputFile(filename, flowFile);
2472 int currentNumFlowCells;
2474 thisSeqNameVector.clear();
2475 thisLengths.clear();
2476 thisFlowDataIntI.clear();
2477 thisNameMap.clear();
2480 flowFile >> numFlowTest;
2482 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); }
2483 else { convert(numFlowTest, numFlowCells); }
2485 if (m->debug) { m->mothurOut("[DEBUG]: numFlowCells = " + toString(numFlowCells) + ".\n"); }
2486 int index = 0;//pcluster
2487 while(!flowFile.eof()){
2489 if (m->control_pressed) { break; }
2491 flowFile >> seqName >> currentNumFlowCells;
2493 thisLengths.push_back(currentNumFlowCells);
2495 thisSeqNameVector.push_back(seqName);
2496 thisNameMap[seqName] = index++;//pcluster
2498 if (m->debug) { m->mothurOut("[DEBUG]: seqName = " + seqName + " length = " + toString(currentNumFlowCells) + " index = " + toString(index) + "\n"); }
2500 for(int i=0;i<numFlowCells;i++){
2501 flowFile >> intensity;
2502 if(intensity > 9.99) { intensity = 9.99; }
2503 int intI = int(100 * intensity + 0.0001);
2504 thisFlowDataIntI.push_back(intI);
2506 m->gobble(flowFile);
2510 int numSeqs = thisSeqNameVector.size();
2512 for(int i=0;i<numSeqs;i++){
2514 if (m->control_pressed) { break; }
2516 int iNumFlowCells = i * numFlowCells;
2517 for(int j=thisLengths[i];j<numFlowCells;j++){
2518 thisFlowDataIntI[iNumFlowCells + j] = 0;
2525 catch(exception& e) {
2526 m->errorOut(e, "ShhherCommand", "getFlowData");
2530 /**************************************************************************************************/
2532 int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int stopSeq, vector<int>& mapUniqueToSeq, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2535 ostringstream outStream;
2536 outStream.setf(ios::fixed, ios::floatfield);
2537 outStream.setf(ios::dec, ios::basefield);
2538 outStream.setf(ios::showpoint);
2539 outStream.precision(6);
2541 int begTime = time(NULL);
2542 double begClock = clock();
2544 for(int i=0;i<stopSeq;i++){
2546 if (m->control_pressed) { break; }
2548 for(int j=0;j<i;j++){
2549 float flowDistance = calcPairwiseDist(numFlowCells, mapUniqueToSeq[i], mapUniqueToSeq[j], mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2551 if(flowDistance < 1e-6){
2552 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
2554 else if(flowDistance <= cutoff){
2555 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
2559 m->mothurOutJustToScreen(toString(i) + "\t" + toString(time(NULL) - begTime));
2560 m->mothurOutJustToScreen("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)+"\n");
2564 ofstream distFile(distFileName.c_str());
2565 distFile << outStream.str();
2568 if (m->control_pressed) {}
2570 m->mothurOutJustToScreen(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
2571 m->mothurOutJustToScreen("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)+"\n");
2576 catch(exception& e) {
2577 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
2581 /**************************************************************************************************/
2583 float ShhherCommand::calcPairwiseDist(int numFlowCells, int seqA, int seqB, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2585 int minLength = lengths[mapSeqToUnique[seqA]];
2586 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
2588 int ANumFlowCells = seqA * numFlowCells;
2589 int BNumFlowCells = seqB * numFlowCells;
2593 for(int i=0;i<minLength;i++){
2595 if (m->control_pressed) { break; }
2597 int flowAIntI = flowDataIntI[ANumFlowCells + i];
2598 float flowAPrI = flowDataPrI[ANumFlowCells + i];
2600 int flowBIntI = flowDataIntI[BNumFlowCells + i];
2601 float flowBPrI = flowDataPrI[BNumFlowCells + i];
2602 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
2605 dist /= (float) minLength;
2608 catch(exception& e) {
2609 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
2614 /**************************************************************************************************/
2616 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){
2619 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
2620 uniqueCount.assign(numSeqs, 0); // anWeights
2621 uniqueLengths.assign(numSeqs, 0);
2622 mapSeqToUnique.assign(numSeqs, -1);
2623 mapUniqueToSeq.assign(numSeqs, -1);
2625 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
2627 for(int i=0;i<numSeqs;i++){
2629 if (m->control_pressed) { break; }
2633 vector<short> current(numFlowCells);
2634 for(int j=0;j<numFlowCells;j++){
2635 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
2638 for(int j=0;j<numUniques;j++){
2639 int offset = j * numFlowCells;
2643 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
2644 else { shorterLength = uniqueLengths[j]; }
2646 for(int k=0;k<shorterLength;k++){
2647 if(current[k] != uniqueFlowgrams[offset + k]){
2654 mapSeqToUnique[i] = j;
2657 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
2663 if(index == numUniques){
2664 uniqueLengths[numUniques] = lengths[i];
2665 uniqueCount[numUniques] = 1;
2666 mapSeqToUnique[i] = numUniques;//anMap
2667 mapUniqueToSeq[numUniques] = i;//anF
2669 for(int k=0;k<numFlowCells;k++){
2670 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
2671 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
2677 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
2678 uniqueLengths.resize(numUniques);
2680 flowDataPrI.resize(numSeqs * numFlowCells, 0);
2681 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
2685 catch(exception& e) {
2686 m->errorOut(e, "ShhherCommand", "getUniques");
2690 /**************************************************************************************************/
2691 int ShhherCommand::createNamesFile(int numSeqs, int numUniques, string filename, vector<string>& seqNameVector, vector<int>& mapSeqToUnique, vector<int>& mapUniqueToSeq){
2694 vector<string> duplicateNames(numUniques, "");
2695 for(int i=0;i<numSeqs;i++){
2696 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
2700 m->openOutputFile(filename, nameFile);
2702 for(int i=0;i<numUniques;i++){
2704 if (m->control_pressed) { break; }
2706 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2707 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2714 catch(exception& e) {
2715 m->errorOut(e, "ShhherCommand", "createNamesFile");
2719 //**********************************************************************************************************************
2721 int ShhherCommand::cluster(string filename, string distFileName, string namesFileName){
2724 ReadMatrix* read = new ReadColumnMatrix(distFileName);
2725 read->setCutoff(cutoff);
2727 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
2728 clusterNameMap->readMap();
2729 read->read(clusterNameMap);
2731 ListVector* list = read->getListVector();
2732 SparseDistanceMatrix* matrix = read->getDMatrix();
2735 delete clusterNameMap;
2737 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
2739 float adjust = -1.0;
2740 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest", adjust);
2741 string tag = cluster->getTag();
2743 double clusterCutoff = cutoff;
2744 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
2746 if (m->control_pressed) { break; }
2748 cluster->update(clusterCutoff);
2751 list->setLabel(toString(cutoff));
2754 m->openOutputFile(filename, listFile);
2755 list->print(listFile);
2758 delete matrix; delete cluster; delete rabund; delete list;
2762 catch(exception& e) {
2763 m->errorOut(e, "ShhherCommand", "cluster");
2767 /**************************************************************************************************/
2769 int ShhherCommand::getOTUData(int numSeqs, string fileName, vector<int>& otuData,
2770 vector<int>& cumNumSeqs,
2771 vector<int>& nSeqsPerOTU,
2772 vector<vector<int> >& aaP, //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2773 vector<vector<int> >& aaI, //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2774 vector<int>& seqNumber, //tMaster->anP: the sequence id number sorted by OTU
2775 vector<int>& seqIndex,
2776 map<string, int>& nameMap){
2780 m->openInputFile(fileName, listFile);
2784 listFile >> label >> numOTUs;
2786 if (m->debug) { m->mothurOut("[DEBUG]: Getting OTU Data...\n"); }
2788 otuData.assign(numSeqs, 0);
2789 cumNumSeqs.assign(numOTUs, 0);
2790 nSeqsPerOTU.assign(numOTUs, 0);
2791 aaP.clear();aaP.resize(numOTUs);
2797 string singleOTU = "";
2799 for(int i=0;i<numOTUs;i++){
2801 if (m->control_pressed) { break; }
2802 if (m->debug) { m->mothurOut("[DEBUG]: processing OTU " + toString(i) + ".\n"); }
2804 listFile >> singleOTU;
2806 istringstream otuString(singleOTU);
2810 string seqName = "";
2812 for(int j=0;j<singleOTU.length();j++){
2813 char letter = otuString.get();
2819 map<string,int>::iterator nmIt = nameMap.find(seqName);
2820 int index = nmIt->second;
2822 nameMap.erase(nmIt);
2826 aaP[i].push_back(index);
2831 map<string,int>::iterator nmIt = nameMap.find(seqName);
2833 int index = nmIt->second;
2834 nameMap.erase(nmIt);
2838 aaP[i].push_back(index);
2843 sort(aaP[i].begin(), aaP[i].end());
2844 for(int j=0;j<nSeqsPerOTU[i];j++){
2845 seqNumber.push_back(aaP[i][j]);
2847 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
2848 aaP[i].push_back(0);
2854 for(int i=1;i<numOTUs;i++){
2855 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
2858 seqIndex = seqNumber;
2865 catch(exception& e) {
2866 m->errorOut(e, "ShhherCommand", "getOTUData");
2870 /**************************************************************************************************/
2872 int ShhherCommand::calcCentroidsDriver(int numOTUs,
2873 vector<int>& cumNumSeqs,
2874 vector<int>& nSeqsPerOTU,
2875 vector<int>& seqIndex,
2876 vector<short>& change, //did the centroid sequence change? 0 = no; 1 = yes
2877 vector<int>& centroids, //the representative flowgram for each cluster m
2878 vector<double>& singleTau, //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2879 vector<int>& mapSeqToUnique,
2880 vector<short>& uniqueFlowgrams,
2881 vector<short>& flowDataIntI,
2882 vector<int>& lengths,
2884 vector<int>& seqNumber){
2886 //this function gets the most likely homopolymer length at a flow position for a group of sequences
2891 for(int i=0;i<numOTUs;i++){
2893 if (m->control_pressed) { break; }
2897 int minFlowGram = 100000000;
2898 double minFlowValue = 1e8;
2899 change[i] = 0; //FALSE
2901 for(int j=0;j<nSeqsPerOTU[i];j++){
2902 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
2905 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
2906 vector<double> adF(nSeqsPerOTU[i]);
2907 vector<int> anL(nSeqsPerOTU[i]);
2909 for(int j=0;j<nSeqsPerOTU[i];j++){
2910 int index = cumNumSeqs[i] + j;
2911 int nI = seqIndex[index];
2912 int nIU = mapSeqToUnique[nI];
2915 for(k=0;k<position;k++){
2921 anL[position] = nIU;
2922 adF[position] = 0.0000;
2927 for(int j=0;j<nSeqsPerOTU[i];j++){
2928 int index = cumNumSeqs[i] + j;
2929 int nI = seqIndex[index];
2931 double tauValue = singleTau[seqNumber[index]];
2933 for(int k=0;k<position;k++){
2934 double dist = getDistToCentroid(anL[k], nI, lengths[nI], uniqueFlowgrams, flowDataIntI, numFlowCells);
2935 adF[k] += dist * tauValue;
2939 for(int j=0;j<position;j++){
2940 if(adF[j] < minFlowValue){
2942 minFlowValue = adF[j];
2946 if(centroids[i] != anL[minFlowGram]){
2948 centroids[i] = anL[minFlowGram];
2951 else if(centroids[i] != -1){
2959 catch(exception& e) {
2960 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
2964 /**************************************************************************************************/
2966 double ShhherCommand::getDistToCentroid(int cent, int flow, int length, vector<short>& uniqueFlowgrams,
2967 vector<short>& flowDataIntI, int numFlowCells){
2970 int flowAValue = cent * numFlowCells;
2971 int flowBValue = flow * numFlowCells;
2975 for(int i=0;i<length;i++){
2976 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
2981 return dist / (double)length;
2983 catch(exception& e) {
2984 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
2988 /**************************************************************************************************/
2990 double ShhherCommand::getNewWeights(int numOTUs, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU, vector<double>& singleTau, vector<int>& seqNumber, vector<double>& weight){
2993 double maxChange = 0;
2995 for(int i=0;i<numOTUs;i++){
2997 if (m->control_pressed) { break; }
2999 double difference = weight[i];
3002 for(int j=0;j<nSeqsPerOTU[i];j++){
3003 int index = cumNumSeqs[i] + j;
3004 double tauValue = singleTau[seqNumber[index]];
3005 weight[i] += tauValue;
3008 difference = fabs(weight[i] - difference);
3009 if(difference > maxChange){ maxChange = difference; }
3013 catch(exception& e) {
3014 m->errorOut(e, "ShhherCommand", "getNewWeights");
3019 /**************************************************************************************************/
3021 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){
3025 vector<long double> P(numSeqs, 0);
3028 for(int i=0;i<numOTUs;i++){
3029 if(weight[i] > MIN_WEIGHT){
3035 for(int i=0;i<numOTUs;i++){
3037 if (m->control_pressed) { break; }
3039 for(int j=0;j<nSeqsPerOTU[i];j++){
3040 int index = cumNumSeqs[i] + j;
3041 int nI = seqIndex[index];
3042 double singleDist = dist[seqNumber[index]];
3044 P[nI] += weight[i] * exp(-singleDist * sigma);
3048 for(int i=0;i<numSeqs;i++){
3049 if(P[i] == 0){ P[i] = DBL_EPSILON; }
3054 nLL = nLL -(double)numSeqs * log(sigma);
3058 catch(exception& e) {
3059 m->errorOut(e, "ShhherCommand", "getNewWeights");
3064 /**************************************************************************************************/
3066 int ShhherCommand::checkCentroids(int numOTUs, vector<int>& centroids, vector<double>& weight){
3068 vector<int> unique(numOTUs, 1);
3070 for(int i=0;i<numOTUs;i++){
3071 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
3076 for(int i=0;i<numOTUs;i++){
3078 if (m->control_pressed) { break; }
3081 for(int j=i+1;j<numOTUs;j++){
3084 if(centroids[j] == centroids[i]){
3088 weight[i] += weight[j];
3098 catch(exception& e) {
3099 m->errorOut(e, "ShhherCommand", "checkCentroids");
3103 /**************************************************************************************************/
3105 void ShhherCommand::calcNewDistances(int numSeqs, int numOTUs, vector<int>& nSeqsPerOTU, vector<double>& dist,
3106 vector<double>& weight, vector<short>& change, vector<int>& centroids,
3107 vector<vector<int> >& aaP, vector<double>& singleTau, vector<vector<int> >& aaI,
3108 vector<int>& seqNumber, vector<int>& seqIndex,
3109 vector<short>& uniqueFlowgrams,
3110 vector<short>& flowDataIntI, int numFlowCells, vector<int>& lengths){
3115 vector<double> newTau(numOTUs,0);
3116 vector<double> norms(numSeqs, 0);
3117 nSeqsPerOTU.assign(numOTUs, 0);
3119 for(int i=0;i<numSeqs;i++){
3121 if (m->control_pressed) { break; }
3123 int indexOffset = i * numOTUs;
3125 double offset = 1e8;
3127 for(int j=0;j<numOTUs;j++){
3129 if(weight[j] > MIN_WEIGHT && change[j] == 1){
3130 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i], uniqueFlowgrams, flowDataIntI, numFlowCells);
3133 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
3134 offset = dist[indexOffset + j];
3138 for(int j=0;j<numOTUs;j++){
3139 if(weight[j] > MIN_WEIGHT){
3140 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
3141 norms[i] += newTau[j];
3148 for(int j=0;j<numOTUs;j++){
3149 newTau[j] /= norms[i];
3152 for(int j=0;j<numOTUs;j++){
3153 if(newTau[j] > MIN_TAU){
3155 int oldTotal = total;
3159 singleTau.resize(total, 0);
3160 seqNumber.resize(total, 0);
3161 seqIndex.resize(total, 0);
3163 singleTau[oldTotal] = newTau[j];
3165 aaP[j][nSeqsPerOTU[j]] = oldTotal;
3166 aaI[j][nSeqsPerOTU[j]] = i;
3174 catch(exception& e) {
3175 m->errorOut(e, "ShhherCommand", "calcNewDistances");
3179 /**************************************************************************************************/
3181 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){
3184 for(int i=0;i<numOTUs;i++){
3186 if (m->control_pressed) { return 0; }
3188 cumNumSeqs[i] = index;
3189 for(int j=0;j<nSeqsPerOTU[i];j++){
3190 seqNumber[index] = aaP[i][j];
3191 seqIndex[index] = aaI[i][j];
3199 catch(exception& e) {
3200 m->errorOut(e, "ShhherCommand", "fill");
3204 /**************************************************************************************************/
3206 void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector<int>& seqNumber, vector<int>& seqIndex, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU,
3207 vector<int>& otuData, vector<double>& singleTau, vector<double>& dist, vector<vector<int> >& aaP, vector<vector<int> >& aaI){
3210 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
3212 for(int i=0;i<numOTUs;i++){
3214 if (m->control_pressed) { break; }
3216 for(int j=0;j<nSeqsPerOTU[i];j++){
3217 int index = cumNumSeqs[i] + j;
3218 double tauValue = singleTau[seqNumber[index]];
3219 int sIndex = seqIndex[index];
3220 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
3224 for(int i=0;i<numSeqs;i++){
3225 double maxTau = -1.0000;
3227 for(int j=0;j<numOTUs;j++){
3228 if(bigTauMatrix[i * numOTUs + j] > maxTau){
3229 maxTau = bigTauMatrix[i * numOTUs + j];
3234 otuData[i] = maxOTU;
3237 nSeqsPerOTU.assign(numOTUs, 0);
3239 for(int i=0;i<numSeqs;i++){
3240 int index = otuData[i];
3242 singleTau[i] = 1.0000;
3245 aaP[index][nSeqsPerOTU[index]] = i;
3246 aaI[index][nSeqsPerOTU[index]] = i;
3248 nSeqsPerOTU[index]++;
3251 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
3253 catch(exception& e) {
3254 m->errorOut(e, "ShhherCommand", "setOTUs");
3258 /**************************************************************************************************/
3260 void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string qualityFileName, vector<int> otuCounts, vector<int>& nSeqsPerOTU, vector<int>& seqNumber,
3261 vector<double>& singleTau, vector<short>& flowDataIntI, vector<short>& uniqueFlowgrams, vector<int>& cumNumSeqs,
3262 vector<int>& mapUniqueToSeq, vector<string>& seqNameVector, vector<int>& centroids, vector<vector<int> >& aaI){
3266 ofstream qualityFile;
3267 m->openOutputFile(qualityFileName, qualityFile);
3269 qualityFile.setf(ios::fixed, ios::floatfield);
3270 qualityFile.setf(ios::showpoint);
3271 qualityFile << setprecision(6);
3273 vector<vector<int> > qualities(numOTUs);
3274 vector<double> pr(HOMOPS, 0);
3277 for(int i=0;i<numOTUs;i++){
3279 if (m->control_pressed) { break; }
3283 if(nSeqsPerOTU[i] > 0){
3285 while(index < numFlowCells){
3287 double maxPrValue = 1e8;
3288 short maxPrIndex = -1;
3289 double count = 0.0000;
3291 pr.assign(HOMOPS, 0);
3293 for(int j=0;j<nSeqsPerOTU[i];j++){
3294 int lIndex = cumNumSeqs[i] + j;
3295 double tauValue = singleTau[seqNumber[lIndex]];
3296 int sequenceIndex = aaI[i][j];
3297 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
3301 for(int s=0;s<HOMOPS;s++){
3302 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
3306 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
3307 maxPrValue = pr[maxPrIndex];
3309 if(count > MIN_COUNT){
3311 double norm = 0.0000;
3313 for(int s=0;s<HOMOPS;s++){
3314 norm += exp(-(pr[s] - maxPrValue));
3317 for(int s=1;s<=maxPrIndex;s++){
3319 double temp = 0.0000;
3321 U += exp(-(pr[s-1]-maxPrValue))/norm;
3329 temp = floor(-10 * temp);
3330 value = (int)floor(temp);
3331 if(value > 100){ value = 100; }
3333 qualities[i].push_back((int)value);
3344 if(otuCounts[i] > 0){
3345 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
3346 //need to get past the first four bases
3347 for (int j = 4; j < qualities[i].size(); j++) { qualityFile << qualities[i][j] << ' '; }
3348 qualityFile << endl;
3351 qualityFile.close();
3352 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
3355 catch(exception& e) {
3356 m->errorOut(e, "ShhherCommand", "writeQualities");
3361 /**************************************************************************************************/
3363 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){
3367 m->openOutputFile(fastaFileName, fastaFile);
3369 vector<string> names(numOTUs, "");
3371 for(int i=0;i<numOTUs;i++){
3373 if (m->control_pressed) { break; }
3375 int index = centroids[i];
3377 if(otuCounts[i] > 0){
3378 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
3382 for(int j=0;j<numFlowCells;j++){
3384 char base = flowOrder[j % flowOrder.length()];
3385 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
3390 if (newSeq.length() >= 4) { fastaFile << newSeq.substr(4) << endl; }
3391 else { fastaFile << "NNNN" << endl; }
3396 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
3398 if(thisCompositeFASTAFileName != ""){
3399 m->appendFiles(fastaFileName, thisCompositeFASTAFileName);
3402 catch(exception& e) {
3403 m->errorOut(e, "ShhherCommand", "writeSequences");
3408 /**************************************************************************************************/
3410 void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string nameFileName, vector<int> otuCounts, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& nSeqsPerOTU){
3414 m->openOutputFile(nameFileName, nameFile);
3416 for(int i=0;i<numOTUs;i++){
3418 if (m->control_pressed) { break; }
3420 if(otuCounts[i] > 0){
3421 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
3423 for(int j=1;j<nSeqsPerOTU[i];j++){
3424 nameFile << ',' << seqNameVector[aaI[i][j]];
3431 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
3434 if(thisCompositeNamesFileName != ""){
3435 m->appendFiles(nameFileName, thisCompositeNamesFileName);
3438 catch(exception& e) {
3439 m->errorOut(e, "ShhherCommand", "writeNames");
3444 /**************************************************************************************************/
3446 void ShhherCommand::writeGroups(string groupFileName, string fileRoot, int numSeqs, vector<string>& seqNameVector){
3449 m->openOutputFile(groupFileName, groupFile);
3451 for(int i=0;i<numSeqs;i++){
3452 if (m->control_pressed) { break; }
3453 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
3456 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
3459 catch(exception& e) {
3460 m->errorOut(e, "ShhherCommand", "writeGroups");
3465 /**************************************************************************************************/
3467 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){
3469 ofstream otuCountsFile;
3470 m->openOutputFile(otuCountsFileName, otuCountsFile);
3472 string bases = flowOrder;
3474 for(int i=0;i<numOTUs;i++){
3476 if (m->control_pressed) {
3479 //output the translated version of the centroid sequence for the otu
3480 if(otuCounts[i] > 0){
3481 int index = centroids[i];
3483 otuCountsFile << "ideal\t";
3484 for(int j=8;j<numFlowCells;j++){
3485 char base = bases[j % bases.length()];
3486 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
3487 otuCountsFile << base;
3490 otuCountsFile << endl;
3492 for(int j=0;j<nSeqsPerOTU[i];j++){
3493 int sequence = aaI[i][j];
3494 otuCountsFile << seqNameVector[sequence] << '\t';
3498 for(int k=0;k<lengths[sequence];k++){
3499 char base = bases[k % bases.length()];
3500 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
3502 for(int s=0;s<freq;s++){
3504 //otuCountsFile << base;
3508 if (newSeq.length() >= 4) { otuCountsFile << newSeq.substr(4) << endl; }
3509 else { otuCountsFile << "NNNN" << endl; }
3511 otuCountsFile << endl;
3514 otuCountsFile.close();
3515 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
3518 catch(exception& e) {
3519 m->errorOut(e, "ShhherCommand", "writeClusters");
3524 /**************************************************************************************************/
3526 void ShhherCommand::getSingleLookUp(){
3528 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
3529 singleLookUp.assign(HOMOPS * NUMBINS, 0);
3532 ifstream lookUpFile;
3533 m->openInputFile(lookupFileName, lookUpFile);
3535 for(int i=0;i<HOMOPS;i++){
3537 if (m->control_pressed) { break; }
3540 lookUpFile >> logFracFreq;
3542 for(int j=0;j<NUMBINS;j++) {
3543 lookUpFile >> singleLookUp[index];
3549 catch(exception& e) {
3550 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
3555 /**************************************************************************************************/
3557 void ShhherCommand::getJointLookUp(){
3560 // the most likely joint probability (-log) that two intenities have the same polymer length
3561 jointLookUp.resize(NUMBINS * NUMBINS, 0);
3563 for(int i=0;i<NUMBINS;i++){
3565 if (m->control_pressed) { break; }
3567 for(int j=0;j<NUMBINS;j++){
3569 double minSum = 100000000;
3571 for(int k=0;k<HOMOPS;k++){
3572 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
3574 if(sum < minSum) { minSum = sum; }
3576 jointLookUp[i * NUMBINS + j] = minSum;
3580 catch(exception& e) {
3581 m->errorOut(e, "ShhherCommand", "getJointLookUp");
3586 /**************************************************************************************************/
3588 double ShhherCommand::getProbIntensity(int intIntensity){
3591 double minNegLogProb = 100000000;
3594 for(int i=0;i<HOMOPS;i++){//loop signal strength
3596 if (m->control_pressed) { break; }
3598 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
3599 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
3602 return minNegLogProb;
3604 catch(exception& e) {
3605 m->errorOut(e, "ShhherCommand", "getProbIntensity");