2 * removerarecommand.cpp
5 * Created by westcott on 1/21/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "removerarecommand.h"
11 #include "sequence.hpp"
13 #include "sharedutilities.h"
14 #include "inputdata.h"
16 //**********************************************************************************************************************
17 vector<string> RemoveRareCommand::setParameters(){
19 CommandParameter plist("list", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(plist);
20 CommandParameter prabund("rabund", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(prabund);
21 CommandParameter psabund("sabund", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(psabund);
22 CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pshared);
23 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup);
24 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
25 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
26 CommandParameter pnseqs("nseqs", "Number", "", "0", "", "", "",false,true); parameters.push_back(pnseqs);
27 CommandParameter pbygroup("bygroup", "Boolean", "", "f", "", "", "",false,true); parameters.push_back(pbygroup);
28 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
29 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
31 vector<string> myArray;
32 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
36 m->errorOut(e, "RemoveRareCommand", "setParameters");
40 //**********************************************************************************************************************
41 string RemoveRareCommand::getHelpString(){
43 string helpString = "";
44 helpString += "The remove.rare command parameters are list, rabund, sabund, shared, group, label, groups, bygroup and nseqs.\n";
45 helpString += "The remove.rare command reads one of the following file types: list, rabund, sabund or shared file. It outputs a new file after removing the rare otus.\n";
46 helpString += "The groups parameter allows you to specify which of the groups you would like analyzed. Default=all. You may separate group names with dashes.\n";
47 helpString += "The label parameter is used to analyze specific labels in your input. default=all. You may separate label names with dashes.\n";
48 helpString += "The bygroup parameter is only valid with the shared file. default=f, meaning remove any OTU that has nseqs or fewer sequences across all groups.\n";
49 helpString += "bygroups=T means remove any OTU that has nseqs or fewer sequences in each group (if groupA has 1 sequence and group B has 100 sequences in OTUZ and nseqs=1, then set the groupA count for OTUZ to 0 and keep groupB's count at 100.) \n";
50 helpString += "The nseqs parameter allows you to set the cutoff for an otu to be deemed rare. It is required.\n";
51 helpString += "The remove.rare command should be in the following format: remove.rare(shared=yourSharedFile, nseqs=yourRareCutoff).\n";
52 helpString += "Example remove.rare(shared=amazon.fn.shared, nseqs=2).\n";
53 helpString += "Note: No spaces between parameter labels (i.e. shared), '=' and parameters (i.e.yourSharedFile).\n";
57 m->errorOut(e, "RemoveRareCommand", "getHelpString");
61 //**********************************************************************************************************************
62 string RemoveRareCommand::getOutputFileNameTag(string type, string inputName=""){
64 string outputFileName = "";
65 map<string, vector<string> >::iterator it;
67 //is this a type this command creates
68 it = outputTypes.find(type);
69 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
71 if (type == "rabund") { outputFileName = "pick" + m->getExtension(inputName); }
72 else if (type == "sabund") { outputFileName = "pick" + m->getExtension(inputName); }
73 else if (type == "shared") { outputFileName = "pick" + m->getExtension(inputName); }
74 else if (type == "group") { outputFileName = "pick" + m->getExtension(inputName); }
75 else if (type == "list") { outputFileName = "pick" + m->getExtension(inputName); }
76 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
78 return outputFileName;
81 m->errorOut(e, "RemoveRareCommand", "getOutputFileNameTag");
86 //**********************************************************************************************************************
87 RemoveRareCommand::RemoveRareCommand(){
89 abort = true; calledHelp = true;
91 vector<string> tempOutNames;
92 outputTypes["rabund"] = tempOutNames;
93 outputTypes["sabund"] = tempOutNames;
94 outputTypes["list"] = tempOutNames;
95 outputTypes["group"] = tempOutNames;
96 outputTypes["shared"] = tempOutNames;
99 m->errorOut(e, "RemoveRareCommand", "RemoveRareCommand");
103 //**********************************************************************************************************************
104 RemoveRareCommand::RemoveRareCommand(string option) {
106 abort = false; calledHelp = false;
109 //allow user to run help
110 if(option == "help") { help(); abort = true; calledHelp = true; }
111 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
114 vector<string> myArray = setParameters();
116 OptionParser parser(option);
117 map<string,string> parameters = parser.getParameters();
119 ValidParameters validParameter;
120 map<string,string>::iterator it;
122 //check to make sure all parameters are valid for command
123 for (it = parameters.begin(); it != parameters.end(); it++) {
124 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
127 //initialize outputTypes
128 vector<string> tempOutNames;
129 outputTypes["rabund"] = tempOutNames;
130 outputTypes["sabund"] = tempOutNames;
131 outputTypes["list"] = tempOutNames;
132 outputTypes["group"] = tempOutNames;
133 outputTypes["shared"] = tempOutNames;
135 //if the user changes the output directory command factory will send this info to us in the output parameter
136 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
138 //if the user changes the input directory command factory will send this info to us in the output parameter
139 string inputDir = validParameter.validFile(parameters, "inputdir", false);
140 if (inputDir == "not found"){ inputDir = ""; }
143 it = parameters.find("list");
144 //user has given a template file
145 if(it != parameters.end()){
146 path = m->hasPath(it->second);
147 //if the user has not given a path then, add inputdir. else leave path alone.
148 if (path == "") { parameters["list"] = inputDir + it->second; }
151 it = parameters.find("group");
152 //user has given a template file
153 if(it != parameters.end()){
154 path = m->hasPath(it->second);
155 //if the user has not given a path then, add inputdir. else leave path alone.
156 if (path == "") { parameters["group"] = inputDir + it->second; }
159 it = parameters.find("sabund");
160 //user has given a template file
161 if(it != parameters.end()){
162 path = m->hasPath(it->second);
163 //if the user has not given a path then, add inputdir. else leave path alone.
164 if (path == "") { parameters["sabund"] = inputDir + it->second; }
167 it = parameters.find("rabund");
168 //user has given a template file
169 if(it != parameters.end()){
170 path = m->hasPath(it->second);
171 //if the user has not given a path then, add inputdir. else leave path alone.
172 if (path == "") { parameters["rabund"] = inputDir + it->second; }
175 it = parameters.find("shared");
176 //user has given a template file
177 if(it != parameters.end()){
178 path = m->hasPath(it->second);
179 //if the user has not given a path then, add inputdir. else leave path alone.
180 if (path == "") { parameters["shared"] = inputDir + it->second; }
185 //check for file parameters
186 listfile = validParameter.validFile(parameters, "list", true);
187 if (listfile == "not open") { abort = true; }
188 else if (listfile == "not found") { listfile = ""; }
189 else { m->setListFile(listfile); }
191 sabundfile = validParameter.validFile(parameters, "sabund", true);
192 if (sabundfile == "not open") { abort = true; }
193 else if (sabundfile == "not found") { sabundfile = ""; }
194 else { m->setSabundFile(sabundfile); }
196 rabundfile = validParameter.validFile(parameters, "rabund", true);
197 if (rabundfile == "not open") { abort = true; }
198 else if (rabundfile == "not found") { rabundfile = ""; }
199 else { m->setRabundFile(rabundfile); }
201 groupfile = validParameter.validFile(parameters, "group", true);
202 if (groupfile == "not open") { groupfile = ""; abort = true; }
203 else if (groupfile == "not found") { groupfile = ""; }
204 else { m->setGroupFile(groupfile); }
206 sharedfile = validParameter.validFile(parameters, "shared", true);
207 if (sharedfile == "not open") { sharedfile = ""; abort = true; }
208 else if (sharedfile == "not found") { sharedfile = ""; }
209 else { m->setSharedFile(sharedfile); }
211 if ((sharedfile == "") && (listfile == "") && (rabundfile == "") && (sabundfile == "")) {
212 //is there are current file available for any of these?
213 //give priority to shared, then list, then rabund, then sabund
214 //if there is a current shared file, use it
215 sharedfile = m->getSharedFile();
216 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
218 listfile = m->getListFile();
219 if (listfile != "") { m->mothurOut("Using " + listfile + " as input file for the list parameter."); m->mothurOutEndLine(); }
221 rabundfile = m->getRabundFile();
222 if (rabundfile != "") { m->mothurOut("Using " + rabundfile + " as input file for the rabund parameter."); m->mothurOutEndLine(); }
224 sabundfile = m->getSabundFile();
225 if (sabundfile != "") { m->mothurOut("Using " + sabundfile + " as input file for the sabund parameter."); m->mothurOutEndLine(); }
227 m->mothurOut("No valid current files. You must provide a list, sabund, rabund or shared file."); m->mothurOutEndLine();
235 groups = validParameter.validFile(parameters, "groups", false);
236 if (groups == "not found") { groups = "all"; }
237 m->splitAtDash(groups, Groups);
239 label = validParameter.validFile(parameters, "label", false);
240 if (label == "not found") { label = ""; }
242 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
243 else { allLines = 1; }
246 string temp = validParameter.validFile(parameters, "nseqs", false);
247 if (temp == "not found") { m->mothurOut("nseqs is a required parameter."); m->mothurOutEndLine(); abort = true; }
248 else { m->mothurConvert(temp, nseqs); }
250 temp = validParameter.validFile(parameters, "bygroup", false); if (temp == "not found") { temp = "f"; }
251 byGroup = m->isTrue(temp);
253 if (byGroup && (sharedfile == "")) { m->mothurOut("The byGroup parameter is only valid with a shared file."); m->mothurOutEndLine(); }
255 if ((groupfile != "") && (listfile == "")) { m->mothurOut("A groupfile is only valid with a list file."); m->mothurOutEndLine(); groupfile = ""; }
259 catch(exception& e) {
260 m->errorOut(e, "RemoveRareCommand", "RemoveRareCommand");
264 //**********************************************************************************************************************
266 int RemoveRareCommand::execute(){
269 if (abort == true) { if (calledHelp) { return 0; } return 2; }
271 if (m->control_pressed) { return 0; }
273 //read through the correct file and output lines you want to keep
274 if (sabundfile != "") { processSabund(); }
275 if (rabundfile != "") { processRabund(); }
276 if (listfile != "") { processList(); }
277 if (sharedfile != "") { processShared(); }
279 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
281 if (outputNames.size() != 0) {
282 m->mothurOutEndLine();
283 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
284 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
285 m->mothurOutEndLine();
287 //set rabund file as new current rabundfile
289 itTypes = outputTypes.find("rabund");
290 if (itTypes != outputTypes.end()) {
291 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRabundFile(current); }
294 itTypes = outputTypes.find("sabund");
295 if (itTypes != outputTypes.end()) {
296 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSabundFile(current); }
299 itTypes = outputTypes.find("group");
300 if (itTypes != outputTypes.end()) {
301 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
304 itTypes = outputTypes.find("list");
305 if (itTypes != outputTypes.end()) {
306 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
309 itTypes = outputTypes.find("shared");
310 if (itTypes != outputTypes.end()) {
311 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSharedFile(current); }
318 catch(exception& e) {
319 m->errorOut(e, "RemoveRareCommand", "execute");
324 //**********************************************************************************************************************
325 int RemoveRareCommand::processList(){
327 string thisOutputDir = outputDir;
328 if (outputDir == "") { thisOutputDir += m->hasPath(listfile); }
329 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(listfile)) + getOutputFileNameTag("list", listfile);
330 string outputGroupFileName = thisOutputDir + m->getRootName(m->getSimpleName(groupfile)) + getOutputFileNameTag("group", groupfile);
331 ofstream out, outGroup;
332 m->openOutputFile(outputFileName, out);
334 bool wroteSomething = false;
336 //you must provide a label because the names in the listfile need to be consistent
337 string thisLabel = "";
338 if (allLines) { m->mothurOut("For the listfile you must select one label, using first label in your listfile."); m->mothurOutEndLine(); }
339 else if (labels.size() > 1) { m->mothurOut("For the listfile you must select one label, using " + (*labels.begin()) + "."); m->mothurOutEndLine(); thisLabel = *labels.begin(); }
340 else { thisLabel = *labels.begin(); }
342 InputData input(listfile, "list");
343 ListVector* list = input.getListVector();
345 //get first one or the one we want
346 if (thisLabel != "") {
347 //use smart distancing
348 set<string> userLabels; userLabels.insert(thisLabel);
349 set<string> processedLabels;
350 string lastLabel = list->getLabel();
351 while((list != NULL) && (userLabels.size() != 0)) {
352 if(userLabels.count(list->getLabel()) == 1){
353 processedLabels.insert(list->getLabel());
354 userLabels.erase(list->getLabel());
358 if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
359 processedLabels.insert(list->getLabel());
360 userLabels.erase(list->getLabel());
362 list = input.getListVector(lastLabel);
365 lastLabel = list->getLabel();
367 list = input.getListVector();
369 if (userLabels.size() != 0) {
370 m->mothurOut("Your file does not include the label " + thisLabel + ". I will use " + lastLabel + "."); m->mothurOutEndLine();
371 list = input.getListVector(lastLabel);
375 //if groupfile is given then use it
377 if (groupfile != "") {
378 groupMap = new GroupMap(groupfile); groupMap->readMap();
380 vector<string> namesGroups = groupMap->getNamesOfGroups();
381 util.setGroups(Groups, namesGroups);
382 m->openOutputFile(outputGroupFileName, outGroup);
387 //make a new list vector
389 newList.setLabel(list->getLabel());
392 for (int i = 0; i < list->getNumBins(); i++) {
393 if (m->control_pressed) { if (groupfile != "") { delete groupMap; outGroup.close(); m->mothurRemove(outputGroupFileName); } out.close(); m->mothurRemove(outputFileName); return 0; }
395 //parse out names that are in accnos file
396 string binnames = list->get(i);
397 vector<string> names;
398 string saveBinNames = binnames;
399 m->splitAtComma(binnames, names);
401 vector<string> newGroupFile;
402 if (groupfile != "") {
403 vector<string> newNames;
405 for(int k = 0; k < names.size(); k++) {
406 string group = groupMap->getGroup(names[k]);
408 if (m->inUsersGroups(group, Groups)) {
409 newGroupFile.push_back(names[k] + "\t" + group);
411 newNames.push_back(names[k]);
412 saveBinNames += names[k] + ",";
416 saveBinNames = saveBinNames.substr(0, saveBinNames.length()-1);
419 if (names.size() > nseqs) { //keep bin
420 newList.push_back(saveBinNames);
421 for(int k = 0; k < newGroupFile.size(); k++) { outGroup << newGroupFile[k] << endl; }
425 //print new listvector
426 if (newList.getNumBins() != 0) {
427 wroteSomething = true;
433 if (groupfile != "") { outGroup.close(); outputTypes["group"].push_back(outputGroupFileName); outputNames.push_back(outputGroupFileName); }
435 if (wroteSomething == false) { m->mothurOut("Your file contains only rare sequences."); m->mothurOutEndLine(); }
436 outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName);
440 catch(exception& e) {
441 m->errorOut(e, "RemoveRareCommand", "processList");
445 //**********************************************************************************************************************
446 int RemoveRareCommand::processSabund(){
448 string thisOutputDir = outputDir;
449 if (outputDir == "") { thisOutputDir += m->hasPath(sabundfile); }
450 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(sabundfile)) + getOutputFileNameTag("sabund", sabundfile);
451 outputTypes["sabund"].push_back(outputFileName); outputNames.push_back(outputFileName);
454 m->openOutputFile(outputFileName, out);
456 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
457 InputData input(sabundfile, "sabund");
458 SAbundVector* sabund = input.getSAbundVector();
459 string lastLabel = sabund->getLabel();
460 set<string> processedLabels;
461 set<string> userLabels = labels;
463 while((sabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
465 if (m->control_pressed) { delete sabund; out.close(); return 0; }
467 if(allLines == 1 || labels.count(sabund->getLabel()) == 1){
469 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
470 processedLabels.insert(sabund->getLabel());
471 userLabels.erase(sabund->getLabel());
473 if (sabund->getMaxRank() > nseqs) {
474 for(int i = 1; i <=nseqs; i++) { sabund->set(i, 0); }
475 }else { sabund->clear(); }
477 if (sabund->getNumBins() > 0) { sabund->print(out); }
480 if ((m->anyLabelsToProcess(sabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
481 string saveLabel = sabund->getLabel();
484 sabund = input.getSAbundVector(lastLabel);
486 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
487 processedLabels.insert(sabund->getLabel());
488 userLabels.erase(sabund->getLabel());
490 if (sabund->getMaxRank() > nseqs) {
491 for(int i = 1; i <=nseqs; i++) { sabund->set(i, 0); }
492 }else { sabund->clear(); }
494 if (sabund->getNumBins() > 0) { sabund->print(out); }
496 //restore real lastlabel to save below
497 sabund->setLabel(saveLabel);
500 lastLabel = sabund->getLabel();
503 sabund = input.getSAbundVector();
506 if (m->control_pressed) { out.close(); return 0; }
508 //output error messages about any remaining user labels
509 set<string>::iterator it;
510 bool needToRun = false;
511 for (it = userLabels.begin(); it != userLabels.end(); it++) {
512 m->mothurOut("Your file does not include the label " + *it);
513 if (processedLabels.count(lastLabel) != 1) {
514 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
517 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
521 //run last label if you need to
522 if (needToRun == true) {
523 if (sabund != NULL) { delete sabund; }
524 sabund = input.getSAbundVector(lastLabel);
526 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
528 if (sabund->getMaxRank() > nseqs) {
529 for(int i = 1; i <=nseqs; i++) { sabund->set(i, 0); }
530 }else { sabund->clear(); }
532 if (sabund->getNumBins() > 0) { sabund->print(out); }
539 catch(exception& e) {
540 m->errorOut(e, "RemoveRareCommand", "processSabund");
544 //**********************************************************************************************************************
545 int RemoveRareCommand::processRabund(){
547 string thisOutputDir = outputDir;
548 if (outputDir == "") { thisOutputDir += m->hasPath(rabundfile); }
549 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(rabundfile)) + getOutputFileNameTag("rabund", rabundfile);
550 outputTypes["rabund"].push_back(outputFileName); outputNames.push_back(outputFileName);
553 m->openOutputFile(outputFileName, out);
555 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
556 InputData input(rabundfile, "rabund");
557 RAbundVector* rabund = input.getRAbundVector();
558 string lastLabel = rabund->getLabel();
559 set<string> processedLabels;
560 set<string> userLabels = labels;
562 while((rabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
564 if (m->control_pressed) { delete rabund; out.close(); return 0; }
566 if(allLines == 1 || labels.count(rabund->getLabel()) == 1){
568 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
569 processedLabels.insert(rabund->getLabel());
570 userLabels.erase(rabund->getLabel());
572 RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
573 for (int i = 0; i < rabund->getNumBins(); i++) {
574 if (rabund->get(i) > nseqs) {
575 newRabund.push_back(rabund->get(i));
578 if (newRabund.getNumBins() > 0) { newRabund.print(out); }
581 if ((m->anyLabelsToProcess(rabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
582 string saveLabel = rabund->getLabel();
585 rabund = input.getRAbundVector(lastLabel);
587 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
588 processedLabels.insert(rabund->getLabel());
589 userLabels.erase(rabund->getLabel());
591 RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
592 for (int i = 0; i < rabund->getNumBins(); i++) {
593 if (rabund->get(i) > nseqs) {
594 newRabund.push_back(rabund->get(i));
597 if (newRabund.getNumBins() > 0) { newRabund.print(out); }
599 //restore real lastlabel to save below
600 rabund->setLabel(saveLabel);
603 lastLabel = rabund->getLabel();
606 rabund = input.getRAbundVector();
609 if (m->control_pressed) { out.close(); return 0; }
611 //output error messages about any remaining user labels
612 set<string>::iterator it;
613 bool needToRun = false;
614 for (it = userLabels.begin(); it != userLabels.end(); it++) {
615 m->mothurOut("Your file does not include the label " + *it);
616 if (processedLabels.count(lastLabel) != 1) {
617 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
620 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
624 //run last label if you need to
625 if (needToRun == true) {
626 if (rabund != NULL) { delete rabund; }
627 rabund = input.getRAbundVector(lastLabel);
629 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
631 RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
632 for (int i = 0; i < rabund->getNumBins(); i++) {
633 if (rabund->get(i) > nseqs) {
634 newRabund.push_back(rabund->get(i));
637 if (newRabund.getNumBins() > 0) { newRabund.print(out); }
644 catch(exception& e) {
645 m->errorOut(e, "RemoveRareCommand", "processRabund");
649 //**********************************************************************************************************************
650 int RemoveRareCommand::processShared(){
652 m->setGroups(Groups);
654 string thisOutputDir = outputDir;
655 if (outputDir == "") { thisOutputDir += m->hasPath(sharedfile); }
656 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(sharedfile)) + getOutputFileNameTag("shared", sharedfile);
657 outputTypes["shared"].push_back(outputFileName); outputNames.push_back(outputFileName);
660 m->openOutputFile(outputFileName, out);
662 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
663 InputData input(sharedfile, "sharedfile");
664 vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
665 string lastLabel = lookup[0]->getLabel();
666 set<string> processedLabels;
667 set<string> userLabels = labels;
669 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
671 if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } out.close(); return 0; }
673 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
675 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
676 processedLabels.insert(lookup[0]->getLabel());
677 userLabels.erase(lookup[0]->getLabel());
679 if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
680 processLookup(lookup, out);
683 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
684 string saveLabel = lookup[0]->getLabel();
686 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
687 lookup = input.getSharedRAbundVectors(lastLabel);
689 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
690 processedLabels.insert(lookup[0]->getLabel());
691 userLabels.erase(lookup[0]->getLabel());
693 if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
694 processLookup(lookup, out);
696 //restore real lastlabel to save below
697 lookup[0]->setLabel(saveLabel);
700 lastLabel = lookup[0]->getLabel();
702 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
703 lookup = input.getSharedRAbundVectors();
706 if (m->control_pressed) { out.close(); return 0; }
708 //output error messages about any remaining user labels
709 set<string>::iterator it;
710 bool needToRun = false;
711 for (it = userLabels.begin(); it != userLabels.end(); it++) {
712 m->mothurOut("Your file does not include the label " + *it);
713 if (processedLabels.count(lastLabel) != 1) {
714 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
717 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
721 //run last label if you need to
722 if (needToRun == true) {
723 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
724 lookup = input.getSharedRAbundVectors(lastLabel);
726 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
728 if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
729 processLookup(lookup, out);
731 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
736 catch(exception& e) {
737 m->errorOut(e, "RemoveRareCommand", "processSabund");
741 //**********************************************************************************************************************
742 int RemoveRareCommand::processLookup(vector<SharedRAbundVector*>& lookup, ofstream& out){
745 vector<SharedRAbundVector> newRabunds; newRabunds.resize(lookup.size());
746 for (int i = 0; i < lookup.size(); i++) {
747 newRabunds[i].setGroup(lookup[i]->getGroup());
748 newRabunds[i].setLabel(lookup[i]->getLabel());
754 for (int i = 0; i < lookup[0]->getNumBins(); i++) {
757 if (m->control_pressed) { return 0; }
760 for (int j = 0; j < lookup.size(); j++) {
763 if (lookup[j]->getAbundance(i) > nseqs) {
764 newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
767 newRabunds[j].push_back(0, newRabunds[j].getGroup());
771 //eliminates zero otus
772 if (allZero) { for (int j = 0; j < newRabunds.size(); j++) { newRabunds[j].pop_back(); } }
776 for (int i = 0; i < lookup[0]->getNumBins(); i++) {
778 if (m->control_pressed) { return 0; }
781 //get total otu abundance
782 for (int j = 0; j < lookup.size(); j++) {
783 newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
784 totalAbund += lookup[j]->getAbundance(i);
787 //eliminates otus below rare cutoff
788 if (totalAbund <= nseqs) { for (int j = 0; j < newRabunds.size(); j++) { newRabunds[j].pop_back(); } }
792 //do we have any otus above the rare cutoff
793 if (newRabunds[0].getNumBins() != 0) {
794 for (int j = 0; j < newRabunds.size(); j++) {
795 out << newRabunds[j].getLabel() << '\t' << newRabunds[j].getGroup() << '\t';
796 newRabunds[j].print(out);
802 catch(exception& e) {
803 m->errorOut(e, "RemoveRareCommand", "processLookup");
807 //**********************************************************************************************************************