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 lookupFileName = "LookUp_Titanium.pat";
273 ableToOpen = m->openInputFile(lookupFileName, in, "noerror");
276 //if you can't open it, try input location
277 if (ableToOpen == 1) {
278 if (inputDir != "") { //default path is set
279 string tryPath = inputDir + lookupFileName;
280 m->mothurOut("Unable to open " + lookupFileName + ". Trying input directory " + tryPath); m->mothurOutEndLine();
282 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
284 lookupFileName = tryPath;
288 //if you can't open it, try default location
289 if (ableToOpen == 1) {
290 if (m->getDefaultPath() != "") { //default path is set
291 string tryPath = m->getDefaultPath() + m->getSimpleName(lookupFileName);
292 m->mothurOut("Unable to open " + lookupFileName + ". Trying default " + tryPath); m->mothurOutEndLine();
294 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
296 lookupFileName = tryPath;
300 //if you can't open it its not in current working directory or inputDir, try mothur excutable location
301 if (ableToOpen == 1) {
302 string exepath = m->argv;
303 string tempPath = exepath;
304 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
305 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
307 string tryPath = m->getFullPathName(exepath) + m->getSimpleName(lookupFileName);
308 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
310 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
312 lookupFileName = tryPath;
315 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
317 else if(temp == "not open") {
319 lookupFileName = validParameter.validFile(parameters, "lookup", false);
321 //if you can't open it its not inputDir, try mothur excutable location
322 string exepath = m->argv;
323 string tempPath = exepath;
324 for (int i = 0; i < exepath.length(); i++) { tempPath[i] = tolower(exepath[i]); }
325 exepath = exepath.substr(0, (tempPath.find_last_of('m')));
327 string tryPath = m->getFullPathName(exepath) + lookupFileName;
328 m->mothurOut("Unable to open " + lookupFileName + ". Trying mothur's executable location " + tryPath); m->mothurOutEndLine();
330 int ableToOpen = m->openInputFile(tryPath, in2, "noerror");
332 lookupFileName = tryPath;
334 if (ableToOpen == 1) { m->mothurOut("Unable to open " + lookupFileName + "."); m->mothurOutEndLine(); abort=true; }
335 }else { lookupFileName = temp; }
337 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
338 m->setProcessors(temp);
339 m->mothurConvert(temp, processors);
341 temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found"){ temp = "0.01"; }
342 m->mothurConvert(temp, cutoff);
344 temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){ temp = "0.000001"; }
345 m->mothurConvert(temp, minDelta);
347 temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; }
348 m->mothurConvert(temp, maxIters);
350 temp = validParameter.validFile(parameters, "large", false); if (temp == "not found"){ temp = "0"; }
351 m->mothurConvert(temp, largeSize);
352 if (largeSize != 0) { large = true; }
353 else { large = false; }
354 if (largeSize < 0) { m->mothurOut("The value of the large cannot be negative.\n"); }
357 if (large) { m->mothurOut("The large parameter is not available with the MPI-Enabled version.\n"); large=false; }
361 temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; }
362 m->mothurConvert(temp, sigma);
364 temp = validParameter.validFile(parameters, "order", false); if (temp == "not found"){ temp = "A"; }
365 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;
368 if (toupper(temp[0]) == 'A') { flowOrder = "TACG"; }
369 else if(toupper(temp[0]) == 'B'){
370 flowOrder = "TACGTACGTACGATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATAGATCGCATGACGATCGCATATCGTCAGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGTAGTCGAGCATCATCTGACGCAGTACGTGCATGATCTCAGTCAGCAGCTATGTCAGTGCATGCATAGATCGCATGACGATCGCATATCGTCAGTGCAGTGACTGATCGTCATCAGCTAGCATCGACTGCATGATCTCAGTCAGCAGC"; }
371 else if(toupper(temp[0]) == 'I'){
372 flowOrder = "TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGC"; }
374 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;
384 catch(exception& e) {
385 m->errorOut(e, "ShhherCommand", "ShhherCommand");
389 //**********************************************************************************************************************
391 int ShhherCommand::execute(){
393 if (abort == true) { if (calledHelp) { return 0; } return 2; }
400 for(int i=1;i<ncpus;i++){
401 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
403 if(abort == 1){ return 0; }
407 m->mothurOut("\nGetting preliminary data...\n");
408 getSingleLookUp(); if (m->control_pressed) { return 0; }
409 getJointLookUp(); if (m->control_pressed) { return 0; }
411 vector<string> flowFileVector;
412 if(flowFilesFileName != "not found"){
415 ifstream flowFilesFile;
416 m->openInputFile(flowFilesFileName, flowFilesFile);
417 while(flowFilesFile){
418 fName = m->getline(flowFilesFile);
419 flowFileVector.push_back(fName);
420 m->gobble(flowFilesFile);
424 flowFileVector.push_back(flowFileName);
427 int numFiles = flowFileVector.size();
429 for(int i=1;i<ncpus;i++){
430 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
433 for(int i=0;i<numFiles;i++){
435 if (m->control_pressed) { break; }
437 double begClock = clock();
438 unsigned long long begTime = time(NULL);
440 flowFileName = flowFileVector[i];
442 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
443 m->mothurOut("Reading flowgrams...\n");
446 if (m->control_pressed) { break; }
448 m->mothurOut("Identifying unique flowgrams...\n");
451 if (m->control_pressed) { break; }
453 m->mothurOut("Calculating distances between flowgrams...\n");
455 strcpy(fileName, flowFileName.c_str());
457 for(int i=1;i<ncpus;i++){
458 MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
460 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
461 MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
462 MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
463 MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
464 MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
465 MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
466 MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
467 MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
468 MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
469 MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
472 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
474 if (m->control_pressed) { break; }
477 for(int i=1;i<ncpus;i++){
478 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
480 m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
481 m->mothurRemove((distFileName + ".temp." + toString(i)));
484 string namesFileName = createNamesFile();
486 if (m->control_pressed) { break; }
488 m->mothurOut("\nClustering flowgrams...\n");
489 string listFileName = cluster(distFileName, namesFileName);
491 if (m->control_pressed) { break; }
495 getOTUData(listFileName);
497 m->mothurRemove(distFileName);
498 m->mothurRemove(namesFileName);
499 m->mothurRemove(listFileName);
501 if (m->control_pressed) { break; }
505 if (m->control_pressed) { break; }
508 for(int i=1;i<ncpus;i++){
509 MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
510 MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
511 MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
512 MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
515 if (m->control_pressed) { break; }
520 int numOTUsOnCPU = numOTUs / ncpus;
521 int numSeqsOnCPU = numSeqs / ncpus;
522 m->mothurOut("\nDenoising flowgrams...\n");
523 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
525 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
527 double cycClock = clock();
528 unsigned long long cycTime = time(NULL);
531 if (m->control_pressed) { break; }
533 int total = singleTau.size();
534 for(int i=1;i<ncpus;i++){
535 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
536 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
537 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
539 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
540 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
541 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
542 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
543 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
546 calcCentroidsDriver(0, numOTUsOnCPU);
548 for(int i=1;i<ncpus;i++){
549 int otuStart = i * numOTUs / ncpus;
550 int otuStop = (i + 1) * numOTUs / ncpus;
552 vector<int> tempCentroids(numOTUs, 0);
553 vector<short> tempChange(numOTUs, 0);
555 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
556 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
558 for(int j=otuStart;j<otuStop;j++){
559 centroids[j] = tempCentroids[j];
560 change[j] = tempChange[j];
564 maxDelta = getNewWeights(); if (m->control_pressed) { break; }
565 double nLL = getLikelihood(); if (m->control_pressed) { break; }
566 checkCentroids(); if (m->control_pressed) { break; }
568 for(int i=1;i<ncpus;i++){
569 MPI_Send(¢roids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
570 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
571 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
574 calcNewDistancesParent(0, numSeqsOnCPU);
576 total = singleTau.size();
578 for(int i=1;i<ncpus;i++){
580 int seqStart = i * numSeqs / ncpus;
581 int seqStop = (i + 1) * numSeqs / ncpus;
583 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
585 vector<int> childSeqIndex(childTotal, 0);
586 vector<double> childSingleTau(childTotal, 0);
587 vector<double> childDist(numSeqs * numOTUs, 0);
588 vector<int> otuIndex(childTotal, 0);
590 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
591 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
592 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
593 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
595 int oldTotal = total;
597 singleTau.resize(total, 0);
598 seqIndex.resize(total, 0);
599 seqNumber.resize(total, 0);
603 for(int j=oldTotal;j<total;j++){
604 int otuI = otuIndex[childIndex];
605 int seqI = childSeqIndex[childIndex];
607 singleTau[j] = childSingleTau[childIndex];
609 aaP[otuI][nSeqsPerOTU[otuI]] = j;
610 aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
615 int index = seqStart * numOTUs;
616 for(int j=seqStart;j<seqStop;j++){
617 for(int k=0;k<numOTUs;k++){
618 dist[index] = childDist[index];
626 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
628 if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
630 for(int i=1;i<ncpus;i++){
631 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
636 for(int i=1;i<ncpus;i++){
637 MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
643 if (m->control_pressed) { break; }
645 m->mothurOut("\nFinalizing...\n");
648 if (m->control_pressed) { break; }
652 vector<int> otuCounts(numOTUs, 0);
653 for(int i=0;i<numSeqs;i++) { otuCounts[otuData[i]]++; }
654 calcCentroidsDriver(0, numOTUs);
656 if (m->control_pressed) { break; }
658 writeQualities(otuCounts); if (m->control_pressed) { break; }
659 writeSequences(otuCounts); if (m->control_pressed) { break; }
660 writeNames(otuCounts); if (m->control_pressed) { break; }
661 writeClusters(otuCounts); if (m->control_pressed) { break; }
662 writeGroups(); if (m->control_pressed) { break; }
665 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
671 MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
672 if(abort){ return 0; }
675 MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
677 for(int i=0;i<numFiles;i++){
679 if (m->control_pressed) { break; }
681 //Now into the pyrodist part
685 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
686 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
687 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
688 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
690 flowDataIntI.resize(numSeqs * numFlowCells);
691 flowDataPrI.resize(numSeqs * numFlowCells);
692 mapUniqueToSeq.resize(numSeqs);
693 mapSeqToUnique.resize(numSeqs);
694 lengths.resize(numSeqs);
695 jointLookUp.resize(NUMBINS * NUMBINS);
697 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
698 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
699 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
700 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
701 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
702 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
703 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
705 flowFileName = string(fileName);
706 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
707 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
709 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
711 if (m->control_pressed) { break; }
714 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
716 //Now into the pyronoise part
717 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
719 singleLookUp.resize(HOMOPS * NUMBINS);
720 uniqueFlowgrams.resize(numUniques * numFlowCells);
721 weight.resize(numOTUs);
722 centroids.resize(numOTUs);
723 change.resize(numOTUs);
724 dist.assign(numOTUs * numSeqs, 0);
725 nSeqsPerOTU.resize(numOTUs);
726 cumNumSeqs.resize(numOTUs);
728 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
729 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
730 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
732 int startOTU = pid * numOTUs / ncpus;
733 int endOTU = (pid + 1) * numOTUs / ncpus;
735 int startSeq = pid * numSeqs / ncpus;
736 int endSeq = (pid + 1) * numSeqs /ncpus;
742 if (m->control_pressed) { break; }
744 MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
745 singleTau.assign(total, 0.0000);
746 seqNumber.assign(total, 0);
747 seqIndex.assign(total, 0);
749 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
750 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
751 MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
752 MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
753 MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
754 MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
755 MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
757 calcCentroidsDriver(startOTU, endOTU);
759 MPI_Send(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
760 MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
762 MPI_Recv(¢roids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
763 MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
764 MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
766 vector<int> otuIndex(total, 0);
767 calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
768 total = otuIndex.size();
770 MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
771 MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
772 MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
773 MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
774 MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
776 MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
781 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
783 MPI_Barrier(MPI_COMM_WORLD);
786 if(compositeFASTAFileName != ""){
787 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
788 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
791 m->mothurOutEndLine();
792 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
793 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
794 m->mothurOutEndLine();
799 catch(exception& e) {
800 m->errorOut(e, "ShhherCommand", "execute");
804 /**************************************************************************************************/
805 string ShhherCommand::createNamesFile(){
808 vector<string> duplicateNames(numUniques, "");
809 for(int i=0;i<numSeqs;i++){
810 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
812 map<string, string> variables;
813 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName));
814 string nameFileName = getOutputFileName("name",variables);
817 m->openOutputFile(nameFileName, nameFile);
819 for(int i=0;i<numUniques;i++){
821 if (m->control_pressed) { break; }
823 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
824 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
830 catch(exception& e) {
831 m->errorOut(e, "ShhherCommand", "createNamesFile");
835 /**************************************************************************************************/
837 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
839 ostringstream outStream;
840 outStream.setf(ios::fixed, ios::floatfield);
841 outStream.setf(ios::dec, ios::basefield);
842 outStream.setf(ios::showpoint);
843 outStream.precision(6);
845 int begTime = time(NULL);
846 double begClock = clock();
848 for(int i=startSeq;i<stopSeq;i++){
850 if (m->control_pressed) { break; }
852 for(int j=0;j<i;j++){
853 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
855 if(flowDistance < 1e-6){
856 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
858 else if(flowDistance <= cutoff){
859 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
863 m->mothurOutJustToScreen(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
867 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
868 if(pid != 0){ fDistFileName += ".temp." + toString(pid); }
870 if (m->control_pressed) { return fDistFileName; }
872 m->mothurOutJustToScreen(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
874 ofstream distFile(fDistFileName.c_str());
875 distFile << outStream.str();
878 return fDistFileName;
880 catch(exception& e) {
881 m->errorOut(e, "ShhherCommand", "flowDistMPI");
885 /**************************************************************************************************/
887 void ShhherCommand::getOTUData(string listFileName){
891 m->openInputFile(listFileName, listFile);
894 listFile >> label >> numOTUs;
896 otuData.assign(numSeqs, 0);
897 cumNumSeqs.assign(numOTUs, 0);
898 nSeqsPerOTU.assign(numOTUs, 0);
899 aaP.clear();aaP.resize(numOTUs);
905 string singleOTU = "";
907 for(int i=0;i<numOTUs;i++){
909 if (m->control_pressed) { break; }
911 listFile >> singleOTU;
913 istringstream otuString(singleOTU);
919 for(int j=0;j<singleOTU.length();j++){
920 char letter = otuString.get();
926 map<string,int>::iterator nmIt = nameMap.find(seqName);
927 int index = nmIt->second;
933 aaP[i].push_back(index);
938 map<string,int>::iterator nmIt = nameMap.find(seqName);
940 int index = nmIt->second;
945 aaP[i].push_back(index);
950 sort(aaP[i].begin(), aaP[i].end());
951 for(int j=0;j<nSeqsPerOTU[i];j++){
952 seqNumber.push_back(aaP[i][j]);
954 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
961 for(int i=1;i<numOTUs;i++){
962 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
965 seqIndex = seqNumber;
970 catch(exception& e) {
971 m->errorOut(e, "ShhherCommand", "getOTUData");
976 /**************************************************************************************************/
978 void ShhherCommand::initPyroCluster(){
980 if (numOTUs < processors) { processors = 1; }
982 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
984 dist.assign(numSeqs * numOTUs, 0);
985 change.assign(numOTUs, 1);
986 centroids.assign(numOTUs, -1);
987 weight.assign(numOTUs, 0);
988 singleTau.assign(numSeqs, 1.0);
990 nSeqsBreaks.assign(processors+1, 0);
991 nOTUsBreaks.assign(processors+1, 0);
993 if (m->debug) { m->mothurOut("[DEBUG]: made it through the memory allocation.\n"); }
996 for(int i=0;i<processors;i++){
997 nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
998 nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
1000 nSeqsBreaks[processors] = numSeqs;
1001 nOTUsBreaks[processors] = numOTUs;
1003 catch(exception& e) {
1004 m->errorOut(e, "ShhherCommand", "initPyroCluster");
1009 /**************************************************************************************************/
1011 void ShhherCommand::fill(){
1014 for(int i=0;i<numOTUs;i++){
1016 if (m->control_pressed) { break; }
1018 cumNumSeqs[i] = index;
1019 for(int j=0;j<nSeqsPerOTU[i];j++){
1020 seqNumber[index] = aaP[i][j];
1021 seqIndex[index] = aaI[i][j];
1027 catch(exception& e) {
1028 m->errorOut(e, "ShhherCommand", "fill");
1033 /**************************************************************************************************/
1035 void ShhherCommand::getFlowData(){
1038 m->openInputFile(flowFileName, flowFile);
1041 seqNameVector.clear();
1043 flowDataIntI.clear();
1047 int currentNumFlowCells;
1052 flowFile >> numFlowTest;
1054 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); }
1055 else { convert(numFlowTest, numFlowCells); }
1057 int index = 0;//pcluster
1058 while(!flowFile.eof()){
1060 if (m->control_pressed) { break; }
1062 flowFile >> seqName >> currentNumFlowCells;
1063 lengths.push_back(currentNumFlowCells);
1065 seqNameVector.push_back(seqName);
1066 nameMap[seqName] = index++;//pcluster
1068 for(int i=0;i<numFlowCells;i++){
1069 flowFile >> intensity;
1070 if(intensity > 9.99) { intensity = 9.99; }
1071 int intI = int(100 * intensity + 0.0001);
1072 flowDataIntI.push_back(intI);
1074 m->gobble(flowFile);
1078 numSeqs = seqNameVector.size();
1080 for(int i=0;i<numSeqs;i++){
1082 if (m->control_pressed) { break; }
1084 int iNumFlowCells = i * numFlowCells;
1085 for(int j=lengths[i];j<numFlowCells;j++){
1086 flowDataIntI[iNumFlowCells + j] = 0;
1091 catch(exception& e) {
1092 m->errorOut(e, "ShhherCommand", "getFlowData");
1096 /**************************************************************************************************/
1097 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1100 vector<double> newTau(numOTUs,0);
1101 vector<double> norms(numSeqs, 0);
1106 for(int i=startSeq;i<stopSeq;i++){
1108 if (m->control_pressed) { break; }
1110 double offset = 1e8;
1111 int indexOffset = i * numOTUs;
1113 for(int j=0;j<numOTUs;j++){
1115 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1116 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1118 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1119 offset = dist[indexOffset + j];
1123 for(int j=0;j<numOTUs;j++){
1124 if(weight[j] > MIN_WEIGHT){
1125 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1126 norms[i] += newTau[j];
1133 for(int j=0;j<numOTUs;j++){
1135 newTau[j] /= norms[i];
1137 if(newTau[j] > MIN_TAU){
1138 otuIndex.push_back(j);
1139 seqIndex.push_back(i);
1140 singleTau.push_back(newTau[j]);
1146 catch(exception& e) {
1147 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1152 /**************************************************************************************************/
1154 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1159 vector<double> newTau(numOTUs,0);
1160 vector<double> norms(numSeqs, 0);
1161 nSeqsPerOTU.assign(numOTUs, 0);
1163 for(int i=startSeq;i<stopSeq;i++){
1165 if (m->control_pressed) { break; }
1167 int indexOffset = i * numOTUs;
1169 double offset = 1e8;
1171 for(int j=0;j<numOTUs;j++){
1173 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1174 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1177 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1178 offset = dist[indexOffset + j];
1182 for(int j=0;j<numOTUs;j++){
1183 if(weight[j] > MIN_WEIGHT){
1184 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1185 norms[i] += newTau[j];
1192 for(int j=0;j<numOTUs;j++){
1193 newTau[j] /= norms[i];
1196 for(int j=0;j<numOTUs;j++){
1197 if(newTau[j] > MIN_TAU){
1199 int oldTotal = total;
1203 singleTau.resize(total, 0);
1204 seqNumber.resize(total, 0);
1205 seqIndex.resize(total, 0);
1207 singleTau[oldTotal] = newTau[j];
1209 aaP[j][nSeqsPerOTU[j]] = oldTotal;
1210 aaI[j][nSeqsPerOTU[j]] = i;
1218 catch(exception& e) {
1219 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1224 /**************************************************************************************************/
1226 void ShhherCommand::setOTUs(){
1229 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1231 for(int i=0;i<numOTUs;i++){
1233 if (m->control_pressed) { break; }
1235 for(int j=0;j<nSeqsPerOTU[i];j++){
1236 int index = cumNumSeqs[i] + j;
1237 double tauValue = singleTau[seqNumber[index]];
1238 int sIndex = seqIndex[index];
1239 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1243 for(int i=0;i<numSeqs;i++){
1244 double maxTau = -1.0000;
1246 for(int j=0;j<numOTUs;j++){
1247 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1248 maxTau = bigTauMatrix[i * numOTUs + j];
1253 otuData[i] = maxOTU;
1256 nSeqsPerOTU.assign(numOTUs, 0);
1258 for(int i=0;i<numSeqs;i++){
1259 int index = otuData[i];
1261 singleTau[i] = 1.0000;
1264 aaP[index][nSeqsPerOTU[index]] = i;
1265 aaI[index][nSeqsPerOTU[index]] = i;
1267 nSeqsPerOTU[index]++;
1271 catch(exception& e) {
1272 m->errorOut(e, "ShhherCommand", "setOTUs");
1277 /**************************************************************************************************/
1279 void ShhherCommand::getUniques(){
1284 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
1285 uniqueCount.assign(numSeqs, 0); // anWeights
1286 uniqueLengths.assign(numSeqs, 0);
1287 mapSeqToUnique.assign(numSeqs, -1);
1288 mapUniqueToSeq.assign(numSeqs, -1);
1290 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
1292 for(int i=0;i<numSeqs;i++){
1294 if (m->control_pressed) { break; }
1298 vector<short> current(numFlowCells);
1299 for(int j=0;j<numFlowCells;j++){
1300 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
1303 for(int j=0;j<numUniques;j++){
1304 int offset = j * numFlowCells;
1308 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
1309 else { shorterLength = uniqueLengths[j]; }
1311 for(int k=0;k<shorterLength;k++){
1312 if(current[k] != uniqueFlowgrams[offset + k]){
1319 mapSeqToUnique[i] = j;
1322 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
1328 if(index == numUniques){
1329 uniqueLengths[numUniques] = lengths[i];
1330 uniqueCount[numUniques] = 1;
1331 mapSeqToUnique[i] = numUniques;//anMap
1332 mapUniqueToSeq[numUniques] = i;//anF
1334 for(int k=0;k<numFlowCells;k++){
1335 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
1336 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
1342 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
1343 uniqueLengths.resize(numUniques);
1345 flowDataPrI.resize(numSeqs * numFlowCells, 0);
1346 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
1348 catch(exception& e) {
1349 m->errorOut(e, "ShhherCommand", "getUniques");
1354 /**************************************************************************************************/
1356 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
1358 int minLength = lengths[mapSeqToUnique[seqA]];
1359 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
1361 int ANumFlowCells = seqA * numFlowCells;
1362 int BNumFlowCells = seqB * numFlowCells;
1366 for(int i=0;i<minLength;i++){
1368 if (m->control_pressed) { break; }
1370 int flowAIntI = flowDataIntI[ANumFlowCells + i];
1371 float flowAPrI = flowDataPrI[ANumFlowCells + i];
1373 int flowBIntI = flowDataIntI[BNumFlowCells + i];
1374 float flowBPrI = flowDataPrI[BNumFlowCells + i];
1375 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
1378 dist /= (float) minLength;
1381 catch(exception& e) {
1382 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
1387 //**********************************************************************************************************************/
1389 string ShhherCommand::cluster(string distFileName, string namesFileName){
1392 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1393 read->setCutoff(cutoff);
1395 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1396 clusterNameMap->readMap();
1397 read->read(clusterNameMap);
1399 ListVector* list = read->getListVector();
1400 SparseDistanceMatrix* matrix = read->getDMatrix();
1403 delete clusterNameMap;
1405 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1407 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1408 string tag = cluster->getTag();
1410 double clusterCutoff = cutoff;
1411 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1413 if (m->control_pressed) { break; }
1415 cluster->update(clusterCutoff);
1418 list->setLabel(toString(cutoff));
1420 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
1422 m->openOutputFile(listFileName, listFile);
1423 list->print(listFile);
1426 delete matrix; delete cluster; delete rabund; delete list;
1428 return listFileName;
1430 catch(exception& e) {
1431 m->errorOut(e, "ShhherCommand", "cluster");
1436 /**************************************************************************************************/
1438 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1440 //this function gets the most likely homopolymer length at a flow position for a group of sequences
1445 for(int i=start;i<finish;i++){
1447 if (m->control_pressed) { break; }
1451 int minFlowGram = 100000000;
1452 double minFlowValue = 1e8;
1453 change[i] = 0; //FALSE
1455 for(int j=0;j<nSeqsPerOTU[i];j++){
1456 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1459 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1460 vector<double> adF(nSeqsPerOTU[i]);
1461 vector<int> anL(nSeqsPerOTU[i]);
1463 for(int j=0;j<nSeqsPerOTU[i];j++){
1464 int index = cumNumSeqs[i] + j;
1465 int nI = seqIndex[index];
1466 int nIU = mapSeqToUnique[nI];
1469 for(k=0;k<position;k++){
1475 anL[position] = nIU;
1476 adF[position] = 0.0000;
1481 for(int j=0;j<nSeqsPerOTU[i];j++){
1482 int index = cumNumSeqs[i] + j;
1483 int nI = seqIndex[index];
1485 double tauValue = singleTau[seqNumber[index]];
1487 for(int k=0;k<position;k++){
1488 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1489 adF[k] += dist * tauValue;
1493 for(int j=0;j<position;j++){
1494 if(adF[j] < minFlowValue){
1496 minFlowValue = adF[j];
1500 if(centroids[i] != anL[minFlowGram]){
1502 centroids[i] = anL[minFlowGram];
1505 else if(centroids[i] != -1){
1511 catch(exception& e) {
1512 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1517 /**************************************************************************************************/
1519 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1522 int flowAValue = cent * numFlowCells;
1523 int flowBValue = flow * numFlowCells;
1527 for(int i=0;i<length;i++){
1528 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1533 return dist / (double)length;
1535 catch(exception& e) {
1536 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1541 /**************************************************************************************************/
1543 double ShhherCommand::getNewWeights(){
1546 double maxChange = 0;
1548 for(int i=0;i<numOTUs;i++){
1550 if (m->control_pressed) { break; }
1552 double difference = weight[i];
1555 for(int j=0;j<nSeqsPerOTU[i];j++){
1556 int index = cumNumSeqs[i] + j;
1557 double tauValue = singleTau[seqNumber[index]];
1558 weight[i] += tauValue;
1561 difference = fabs(weight[i] - difference);
1562 if(difference > maxChange){ maxChange = difference; }
1566 catch(exception& e) {
1567 m->errorOut(e, "ShhherCommand", "getNewWeights");
1572 /**************************************************************************************************/
1574 double ShhherCommand::getLikelihood(){
1578 vector<long double> P(numSeqs, 0);
1581 for(int i=0;i<numOTUs;i++){
1582 if(weight[i] > MIN_WEIGHT){
1588 for(int i=0;i<numOTUs;i++){
1590 if (m->control_pressed) { break; }
1592 for(int j=0;j<nSeqsPerOTU[i];j++){
1593 int index = cumNumSeqs[i] + j;
1594 int nI = seqIndex[index];
1595 double singleDist = dist[seqNumber[index]];
1597 P[nI] += weight[i] * exp(-singleDist * sigma);
1601 for(int i=0;i<numSeqs;i++){
1602 if(P[i] == 0){ P[i] = DBL_EPSILON; }
1607 nLL = nLL -(double)numSeqs * log(sigma);
1611 catch(exception& e) {
1612 m->errorOut(e, "ShhherCommand", "getNewWeights");
1617 /**************************************************************************************************/
1619 void ShhherCommand::checkCentroids(){
1621 vector<int> unique(numOTUs, 1);
1623 for(int i=0;i<numOTUs;i++){
1624 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1629 for(int i=0;i<numOTUs;i++){
1631 if (m->control_pressed) { break; }
1634 for(int j=i+1;j<numOTUs;j++){
1637 if(centroids[j] == centroids[i]){
1641 weight[i] += weight[j];
1649 catch(exception& e) {
1650 m->errorOut(e, "ShhherCommand", "checkCentroids");
1654 /**************************************************************************************************/
1658 void ShhherCommand::writeQualities(vector<int> otuCounts){
1661 string thisOutputDir = outputDir;
1662 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1663 map<string, string> variables;
1664 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(flowFileName));
1665 string qualityFileName = getOutputFileName("qfile",variables);
1667 ofstream qualityFile;
1668 m->openOutputFile(qualityFileName, qualityFile);
1670 qualityFile.setf(ios::fixed, ios::floatfield);
1671 qualityFile.setf(ios::showpoint);
1672 qualityFile << setprecision(6);
1674 vector<vector<int> > qualities(numOTUs);
1675 vector<double> pr(HOMOPS, 0);
1678 for(int i=0;i<numOTUs;i++){
1680 if (m->control_pressed) { break; }
1685 if(nSeqsPerOTU[i] > 0){
1686 qualities[i].assign(1024, -1);
1688 while(index < numFlowCells){
1689 double maxPrValue = 1e8;
1690 short maxPrIndex = -1;
1691 double count = 0.0000;
1693 pr.assign(HOMOPS, 0);
1695 for(int j=0;j<nSeqsPerOTU[i];j++){
1696 int lIndex = cumNumSeqs[i] + j;
1697 double tauValue = singleTau[seqNumber[lIndex]];
1698 int sequenceIndex = aaI[i][j];
1699 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1703 for(int s=0;s<HOMOPS;s++){
1704 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1708 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1709 maxPrValue = pr[maxPrIndex];
1711 if(count > MIN_COUNT){
1713 double norm = 0.0000;
1715 for(int s=0;s<HOMOPS;s++){
1716 norm += exp(-(pr[s] - maxPrValue));
1719 for(int s=1;s<=maxPrIndex;s++){
1721 double temp = 0.0000;
1723 U += exp(-(pr[s-1]-maxPrValue))/norm;
1731 temp = floor(-10 * temp);
1732 value = (int)floor(temp);
1733 if(value > 100){ value = 100; }
1735 qualities[i][base] = (int)value;
1745 if(otuCounts[i] > 0){
1746 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1748 int j=4; //need to get past the first four bases
1749 while(qualities[i][j] != -1){
1750 qualityFile << qualities[i][j] << ' ';
1753 qualityFile << endl;
1756 qualityFile.close();
1757 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
1760 catch(exception& e) {
1761 m->errorOut(e, "ShhherCommand", "writeQualities");
1766 /**************************************************************************************************/
1768 void ShhherCommand::writeSequences(vector<int> otuCounts){
1770 string thisOutputDir = outputDir;
1771 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1772 map<string, string> variables;
1773 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1774 string fastaFileName = getOutputFileName("fasta",variables);
1776 m->openOutputFile(fastaFileName, fastaFile);
1778 vector<string> names(numOTUs, "");
1780 for(int i=0;i<numOTUs;i++){
1782 if (m->control_pressed) { break; }
1784 int index = centroids[i];
1786 if(otuCounts[i] > 0){
1787 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
1791 for(int j=0;j<numFlowCells;j++){
1793 char base = flowOrder[j % flowOrder.length()];
1794 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
1799 fastaFile << newSeq.substr(4) << endl;
1804 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
1806 if(compositeFASTAFileName != ""){
1807 m->appendFiles(fastaFileName, compositeFASTAFileName);
1810 catch(exception& e) {
1811 m->errorOut(e, "ShhherCommand", "writeSequences");
1816 /**************************************************************************************************/
1818 void ShhherCommand::writeNames(vector<int> otuCounts){
1820 string thisOutputDir = outputDir;
1821 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1822 map<string, string> variables;
1823 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1824 string nameFileName = getOutputFileName("name",variables);
1826 m->openOutputFile(nameFileName, nameFile);
1828 for(int i=0;i<numOTUs;i++){
1830 if (m->control_pressed) { break; }
1832 if(otuCounts[i] > 0){
1833 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
1835 for(int j=1;j<nSeqsPerOTU[i];j++){
1836 nameFile << ',' << seqNameVector[aaI[i][j]];
1843 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
1846 if(compositeNamesFileName != ""){
1847 m->appendFiles(nameFileName, compositeNamesFileName);
1850 catch(exception& e) {
1851 m->errorOut(e, "ShhherCommand", "writeNames");
1856 /**************************************************************************************************/
1858 void ShhherCommand::writeGroups(){
1860 string thisOutputDir = outputDir;
1861 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1862 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
1863 int pos = fileRoot.find_first_of('.');
1864 string fileGroup = fileRoot;
1865 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
1866 map<string, string> variables;
1867 variables["[filename]"] = thisOutputDir + fileRoot;
1868 string groupFileName = getOutputFileName("group",variables);
1870 m->openOutputFile(groupFileName, groupFile);
1872 for(int i=0;i<numSeqs;i++){
1873 if (m->control_pressed) { break; }
1874 groupFile << seqNameVector[i] << '\t' << fileGroup << endl;
1877 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
1880 catch(exception& e) {
1881 m->errorOut(e, "ShhherCommand", "writeGroups");
1886 /**************************************************************************************************/
1888 void ShhherCommand::writeClusters(vector<int> otuCounts){
1890 string thisOutputDir = outputDir;
1891 if (outputDir == "") { thisOutputDir += m->hasPath(flowFileName); }
1892 map<string, string> variables;
1893 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
1894 string otuCountsFileName = getOutputFileName("counts",variables);
1895 ofstream otuCountsFile;
1896 m->openOutputFile(otuCountsFileName, otuCountsFile);
1898 string bases = flowOrder;
1900 for(int i=0;i<numOTUs;i++){
1902 if (m->control_pressed) {
1905 //output the translated version of the centroid sequence for the otu
1906 if(otuCounts[i] > 0){
1907 int index = centroids[i];
1909 otuCountsFile << "ideal\t";
1910 for(int j=8;j<numFlowCells;j++){
1911 char base = bases[j % bases.length()];
1912 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
1913 otuCountsFile << base;
1916 otuCountsFile << endl;
1918 for(int j=0;j<nSeqsPerOTU[i];j++){
1919 int sequence = aaI[i][j];
1920 otuCountsFile << seqNameVector[sequence] << '\t';
1924 for(int k=0;k<lengths[sequence];k++){
1925 char base = bases[k % bases.length()];
1926 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
1928 for(int s=0;s<freq;s++){
1930 //otuCountsFile << base;
1933 otuCountsFile << newSeq.substr(4) << endl;
1935 otuCountsFile << endl;
1938 otuCountsFile.close();
1939 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
1942 catch(exception& e) {
1943 m->errorOut(e, "ShhherCommand", "writeClusters");
1949 //**********************************************************************************************************************
1951 int ShhherCommand::execute(){
1953 if (abort == true) { if (calledHelp) { return 0; } return 2; }
1955 getSingleLookUp(); if (m->control_pressed) { return 0; }
1956 getJointLookUp(); if (m->control_pressed) { return 0; }
1958 int numFiles = flowFileVector.size();
1960 if (numFiles < processors) { processors = numFiles; }
1962 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1963 if (processors == 1) { driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName); }
1964 else { createProcesses(flowFileVector); } //each processor processes one file
1966 driver(flowFileVector, compositeFASTAFileName, compositeNamesFileName);
1969 if(compositeFASTAFileName != ""){
1970 outputNames.push_back(compositeFASTAFileName); outputTypes["fasta"].push_back(compositeFASTAFileName);
1971 outputNames.push_back(compositeNamesFileName); outputTypes["name"].push_back(compositeNamesFileName);
1974 m->mothurOutEndLine();
1975 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
1976 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
1977 m->mothurOutEndLine();
1981 catch(exception& e) {
1982 m->errorOut(e, "ShhherCommand", "execute");
1987 //********************************************************************************************************************
1988 //sorts biggest to smallest
1989 inline bool compareFileSizes(string left, string right){
1994 //get num bytes in file
1995 string filename = left;
1996 pFile = fopen (filename.c_str(),"rb");
1997 string error = "Error opening " + filename;
1998 if (pFile==NULL) perror (error.c_str());
2000 fseek (pFile, 0, SEEK_END);
2001 leftsize=ftell (pFile);
2008 //get num bytes in file
2010 pFile2 = fopen (filename.c_str(),"rb");
2011 error = "Error opening " + filename;
2012 if (pFile2==NULL) perror (error.c_str());
2014 fseek (pFile2, 0, SEEK_END);
2015 rightsize=ftell (pFile2);
2019 return (leftsize > rightsize);
2021 /**************************************************************************************************/
2023 int ShhherCommand::createProcesses(vector<string> filenames){
2025 vector<int> processIDS;
2030 if (filenames.size() < processors) { processors = filenames.size(); }
2032 //sort file names by size to divide load better
2033 sort(filenames.begin(), filenames.end(), compareFileSizes);
2035 vector < vector <string> > dividedFiles; //dividedFiles[1] = vector of filenames for process 1...
2036 dividedFiles.resize(processors);
2038 //for each file, figure out which process will complete it
2039 //want to divide the load intelligently so the big files are spread between processes
2040 for (int i = 0; i < filenames.size(); i++) {
2041 int processToAssign = (i+1) % processors;
2042 if (processToAssign == 0) { processToAssign = processors; }
2044 dividedFiles[(processToAssign-1)].push_back(filenames[i]);
2047 //now lets reverse the order of ever other process, so we balance big files running with little ones
2048 for (int i = 0; i < processors; i++) {
2049 int remainder = ((i+1) % processors);
2050 if (remainder) { reverse(dividedFiles[i].begin(), dividedFiles[i].end()); }
2054 //divide the groups between the processors
2055 /*vector<linePair> lines;
2056 vector<int> numFilesToComplete;
2057 int numFilesPerProcessor = filenames.size() / processors;
2058 for (int i = 0; i < processors; i++) {
2059 int startIndex = i * numFilesPerProcessor;
2060 int endIndex = (i+1) * numFilesPerProcessor;
2061 if(i == (processors - 1)){ endIndex = filenames.size(); }
2062 lines.push_back(linePair(startIndex, endIndex));
2063 numFilesToComplete.push_back((endIndex-startIndex));
2066 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
2068 //loop through and create all the processes you want
2069 while (process != processors) {
2073 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
2075 }else if (pid == 0){
2076 num = driver(dividedFiles[process], compositeFASTAFileName + toString(getpid()) + ".temp", compositeNamesFileName + toString(getpid()) + ".temp");
2078 //pass numSeqs to parent
2080 string tempFile = compositeFASTAFileName + toString(getpid()) + ".num.temp";
2081 m->openOutputFile(tempFile, out);
2087 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
2088 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
2094 driver(dividedFiles[0], compositeFASTAFileName, compositeNamesFileName);
2096 //force parent to wait until all the processes are done
2097 for (int i=0;i<processIDS.size();i++) {
2098 int temp = processIDS[i];
2104 //////////////////////////////////////////////////////////////////////////////////////////////////////
2106 /////////////////////// NOT WORKING, ACCESS VIOLATION ON READ OF FLOWGRAMS IN THREAD /////////////////
2108 //////////////////////////////////////////////////////////////////////////////////////////////////////
2109 //Windows version shared memory, so be careful when passing variables through the shhhFlowsData struct.
2110 //Above fork() will clone, so memory is separate, but that's not the case with windows,
2111 //////////////////////////////////////////////////////////////////////////////////////////////////////
2113 vector<shhhFlowsData*> pDataArray;
2114 DWORD dwThreadIdArray[processors-1];
2115 HANDLE hThreadArray[processors-1];
2117 //Create processor worker threads.
2118 for( int i=0; i<processors-1; i++ ){
2119 // Allocate memory for thread data.
2120 string extension = "";
2121 if (i != 0) { extension = toString(i) + ".temp"; }
2123 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);
2124 pDataArray.push_back(tempFlow);
2125 processIDS.push_back(i);
2127 hThreadArray[i] = CreateThread(NULL, 0, ShhhFlowsThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
2130 //using the main process as a worker saves time and memory
2132 driver(filenames, compositeFASTAFileName, compositeNamesFileName, lines[processors-1].start, lines[processors-1].end);
2134 //Wait until all threads have terminated.
2135 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
2137 //Close all thread handles and free memory allocations.
2138 for(int i=0; i < pDataArray.size(); i++){
2139 for(int j=0; j < pDataArray[i]->outputNames.size(); j++){ outputNames.push_back(pDataArray[i]->outputNames[j]); }
2140 CloseHandle(hThreadArray[i]);
2141 delete pDataArray[i];
2146 for (int i=0;i<processIDS.size();i++) {
2148 string tempFile = compositeFASTAFileName + toString(processIDS[i]) + ".num.temp";
2149 m->openInputFile(tempFile, in);
2153 if (tempNum != dividedFiles[i+1].size()) {
2154 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");
2157 in.close(); m->mothurRemove(tempFile);
2159 if (compositeFASTAFileName != "") {
2160 m->appendFiles((compositeFASTAFileName + toString(processIDS[i]) + ".temp"), compositeFASTAFileName);
2161 m->appendFiles((compositeNamesFileName + toString(processIDS[i]) + ".temp"), compositeNamesFileName);
2162 m->mothurRemove((compositeFASTAFileName + toString(processIDS[i]) + ".temp"));
2163 m->mothurRemove((compositeNamesFileName + toString(processIDS[i]) + ".temp"));
2170 catch(exception& e) {
2171 m->errorOut(e, "ShhherCommand", "createProcesses");
2175 /**************************************************************************************************/
2177 vector<string> ShhherCommand::parseFlowFiles(string filename){
2179 vector<string> files;
2183 m->openInputFile(filename, in);
2185 int thisNumFLows = 0;
2186 in >> thisNumFLows; m->gobble(in);
2189 if (m->control_pressed) { break; }
2192 string outputFileName = filename + toString(count) + ".temp";
2193 m->openOutputFile(outputFileName, out);
2194 out << thisNumFLows << endl;
2195 files.push_back(outputFileName);
2197 int numLinesWrote = 0;
2198 for (int i = 0; i < largeSize; i++) {
2199 if (in.eof()) { break; }
2200 string line = m->getline(in); m->gobble(in);
2201 out << line << endl;
2206 if (numLinesWrote == 0) { m->mothurRemove(outputFileName); files.pop_back(); }
2211 if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } files.clear(); }
2213 m->mothurOut("\nDivided " + filename + " into " + toString(files.size()) + " files.\n\n");
2217 catch(exception& e) {
2218 m->errorOut(e, "ShhherCommand", "parseFlowFiles");
2222 /**************************************************************************************************/
2224 int ShhherCommand::driver(vector<string> filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName){
2227 int numCompleted = 0;
2229 for(int i=0;i<filenames.size();i++){
2231 if (m->control_pressed) { break; }
2233 vector<string> theseFlowFileNames; theseFlowFileNames.push_back(filenames[i]);
2234 if (large) { theseFlowFileNames = parseFlowFiles(filenames[i]); }
2236 if (m->control_pressed) { break; }
2238 double begClock = clock();
2239 unsigned long long begTime;
2241 string fileNameForOutput = filenames[i];
2243 for (int g = 0; g < theseFlowFileNames.size(); g++) {
2245 string flowFileName = theseFlowFileNames[g];
2246 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n");
2247 m->mothurOut("Reading flowgrams...\n");
2249 vector<string> seqNameVector;
2250 vector<int> lengths;
2251 vector<short> flowDataIntI;
2252 vector<double> flowDataPrI;
2253 map<string, int> nameMap;
2254 vector<short> uniqueFlowgrams;
2255 vector<int> uniqueCount;
2256 vector<int> mapSeqToUnique;
2257 vector<int> mapUniqueToSeq;
2258 vector<int> uniqueLengths;
2261 if (m->debug) { m->mothurOut("[DEBUG]: About to read flowgrams.\n"); }
2262 int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells);
2264 if (m->control_pressed) { break; }
2266 m->mothurOut("Identifying unique flowgrams...\n");
2267 int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI);
2269 if (m->control_pressed) { break; }
2271 m->mothurOut("Calculating distances between flowgrams...\n");
2272 string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist";
2273 begTime = time(NULL);
2276 flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2278 m->mothurOutEndLine();
2279 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
2282 string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names";
2283 createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq);
2285 if (m->control_pressed) { break; }
2287 m->mothurOut("\nClustering flowgrams...\n");
2288 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list";
2289 cluster(listFileName, distFileName, namesFileName);
2291 if (m->control_pressed) { break; }
2293 vector<int> otuData;
2294 vector<int> cumNumSeqs;
2295 vector<int> nSeqsPerOTU;
2296 vector<vector<int> > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2297 vector<vector<int> > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2298 vector<int> seqNumber; //tMaster->anP: the sequence id number sorted by OTU
2299 vector<int> seqIndex; //tMaster->anI; the index that corresponds to seqNumber
2302 int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap);
2304 if (m->control_pressed) { break; }
2306 m->mothurRemove(distFileName);
2307 m->mothurRemove(namesFileName);
2308 m->mothurRemove(listFileName);
2310 vector<double> dist; //adDist - distance of sequences to centroids
2311 vector<short> change; //did the centroid sequence change? 0 = no; 1 = yes
2312 vector<int> centroids; //the representative flowgram for each cluster m
2313 vector<double> weight;
2314 vector<double> singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2315 vector<int> nSeqsBreaks;
2316 vector<int> nOTUsBreaks;
2318 if (m->debug) { m->mothurOut("[DEBUG]: numSeqs = " + toString(numSeqs) + " numOTUS = " + toString(numOTUs) + " about to alloc a dist vector with size = " + toString((numSeqs * numOTUs)) + ".\n"); }
2320 dist.assign(numSeqs * numOTUs, 0);
2321 change.assign(numOTUs, 1);
2322 centroids.assign(numOTUs, -1);
2323 weight.assign(numOTUs, 0);
2324 singleTau.assign(numSeqs, 1.0);
2326 nSeqsBreaks.assign(2, 0);
2327 nOTUsBreaks.assign(2, 0);
2330 nSeqsBreaks[1] = numSeqs;
2331 nOTUsBreaks[1] = numOTUs;
2333 if (m->debug) { m->mothurOut("[DEBUG]: done allocating memory, about to denoise.\n"); }
2335 if (m->control_pressed) { break; }
2337 double maxDelta = 0;
2341 begTime = time(NULL);
2343 m->mothurOut("\nDenoising flowgrams...\n");
2344 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
2346 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
2348 if (m->control_pressed) { break; }
2350 double cycClock = clock();
2351 unsigned long long cycTime = time(NULL);
2352 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2354 if (m->control_pressed) { break; }
2356 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2358 if (m->control_pressed) { break; }
2360 maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight);
2362 if (m->control_pressed) { break; }
2364 double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight);
2366 if (m->control_pressed) { break; }
2368 checkCentroids(numOTUs, centroids, weight);
2370 if (m->control_pressed) { break; }
2372 calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths);
2374 if (m->control_pressed) { break; }
2378 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
2382 if (m->control_pressed) { break; }
2384 m->mothurOut("\nFinalizing...\n");
2385 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
2387 if (m->control_pressed) { break; }
2389 setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI);
2391 if (m->control_pressed) { break; }
2393 vector<int> otuCounts(numOTUs, 0);
2394 for(int j=0;j<numSeqs;j++) { otuCounts[otuData[j]]++; }
2396 calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber);
2398 if (m->control_pressed) { break; }
2400 if ((large) && (g == 0)) { flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; }
2401 string thisOutputDir = outputDir;
2402 if (outputDir == "") { thisOutputDir = m->hasPath(flowFileName); }
2403 map<string, string> variables;
2404 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName));
2405 string qualityFileName = getOutputFileName("qfile",variables);
2406 string fastaFileName = getOutputFileName("fasta",variables);
2407 string nameFileName = getOutputFileName("name",variables);
2408 string otuCountsFileName = getOutputFileName("counts",variables);
2409 string fileRoot = m->getRootName(m->getSimpleName(flowFileName));
2410 int pos = fileRoot.find_first_of('.');
2411 string fileGroup = fileRoot;
2412 if (pos != string::npos) { fileGroup = fileRoot.substr(pos+1, (fileRoot.length()-1-(pos+1))); }
2413 string groupFileName = getOutputFileName("group",variables);
2416 writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; }
2417 writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; }
2418 writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; }
2419 writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; }
2420 writeGroups(groupFileName, fileGroup, numSeqs, seqNameVector); if (m->control_pressed) { break; }
2424 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0]));
2425 m->appendFiles(qualityFileName, getOutputFileName("qfile",variables));
2426 m->mothurRemove(qualityFileName);
2427 m->appendFiles(fastaFileName, getOutputFileName("fasta",variables));
2428 m->mothurRemove(fastaFileName);
2429 m->appendFiles(nameFileName, getOutputFileName("name",variables));
2430 m->mothurRemove(nameFileName);
2431 m->appendFiles(otuCountsFileName, getOutputFileName("counts",variables));
2432 m->mothurRemove(otuCountsFileName);
2433 m->appendFiles(groupFileName, getOutputFileName("group",variables));
2434 m->mothurRemove(groupFileName);
2436 m->mothurRemove(theseFlowFileNames[g]);
2441 m->mothurOut("Total time to process " + fileNameForOutput + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
2444 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
2446 return numCompleted;
2448 }catch(exception& e) {
2449 m->errorOut(e, "ShhherCommand", "driver");
2454 /**************************************************************************************************/
2455 int ShhherCommand::getFlowData(string filename, vector<string>& thisSeqNameVector, vector<int>& thisLengths, vector<short>& thisFlowDataIntI, map<string, int>& thisNameMap, int& numFlowCells){
2460 m->openInputFile(filename, flowFile);
2463 int currentNumFlowCells;
2465 thisSeqNameVector.clear();
2466 thisLengths.clear();
2467 thisFlowDataIntI.clear();
2468 thisNameMap.clear();
2471 flowFile >> numFlowTest;
2473 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); }
2474 else { convert(numFlowTest, numFlowCells); }
2476 if (m->debug) { m->mothurOut("[DEBUG]: numFlowCells = " + toString(numFlowCells) + ".\n"); }
2477 int index = 0;//pcluster
2478 while(!flowFile.eof()){
2480 if (m->control_pressed) { break; }
2482 flowFile >> seqName >> currentNumFlowCells;
2484 thisLengths.push_back(currentNumFlowCells);
2486 thisSeqNameVector.push_back(seqName);
2487 thisNameMap[seqName] = index++;//pcluster
2489 if (m->debug) { m->mothurOut("[DEBUG]: seqName = " + seqName + " length = " + toString(currentNumFlowCells) + " index = " + toString(index) + "\n"); }
2491 for(int i=0;i<numFlowCells;i++){
2492 flowFile >> intensity;
2493 if(intensity > 9.99) { intensity = 9.99; }
2494 int intI = int(100 * intensity + 0.0001);
2495 thisFlowDataIntI.push_back(intI);
2497 m->gobble(flowFile);
2501 int numSeqs = thisSeqNameVector.size();
2503 for(int i=0;i<numSeqs;i++){
2505 if (m->control_pressed) { break; }
2507 int iNumFlowCells = i * numFlowCells;
2508 for(int j=thisLengths[i];j<numFlowCells;j++){
2509 thisFlowDataIntI[iNumFlowCells + j] = 0;
2516 catch(exception& e) {
2517 m->errorOut(e, "ShhherCommand", "getFlowData");
2521 /**************************************************************************************************/
2523 int ShhherCommand::flowDistParentFork(int numFlowCells, string distFileName, int stopSeq, vector<int>& mapUniqueToSeq, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2526 ostringstream outStream;
2527 outStream.setf(ios::fixed, ios::floatfield);
2528 outStream.setf(ios::dec, ios::basefield);
2529 outStream.setf(ios::showpoint);
2530 outStream.precision(6);
2532 int begTime = time(NULL);
2533 double begClock = clock();
2535 for(int i=0;i<stopSeq;i++){
2537 if (m->control_pressed) { break; }
2539 for(int j=0;j<i;j++){
2540 float flowDistance = calcPairwiseDist(numFlowCells, mapUniqueToSeq[i], mapUniqueToSeq[j], mapSeqToUnique, lengths, flowDataPrI, flowDataIntI);
2542 if(flowDistance < 1e-6){
2543 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
2545 else if(flowDistance <= cutoff){
2546 outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
2550 m->mothurOutJustToScreen(toString(i) + "\t" + toString(time(NULL) - begTime));
2551 m->mothurOutJustToScreen("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)+"\n");
2555 ofstream distFile(distFileName.c_str());
2556 distFile << outStream.str();
2559 if (m->control_pressed) {}
2561 m->mothurOutJustToScreen(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
2562 m->mothurOutJustToScreen("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC)+"\n");
2567 catch(exception& e) {
2568 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
2572 /**************************************************************************************************/
2574 float ShhherCommand::calcPairwiseDist(int numFlowCells, int seqA, int seqB, vector<int>& mapSeqToUnique, vector<int>& lengths, vector<double>& flowDataPrI, vector<short>& flowDataIntI){
2576 int minLength = lengths[mapSeqToUnique[seqA]];
2577 if(lengths[seqB] < minLength){ minLength = lengths[mapSeqToUnique[seqB]]; }
2579 int ANumFlowCells = seqA * numFlowCells;
2580 int BNumFlowCells = seqB * numFlowCells;
2584 for(int i=0;i<minLength;i++){
2586 if (m->control_pressed) { break; }
2588 int flowAIntI = flowDataIntI[ANumFlowCells + i];
2589 float flowAPrI = flowDataPrI[ANumFlowCells + i];
2591 int flowBIntI = flowDataIntI[BNumFlowCells + i];
2592 float flowBPrI = flowDataPrI[BNumFlowCells + i];
2593 dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
2596 dist /= (float) minLength;
2599 catch(exception& e) {
2600 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
2605 /**************************************************************************************************/
2607 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){
2610 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
2611 uniqueCount.assign(numSeqs, 0); // anWeights
2612 uniqueLengths.assign(numSeqs, 0);
2613 mapSeqToUnique.assign(numSeqs, -1);
2614 mapUniqueToSeq.assign(numSeqs, -1);
2616 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
2618 for(int i=0;i<numSeqs;i++){
2620 if (m->control_pressed) { break; }
2624 vector<short> current(numFlowCells);
2625 for(int j=0;j<numFlowCells;j++){
2626 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
2629 for(int j=0;j<numUniques;j++){
2630 int offset = j * numFlowCells;
2634 if(lengths[i] < uniqueLengths[j]) { shorterLength = lengths[i]; }
2635 else { shorterLength = uniqueLengths[j]; }
2637 for(int k=0;k<shorterLength;k++){
2638 if(current[k] != uniqueFlowgrams[offset + k]){
2645 mapSeqToUnique[i] = j;
2648 if(lengths[i] > uniqueLengths[j]) { uniqueLengths[j] = lengths[i]; }
2654 if(index == numUniques){
2655 uniqueLengths[numUniques] = lengths[i];
2656 uniqueCount[numUniques] = 1;
2657 mapSeqToUnique[i] = numUniques;//anMap
2658 mapUniqueToSeq[numUniques] = i;//anF
2660 for(int k=0;k<numFlowCells;k++){
2661 uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
2662 uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
2668 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
2669 uniqueLengths.resize(numUniques);
2671 flowDataPrI.resize(numSeqs * numFlowCells, 0);
2672 for(int i=0;i<flowDataPrI.size();i++) { if (m->control_pressed) { break; } flowDataPrI[i] = getProbIntensity(flowDataIntI[i]); }
2676 catch(exception& e) {
2677 m->errorOut(e, "ShhherCommand", "getUniques");
2681 /**************************************************************************************************/
2682 int ShhherCommand::createNamesFile(int numSeqs, int numUniques, string filename, vector<string>& seqNameVector, vector<int>& mapSeqToUnique, vector<int>& mapUniqueToSeq){
2685 vector<string> duplicateNames(numUniques, "");
2686 for(int i=0;i<numSeqs;i++){
2687 duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
2691 m->openOutputFile(filename, nameFile);
2693 for(int i=0;i<numUniques;i++){
2695 if (m->control_pressed) { break; }
2697 // nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2698 nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
2705 catch(exception& e) {
2706 m->errorOut(e, "ShhherCommand", "createNamesFile");
2710 //**********************************************************************************************************************
2712 int ShhherCommand::cluster(string filename, string distFileName, string namesFileName){
2715 ReadMatrix* read = new ReadColumnMatrix(distFileName);
2716 read->setCutoff(cutoff);
2718 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
2719 clusterNameMap->readMap();
2720 read->read(clusterNameMap);
2722 ListVector* list = read->getListVector();
2723 SparseDistanceMatrix* matrix = read->getDMatrix();
2726 delete clusterNameMap;
2728 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
2730 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
2731 string tag = cluster->getTag();
2733 double clusterCutoff = cutoff;
2734 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
2736 if (m->control_pressed) { break; }
2738 cluster->update(clusterCutoff);
2741 list->setLabel(toString(cutoff));
2744 m->openOutputFile(filename, listFile);
2745 list->print(listFile);
2748 delete matrix; delete cluster; delete rabund; delete list;
2752 catch(exception& e) {
2753 m->errorOut(e, "ShhherCommand", "cluster");
2757 /**************************************************************************************************/
2759 int ShhherCommand::getOTUData(int numSeqs, string fileName, vector<int>& otuData,
2760 vector<int>& cumNumSeqs,
2761 vector<int>& nSeqsPerOTU,
2762 vector<vector<int> >& aaP, //tMaster->aanP: each row is a different otu / each col contains the sequence indices
2763 vector<vector<int> >& aaI, //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI
2764 vector<int>& seqNumber, //tMaster->anP: the sequence id number sorted by OTU
2765 vector<int>& seqIndex,
2766 map<string, int>& nameMap){
2770 m->openInputFile(fileName, listFile);
2774 listFile >> label >> numOTUs;
2776 if (m->debug) { m->mothurOut("[DEBUG]: Getting OTU Data...\n"); }
2778 otuData.assign(numSeqs, 0);
2779 cumNumSeqs.assign(numOTUs, 0);
2780 nSeqsPerOTU.assign(numOTUs, 0);
2781 aaP.clear();aaP.resize(numOTUs);
2787 string singleOTU = "";
2789 for(int i=0;i<numOTUs;i++){
2791 if (m->control_pressed) { break; }
2792 if (m->debug) { m->mothurOut("[DEBUG]: processing OTU " + toString(i) + ".\n"); }
2794 listFile >> singleOTU;
2796 istringstream otuString(singleOTU);
2800 string seqName = "";
2802 for(int j=0;j<singleOTU.length();j++){
2803 char letter = otuString.get();
2809 map<string,int>::iterator nmIt = nameMap.find(seqName);
2810 int index = nmIt->second;
2812 nameMap.erase(nmIt);
2816 aaP[i].push_back(index);
2821 map<string,int>::iterator nmIt = nameMap.find(seqName);
2823 int index = nmIt->second;
2824 nameMap.erase(nmIt);
2828 aaP[i].push_back(index);
2833 sort(aaP[i].begin(), aaP[i].end());
2834 for(int j=0;j<nSeqsPerOTU[i];j++){
2835 seqNumber.push_back(aaP[i][j]);
2837 for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
2838 aaP[i].push_back(0);
2844 for(int i=1;i<numOTUs;i++){
2845 cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
2848 seqIndex = seqNumber;
2855 catch(exception& e) {
2856 m->errorOut(e, "ShhherCommand", "getOTUData");
2860 /**************************************************************************************************/
2862 int ShhherCommand::calcCentroidsDriver(int numOTUs,
2863 vector<int>& cumNumSeqs,
2864 vector<int>& nSeqsPerOTU,
2865 vector<int>& seqIndex,
2866 vector<short>& change, //did the centroid sequence change? 0 = no; 1 = yes
2867 vector<int>& centroids, //the representative flowgram for each cluster m
2868 vector<double>& singleTau, //tMaster->adTau: 1-D Tau vector (1xnumSeqs)
2869 vector<int>& mapSeqToUnique,
2870 vector<short>& uniqueFlowgrams,
2871 vector<short>& flowDataIntI,
2872 vector<int>& lengths,
2874 vector<int>& seqNumber){
2876 //this function gets the most likely homopolymer length at a flow position for a group of sequences
2881 for(int i=0;i<numOTUs;i++){
2883 if (m->control_pressed) { break; }
2887 int minFlowGram = 100000000;
2888 double minFlowValue = 1e8;
2889 change[i] = 0; //FALSE
2891 for(int j=0;j<nSeqsPerOTU[i];j++){
2892 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
2895 if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
2896 vector<double> adF(nSeqsPerOTU[i]);
2897 vector<int> anL(nSeqsPerOTU[i]);
2899 for(int j=0;j<nSeqsPerOTU[i];j++){
2900 int index = cumNumSeqs[i] + j;
2901 int nI = seqIndex[index];
2902 int nIU = mapSeqToUnique[nI];
2905 for(k=0;k<position;k++){
2911 anL[position] = nIU;
2912 adF[position] = 0.0000;
2917 for(int j=0;j<nSeqsPerOTU[i];j++){
2918 int index = cumNumSeqs[i] + j;
2919 int nI = seqIndex[index];
2921 double tauValue = singleTau[seqNumber[index]];
2923 for(int k=0;k<position;k++){
2924 double dist = getDistToCentroid(anL[k], nI, lengths[nI], uniqueFlowgrams, flowDataIntI, numFlowCells);
2925 adF[k] += dist * tauValue;
2929 for(int j=0;j<position;j++){
2930 if(adF[j] < minFlowValue){
2932 minFlowValue = adF[j];
2936 if(centroids[i] != anL[minFlowGram]){
2938 centroids[i] = anL[minFlowGram];
2941 else if(centroids[i] != -1){
2949 catch(exception& e) {
2950 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
2954 /**************************************************************************************************/
2956 double ShhherCommand::getDistToCentroid(int cent, int flow, int length, vector<short>& uniqueFlowgrams,
2957 vector<short>& flowDataIntI, int numFlowCells){
2960 int flowAValue = cent * numFlowCells;
2961 int flowBValue = flow * numFlowCells;
2965 for(int i=0;i<length;i++){
2966 dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
2971 return dist / (double)length;
2973 catch(exception& e) {
2974 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
2978 /**************************************************************************************************/
2980 double ShhherCommand::getNewWeights(int numOTUs, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU, vector<double>& singleTau, vector<int>& seqNumber, vector<double>& weight){
2983 double maxChange = 0;
2985 for(int i=0;i<numOTUs;i++){
2987 if (m->control_pressed) { break; }
2989 double difference = weight[i];
2992 for(int j=0;j<nSeqsPerOTU[i];j++){
2993 int index = cumNumSeqs[i] + j;
2994 double tauValue = singleTau[seqNumber[index]];
2995 weight[i] += tauValue;
2998 difference = fabs(weight[i] - difference);
2999 if(difference > maxChange){ maxChange = difference; }
3003 catch(exception& e) {
3004 m->errorOut(e, "ShhherCommand", "getNewWeights");
3009 /**************************************************************************************************/
3011 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){
3015 vector<long double> P(numSeqs, 0);
3018 for(int i=0;i<numOTUs;i++){
3019 if(weight[i] > MIN_WEIGHT){
3025 for(int i=0;i<numOTUs;i++){
3027 if (m->control_pressed) { break; }
3029 for(int j=0;j<nSeqsPerOTU[i];j++){
3030 int index = cumNumSeqs[i] + j;
3031 int nI = seqIndex[index];
3032 double singleDist = dist[seqNumber[index]];
3034 P[nI] += weight[i] * exp(-singleDist * sigma);
3038 for(int i=0;i<numSeqs;i++){
3039 if(P[i] == 0){ P[i] = DBL_EPSILON; }
3044 nLL = nLL -(double)numSeqs * log(sigma);
3048 catch(exception& e) {
3049 m->errorOut(e, "ShhherCommand", "getNewWeights");
3054 /**************************************************************************************************/
3056 int ShhherCommand::checkCentroids(int numOTUs, vector<int>& centroids, vector<double>& weight){
3058 vector<int> unique(numOTUs, 1);
3060 for(int i=0;i<numOTUs;i++){
3061 if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
3066 for(int i=0;i<numOTUs;i++){
3068 if (m->control_pressed) { break; }
3071 for(int j=i+1;j<numOTUs;j++){
3074 if(centroids[j] == centroids[i]){
3078 weight[i] += weight[j];
3088 catch(exception& e) {
3089 m->errorOut(e, "ShhherCommand", "checkCentroids");
3093 /**************************************************************************************************/
3095 void ShhherCommand::calcNewDistances(int numSeqs, int numOTUs, vector<int>& nSeqsPerOTU, vector<double>& dist,
3096 vector<double>& weight, vector<short>& change, vector<int>& centroids,
3097 vector<vector<int> >& aaP, vector<double>& singleTau, vector<vector<int> >& aaI,
3098 vector<int>& seqNumber, vector<int>& seqIndex,
3099 vector<short>& uniqueFlowgrams,
3100 vector<short>& flowDataIntI, int numFlowCells, vector<int>& lengths){
3105 vector<double> newTau(numOTUs,0);
3106 vector<double> norms(numSeqs, 0);
3107 nSeqsPerOTU.assign(numOTUs, 0);
3109 for(int i=0;i<numSeqs;i++){
3111 if (m->control_pressed) { break; }
3113 int indexOffset = i * numOTUs;
3115 double offset = 1e8;
3117 for(int j=0;j<numOTUs;j++){
3119 if(weight[j] > MIN_WEIGHT && change[j] == 1){
3120 dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i], uniqueFlowgrams, flowDataIntI, numFlowCells);
3123 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
3124 offset = dist[indexOffset + j];
3128 for(int j=0;j<numOTUs;j++){
3129 if(weight[j] > MIN_WEIGHT){
3130 newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
3131 norms[i] += newTau[j];
3138 for(int j=0;j<numOTUs;j++){
3139 newTau[j] /= norms[i];
3142 for(int j=0;j<numOTUs;j++){
3143 if(newTau[j] > MIN_TAU){
3145 int oldTotal = total;
3149 singleTau.resize(total, 0);
3150 seqNumber.resize(total, 0);
3151 seqIndex.resize(total, 0);
3153 singleTau[oldTotal] = newTau[j];
3155 aaP[j][nSeqsPerOTU[j]] = oldTotal;
3156 aaI[j][nSeqsPerOTU[j]] = i;
3164 catch(exception& e) {
3165 m->errorOut(e, "ShhherCommand", "calcNewDistances");
3169 /**************************************************************************************************/
3171 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){
3174 for(int i=0;i<numOTUs;i++){
3176 if (m->control_pressed) { return 0; }
3178 cumNumSeqs[i] = index;
3179 for(int j=0;j<nSeqsPerOTU[i];j++){
3180 seqNumber[index] = aaP[i][j];
3181 seqIndex[index] = aaI[i][j];
3189 catch(exception& e) {
3190 m->errorOut(e, "ShhherCommand", "fill");
3194 /**************************************************************************************************/
3196 void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector<int>& seqNumber, vector<int>& seqIndex, vector<int>& cumNumSeqs, vector<int>& nSeqsPerOTU,
3197 vector<int>& otuData, vector<double>& singleTau, vector<double>& dist, vector<vector<int> >& aaP, vector<vector<int> >& aaI){
3200 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
3202 for(int i=0;i<numOTUs;i++){
3204 if (m->control_pressed) { break; }
3206 for(int j=0;j<nSeqsPerOTU[i];j++){
3207 int index = cumNumSeqs[i] + j;
3208 double tauValue = singleTau[seqNumber[index]];
3209 int sIndex = seqIndex[index];
3210 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
3214 for(int i=0;i<numSeqs;i++){
3215 double maxTau = -1.0000;
3217 for(int j=0;j<numOTUs;j++){
3218 if(bigTauMatrix[i * numOTUs + j] > maxTau){
3219 maxTau = bigTauMatrix[i * numOTUs + j];
3224 otuData[i] = maxOTU;
3227 nSeqsPerOTU.assign(numOTUs, 0);
3229 for(int i=0;i<numSeqs;i++){
3230 int index = otuData[i];
3232 singleTau[i] = 1.0000;
3235 aaP[index][nSeqsPerOTU[index]] = i;
3236 aaI[index][nSeqsPerOTU[index]] = i;
3238 nSeqsPerOTU[index]++;
3241 fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI);
3243 catch(exception& e) {
3244 m->errorOut(e, "ShhherCommand", "setOTUs");
3248 /**************************************************************************************************/
3250 void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string qualityFileName, vector<int> otuCounts, vector<int>& nSeqsPerOTU, vector<int>& seqNumber,
3251 vector<double>& singleTau, vector<short>& flowDataIntI, vector<short>& uniqueFlowgrams, vector<int>& cumNumSeqs,
3252 vector<int>& mapUniqueToSeq, vector<string>& seqNameVector, vector<int>& centroids, vector<vector<int> >& aaI){
3256 ofstream qualityFile;
3257 m->openOutputFile(qualityFileName, qualityFile);
3259 qualityFile.setf(ios::fixed, ios::floatfield);
3260 qualityFile.setf(ios::showpoint);
3261 qualityFile << setprecision(6);
3263 vector<vector<int> > qualities(numOTUs);
3264 vector<double> pr(HOMOPS, 0);
3267 for(int i=0;i<numOTUs;i++){
3269 if (m->control_pressed) { break; }
3274 if(nSeqsPerOTU[i] > 0){
3275 qualities[i].assign(1024, -1);
3277 while(index < numFlowCells){
3278 double maxPrValue = 1e8;
3279 short maxPrIndex = -1;
3280 double count = 0.0000;
3282 pr.assign(HOMOPS, 0);
3284 for(int j=0;j<nSeqsPerOTU[i];j++){
3285 int lIndex = cumNumSeqs[i] + j;
3286 double tauValue = singleTau[seqNumber[lIndex]];
3287 int sequenceIndex = aaI[i][j];
3288 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
3292 for(int s=0;s<HOMOPS;s++){
3293 pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
3297 maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
3298 maxPrValue = pr[maxPrIndex];
3300 if(count > MIN_COUNT){
3302 double norm = 0.0000;
3304 for(int s=0;s<HOMOPS;s++){
3305 norm += exp(-(pr[s] - maxPrValue));
3308 for(int s=1;s<=maxPrIndex;s++){
3310 double temp = 0.0000;
3312 U += exp(-(pr[s-1]-maxPrValue))/norm;
3320 temp = floor(-10 * temp);
3321 value = (int)floor(temp);
3322 if(value > 100){ value = 100; }
3324 qualities[i][base] = (int)value;
3334 if(otuCounts[i] > 0){
3335 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
3337 int j=4; //need to get past the first four bases
3338 while(qualities[i][j] != -1){
3339 qualityFile << qualities[i][j] << ' ';
3340 if (j > qualities[i].size()) { break; }
3343 qualityFile << endl;
3346 qualityFile.close();
3347 outputNames.push_back(qualityFileName); outputTypes["qfile"].push_back(qualityFileName);
3350 catch(exception& e) {
3351 m->errorOut(e, "ShhherCommand", "writeQualities");
3356 /**************************************************************************************************/
3358 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){
3362 m->openOutputFile(fastaFileName, fastaFile);
3364 vector<string> names(numOTUs, "");
3366 for(int i=0;i<numOTUs;i++){
3368 if (m->control_pressed) { break; }
3370 int index = centroids[i];
3372 if(otuCounts[i] > 0){
3373 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
3377 for(int j=0;j<numFlowCells;j++){
3379 char base = flowOrder[j % flowOrder.length()];
3380 for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
3385 if (newSeq.length() >= 4) { fastaFile << newSeq.substr(4) << endl; }
3386 else { fastaFile << "NNNN" << endl; }
3391 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
3393 if(thisCompositeFASTAFileName != ""){
3394 m->appendFiles(fastaFileName, thisCompositeFASTAFileName);
3397 catch(exception& e) {
3398 m->errorOut(e, "ShhherCommand", "writeSequences");
3403 /**************************************************************************************************/
3405 void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string nameFileName, vector<int> otuCounts, vector<string>& seqNameVector, vector<vector<int> >& aaI, vector<int>& nSeqsPerOTU){
3409 m->openOutputFile(nameFileName, nameFile);
3411 for(int i=0;i<numOTUs;i++){
3413 if (m->control_pressed) { break; }
3415 if(otuCounts[i] > 0){
3416 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
3418 for(int j=1;j<nSeqsPerOTU[i];j++){
3419 nameFile << ',' << seqNameVector[aaI[i][j]];
3426 outputNames.push_back(nameFileName); outputTypes["name"].push_back(nameFileName);
3429 if(thisCompositeNamesFileName != ""){
3430 m->appendFiles(nameFileName, thisCompositeNamesFileName);
3433 catch(exception& e) {
3434 m->errorOut(e, "ShhherCommand", "writeNames");
3439 /**************************************************************************************************/
3441 void ShhherCommand::writeGroups(string groupFileName, string fileRoot, int numSeqs, vector<string>& seqNameVector){
3444 m->openOutputFile(groupFileName, groupFile);
3446 for(int i=0;i<numSeqs;i++){
3447 if (m->control_pressed) { break; }
3448 groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
3451 outputNames.push_back(groupFileName); outputTypes["group"].push_back(groupFileName);
3454 catch(exception& e) {
3455 m->errorOut(e, "ShhherCommand", "writeGroups");
3460 /**************************************************************************************************/
3462 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){
3464 ofstream otuCountsFile;
3465 m->openOutputFile(otuCountsFileName, otuCountsFile);
3467 string bases = flowOrder;
3469 for(int i=0;i<numOTUs;i++){
3471 if (m->control_pressed) {
3474 //output the translated version of the centroid sequence for the otu
3475 if(otuCounts[i] > 0){
3476 int index = centroids[i];
3478 otuCountsFile << "ideal\t";
3479 for(int j=8;j<numFlowCells;j++){
3480 char base = bases[j % bases.length()];
3481 for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
3482 otuCountsFile << base;
3485 otuCountsFile << endl;
3487 for(int j=0;j<nSeqsPerOTU[i];j++){
3488 int sequence = aaI[i][j];
3489 otuCountsFile << seqNameVector[sequence] << '\t';
3493 for(int k=0;k<lengths[sequence];k++){
3494 char base = bases[k % bases.length()];
3495 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
3497 for(int s=0;s<freq;s++){
3499 //otuCountsFile << base;
3503 if (newSeq.length() >= 4) { otuCountsFile << newSeq.substr(4) << endl; }
3504 else { otuCountsFile << "NNNN" << endl; }
3506 otuCountsFile << endl;
3509 otuCountsFile.close();
3510 outputNames.push_back(otuCountsFileName); outputTypes["counts"].push_back(otuCountsFileName);
3513 catch(exception& e) {
3514 m->errorOut(e, "ShhherCommand", "writeClusters");
3519 /**************************************************************************************************/
3521 void ShhherCommand::getSingleLookUp(){
3523 // these are the -log probabilities that a signal corresponds to a particular homopolymer length
3524 singleLookUp.assign(HOMOPS * NUMBINS, 0);
3527 ifstream lookUpFile;
3528 m->openInputFile(lookupFileName, lookUpFile);
3530 for(int i=0;i<HOMOPS;i++){
3532 if (m->control_pressed) { break; }
3535 lookUpFile >> logFracFreq;
3537 for(int j=0;j<NUMBINS;j++) {
3538 lookUpFile >> singleLookUp[index];
3544 catch(exception& e) {
3545 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
3550 /**************************************************************************************************/
3552 void ShhherCommand::getJointLookUp(){
3555 // the most likely joint probability (-log) that two intenities have the same polymer length
3556 jointLookUp.resize(NUMBINS * NUMBINS, 0);
3558 for(int i=0;i<NUMBINS;i++){
3560 if (m->control_pressed) { break; }
3562 for(int j=0;j<NUMBINS;j++){
3564 double minSum = 100000000;
3566 for(int k=0;k<HOMOPS;k++){
3567 double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
3569 if(sum < minSum) { minSum = sum; }
3571 jointLookUp[i * NUMBINS + j] = minSum;
3575 catch(exception& e) {
3576 m->errorOut(e, "ShhherCommand", "getJointLookUp");
3581 /**************************************************************************************************/
3583 double ShhherCommand::getProbIntensity(int intIntensity){
3586 double minNegLogProb = 100000000;
3589 for(int i=0;i<HOMOPS;i++){//loop signal strength
3591 if (m->control_pressed) { break; }
3593 float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
3594 if(negLogProb < minNegLogProb) { minNegLogProb = negLogProb; }
3597 return minNegLogProb;
3599 catch(exception& e) {
3600 m->errorOut(e, "ShhherCommand", "getProbIntensity");