return data;
}
catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "Standard Error: " << e.what() << " has occurred in the BergerParker class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
catch(...) {
- cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "An unknown error has occurred in the BergerParker class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
}
int nb = 0;
int ns = 0;
- for(int i = rank->size()-1; i > 0; i--)
- {
+ for(int i = rank->size()-1; i > 0; i--) {
int cur = rank->get(i);
if(mr == 1 && cur > 0)
mr = i;
double sumObs = 0;
double maxDiff = 0;
- for(int i = 0; i < rdata.size(); i++)
- {
+ for(int i = 0; i < rdata.size(); i++) {
sumObs += rdata.get(i);
sumExp += numInd/numSpec*invSum(i+1,numSpec);
double diff = fabs(sumObs-sumExp);
maxDiff = diff;
}
+
data[0] = maxDiff/numInd;
data[1] = 0.886/sqrt(rdata.size());
data[2] = 1.031/sqrt(rdata.size());
/*cout << critVal << "\n";
cout << "If D-Statistic is less than the critical value then the data fits the Broken Stick model w/ 95% confidence.\n\n";*/
+
if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
return data;
}
catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "Standard Error: " << e.what() << " has occurred in the BStick class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
catch(...) {
- cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "An unknown error has occurred in the BStick class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
}
//Found on http://www.vgtu.lt/leidiniai/elektroniniai/Probability.pdf/Table%203.pdf
//Confidence Level .90 .95 .975 .99 .995 .999 .9995
- double values[33][7] = {{3.078, 6.314, 12.706, 31.821, 63.656, 318.289, 636.578},
+ double values[30][7] = {{3.078, 6.314, 12.706, 31.821, 63.656, 318.289, 636.578},
{1.886, 2.920, 4.303, 6.965, 9.925, 22.328, 31.600},
{1.638, 2.353, 3.182, 4.541, 5.841, 10.214, 12.924},
{1.533, 2.132, 2.776, 3.747, 4.604, 7.173, 8.610},
{1.314, 1.703, 2.052, 2.473, 2.771, 3.421, 3.689},
{1.313, 1.701, 2.048, 2.467, 2.763, 3.408, 3.674},
{1.311, 1.699, 2.045, 2.462, 2.756, 3.396, 3.660},
- {1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646},
- {1.296, 1.671, 2.000, 2.390, 2.660, 3.232, 3.460},
- {1.289, 1.658, 1.980, 2.358, 2.617, 3.160, 3.373},
- {1.282, 1.645, 1.960, 2.326, 2.576, 3.091, 3.291}};
+ {1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646}};
+
return values[row][col];
}
catch(exception& e) {
//calculate at 0 and the given increment
if((i == 0) || (i+1) % increment == 0){
+
//how many comparisons to make i.e. for group a, b, c = ab, ac, bc.
+
int n = 1;
for (int k = 0; k < (lookup.size() - 1); k++) { // pass cdd each set of groups to commpare
for (int l = n; l < lookup.size(); l++) {
#include "logsd.h"
#include "bergerparker.h"
#include "bstick.h"
+#include "goodscoverage.h"
+
+
#include "coverage.h"
+
//**********************************************************************************************************************
CollectCommand::CollectCommand(){
try {
}else if (globaldata->Estimators[i] == "bergerparker") {
cDisplays.push_back(new CollectDisplay(new BergerParker(), new OneColumnFile(fileNameRoot+"bergerparker")));
}else if (globaldata->Estimators[i] == "bstick") {
- cDisplays.push_back(new CollectDisplay(new BStick(), new OneColumnFile(fileNameRoot+"bstick")));
+ cDisplays.push_back(new CollectDisplay(new BStick(), new ThreeColumnFile(fileNameRoot+"bstick")));
+ }else if (globaldata->Estimators[i] == "goodscoverage") {
+ cDisplays.push_back(new CollectDisplay(new GoodsCoverage(), new OneColumnFile(fileNameRoot+"goodscoverage")));
}
}
}
#include "sharedlennon.h"
#include "sharedmorisitahorn.h"
#include "sharedbraycurtis.h"
+#include "sharedjackknife.h"
+#include "sharedwhittaker.h"
cDisplays.push_back(new CollectDisplay(new Whittaker(), new SharedOneColumnFile(fileNameRoot+"whittaker")));
}else if (globaldata->Estimators[i] == "sharednseqs") {
cDisplays.push_back(new CollectDisplay(new SharedNSeqs(), new SharedOneColumnFile(fileNameRoot+"shared.nseqs")));
+
}else if (globaldata->Estimators[i] == "ochiai") {
cDisplays.push_back(new CollectDisplay(new Ochiai(), new SharedOneColumnFile(fileNameRoot+"ochiai")));
}else if (globaldata->Estimators[i] == "anderberg") {
#include "readdistcommand.h"
#include "readtreecommand.h"
#include "readotucommand.h"
+#include "readseqscommand.h"
#include "clustercommand.h"
#include "parselistcommand.h"
#include "collectcommand.h"
#include "unifracweightedcommand.h"
#include "libshuffcommand.h"
#include "heatmapcommand.h"
+#include "filterseqscommand.h"
+#include "mothur.h"
#include "venncommand.h"
#include "nocommands.h"
#include "binsequencecommand.h"
if(commandName == "read.dist") { command = new ReadDistCommand(); }
else if(commandName == "read.otu") { command = new ReadOtuCommand(); }
+ else if(commandName == "read.seqs") { command = new ReadSeqsCommand(); }
else if(commandName == "read.tree") { command = new ReadTreeCommand(); }
else if(commandName == "cluster") { command = new ClusterCommand(); }
else if(commandName == "deconvolute") { command = new DeconvoluteCommand(); }
else if(commandName == "get.line") { command = new GetlineCommand(); }
else if(commandName == "libshuff") { command = new LibShuffCommand(); }
else if(commandName == "heatmap") { command = new HeatMapCommand(); }
+ else if(commandName == "filter.seqs") { command = new FilterSeqsCommand(); }
else if(commandName == "venn") { command = new VennCommand(); }
else if(commandName == "bin.seqs") { command = new BinSeqCommand(); }
else if(commandName == "get.oturep") { command = new GetOTURepCommand(); }
/******************************************************/
void ErrorCheck::refresh() {
+
//columnfile = globaldata->getColumnFile();
//phylipfile = globaldata->getPhylipFile();
//listfile = globaldata->getListFile();
//method = globaldata->getMethod();
//randomtree = globaldata->getRandomTree();
//sharedfile = globaldata->getSharedFile();
+
}
/*******************************************************/
//is it a valid command
if (validCommand->isValidCommand(commandName) != true) { return false; }
-
string parameter, value;
//reads in parameters and values
if (parameter == "name" ) { namefile = value; }
if (parameter == "order" ) { orderfile = value; }
if (parameter == "fasta" ) { fastafile = value; }
+ if (parameter == "nexus" ) { nexusfile = value; }
+ if (parameter == "clustal" ) { clustalfile = value; }
if (parameter == "tree" ) { treefile = value; }
- if (parameter == "group" ) { groupfile = value; }
- if (parameter == "shared" ) { sharedfile = value; }
+ if (parameter == "group" ) { groupfile = value; }
+ if (parameter == "shared" ) { sharedfile = value; }
if (parameter == "cutoff" ) { cutoff = value; }
if (parameter == "precision" ) { precision = value; }
if (parameter == "iters" ) { iters = value; }
if (parameter == "line" ) { line = value; }
if (parameter == "label" ) { label = value; }
if (parameter == "abund" ) { abund = value; }
- if (parameter == "random" ) { randomtree = value; }
- if (parameter == "sorted" ) { sorted = value; }
+ if (parameter == "random" ) { randomtree = value; }
+ if (parameter == "sorted" ) { sorted = value; }
+ if (parameter == "trump" ) { trump = value; }
+ if (parameter == "soft" ) { soft = value; }
+ if (parameter == "filter" ) { filter = value; }
if (parameter == "scale" ) { scale = value; }
-
+
}
//gets the last parameter and value
if (parameter == "group" ) { groupfile = value; }
if (parameter == "shared" ) { sharedfile = value; }
if (parameter == "fasta" ) { fastafile = value; }
+ if (parameter == "nexus" ) { nexusfile = value; }
+ if (parameter == "clustal" ) { clustalfile = value; }
if (parameter == "tree" ) { treefile = value; }
if (parameter == "cutoff" ) { cutoff = value; }
if (parameter == "precision" ) { precision = value; }
if (parameter == "random" ) { randomtree = value; }
if (parameter == "abund" ) { abund = value; }
if (parameter == "sorted" ) { sorted = value; }
+ if (parameter == "trump" ) { trump = value; }
+ if (parameter == "soft" ) { soft = value; }
+ if (parameter == "filter" ) { filter = value; }
if (parameter == "scale" ) { scale = value; }
+
}
}
}
}else if (commandName == "read.tree") {
validateTreeFiles(); //checks the treefile and groupfile parameters
+ }else if (commandName == "read.seqs") {
+ if ((fastafile == "") && (nexusfile == "") && (clustalfile == "") && (phylipfile == "")) { cout << "You must enter a fastafile, nexusfile, or clustalfile with the read.seqs() command." << endl; return false; }
}else if (commandName == "deconvolute") {
if (fastafile == "") { cout << "You must enter a fastafile with the deconvolute() command." << endl; return false; }
validateReadFiles();
}
}
+ if (commandName == "filter.seqs"){
+ if ((globaldata->getFastaFile() == "") && (globaldata->getNexusFile() == "") && (globaldata->getClustalFile() == "") && (globaldata->getPhylipFile() == "")) {
+ cout << "You must read either a fasta, nexus, clustal, or phylip file before you can use the filter.seqs command." << endl; return false;
+ }
+ }
+
if ((commandName == "bin.seqs")) {
if ((globaldata->getListFile() == "")) { cout << "You must read a list file before you can use the bin.seqs command." << endl; return false; }
validateBinFiles();
groupfile = "";
orderfile = "";
sharedfile = "";
+ fastafile = "";
+ nexusfile = "";
+ clustalfile = "";
line = "";
label = "";
method = "furthest";
void validateBinFiles();
void clear();
void refresh();
- string phylipfile, columnfile, listfile, rabundfile, sabundfile, namefile, groupfile, orderfile, fastafile, treefile, sharedfile, cutoff, format;
- string precision, method, fileroot, label, line, iters, jumble, freq, single, rarefaction, shared, summary, randomtree, abund, sorted, scale;
+ string phylipfile, columnfile, listfile, rabundfile, sabundfile, namefile, groupfile, orderfile, fastafile, nexusfile, clustalfile, treefile, sharedfile, cutoff, format;
+ string precision, method, fileroot, label, line, iters, jumble, freq, single, rarefaction, shared, summary, randomtree, abund, sorted, trump, soft, filter, scale;
+
string commandName, optionText;
bool errorFree;
--- /dev/null
+/*
+ * filterseqscommand.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 5/4/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "filterseqscommand.h"
+#include <iostream>
+#include <fstream>
+
+/**************************************************************************************/
+void FilterSeqsCommand::doTrump() {
+ //trump = globaldata->getTrump();
+//
+// for(int i = 0; i < db->size(); i++) {
+// Sequence cur = db->get(i);
+// string curAligned = cur.getAligned();
+//
+// for(int j = 0; j < curAligned.length-1; j++) {
+// string curChar = curAligned.substr(j, j+1);
+//
+// if(curChar.compare(trump) == 0)
+// columnsToRemove[j] = true;
+// }
+// }
+}
+
+/**************************************************************************************/
+void FilterSeqsCommand::doSoft() {
+ //soft = atoi(globaldata->getSoft().c_str());
+// vector<vector<int> > columnSymbolSums;
+// vector<vector<string> > columnSymbols;
+// for(int i = 0; i < db->get(0).getLength(); i++) {
+// vector<string> symbols;
+// vector<int> sums;
+// columnSymbols[i] = symbols;
+// columnSymbolSums[i] = sums;
+// }
+//
+// for(int i = 0; i < db->size(); i++) {
+// Sequence cur = db->get(i);
+// string curAligned = cur.getAligned();
+//
+// for(int j = 0; j < curAligned.length-1; j++) {
+// string curChar = curAligned.substr(j, j+1);
+// vector<string> curColumnSymbols = columnSymbols[j];
+//
+// bool newSymbol = true;
+//
+// for(int k = 0; j < curColumnSymbols.size(); j++)
+// if(curChar.compare(curColumnSymbols[k]) == 0) {
+// newSymbol = false;
+// columnSymbolSums[j][k]++;
+// }
+//
+// if(newSymbol) {
+// columnSymbols.push_back(curChar);
+// columnSymbolSums[j].push_back(1);
+// }
+// }
+// }
+//
+// for(int i = 0; i < columnSymbolSums.size(); i++) {
+// int totalSum = 0;
+// int max = 0;
+// vector<int> curColumn = columnSymbolSums[i];
+//
+// for(int j = 0; j < curColumn.size(); j++) {
+// int curSum = curColumn[j];
+// if(curSum > max)
+// max = curSum;
+// totalSum += curSum;
+// }
+//
+// if((double)max/(double)totalSum * 100 < soft)
+// columnsToRemove[i] = true;
+// }
+}
+void FilterSeqsCommand::doFilter() {}
+/**************************************************************************************/
+int FilterSeqsCommand::execute() {
+ try {
+ globaldata = GlobalData::getInstance();
+ filename = globaldata->inputFileName;
+
+ if(globaldata->getFastaFile().compare("") != 0) {
+ readFasta = new ReadFasta(filename);
+ readFasta->read();
+ db = readFasta->getDB();
+ }
+
+ else if(globaldata->getNexusFile().compare("") != 0) {
+ readNexus = new ReadNexus(filename);
+ readNexus->read();
+ db = readNexus->getDB();
+ }
+
+ else if(globaldata->getClustalFile().compare("") != 0) {
+ readClustal = new ReadClustal(filename);
+ readClustal->read();
+ db = readClustal->getDB();
+ }
+
+ else if(globaldata->getPhylipFile().compare("") != 0) {
+ readPhylip = new ReadPhylip(filename);
+ readPhylip->read();
+ db = readPhylip->getDB();
+ }
+
+ for(int i = 0; i < db->get(0).getLength(); i++)
+ columnsToRemove[i] = false;
+
+ // Trump
+ if(globaldata->getTrump().compare("") != 0) {
+
+
+ }
+
+ // Soft
+ if(globaldata->getSoft().compare("") != 0) {}
+
+
+
+
+ // Filter
+ //if(globaldata->getFilter().compare("") != 0) {
+//
+// filter = globaldata->getFilter();
+// ifstream filehandle;
+// openInputFile(filter, filehandle);
+//
+// char c;
+// int count = 0;
+// while(!filehandle.eof()) {
+// c = filehandle.get();
+// if(c == '0')
+// columnsToRemove[count] = true;
+// count++;
+// }
+// }
+
+
+
+
+
+ return 0;
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the DeconvoluteCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the DeconvoluteCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/**************************************************************************************/
--- /dev/null
+#ifndef FILTERSEQSCOMMAND_H
+#define FILTERSEQSCOMMAND_H
+
+/*
+ * filterseqscommand.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 5/4/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "command.hpp"
+#include "mothur.h"
+#include "globaldata.hpp"
+#include "readfasta.h"
+#include "readnexus.h"
+#include "readclustal.h"
+#include "readseqsphylip.h"
+
+using namespace std;
+
+class FilterSeqsCommand : public Command {
+
+public:
+ FilterSeqsCommand() {};
+ ~FilterSeqsCommand() {};
+ int execute();
+
+private:
+ void doTrump();
+ void doSoft();
+ void doFilter();
+
+ GlobalData* globaldata;
+ string filename, trump, filter;
+
+ ReadFasta* readFasta;
+ ReadNexus* readNexus;
+ ReadClustal* readClustal;
+ ReadPhylip* readPhylip;
+
+ vector<bool> columnsToRemove;
+ SequenceDB* db;
+ double soft;
+};
+
+#endif
\ No newline at end of file
int nb = 0;
int ns = 0;
- for(int i = rank->size()-1; i > 0; i--)
- {
+ for(int i = rank->size()-1; i > 0; i--) {
int cur = rank->get(i);
if(mr == 1 && cur > 0)
mr = i;
double k = .5;
double step = .49999;
- while(fabs(min - numInd*kEq(k, (double)numSpec)) > .0001) //This uses a binary search to find the value of k.
- {
+ while(fabs(min - numInd*kEq(k, (double)numSpec)) > .0001) { //This uses a binary search to find the value of k.
if(numInd*kEq(k, numSpec) > min)
k += step;
else
/*cout << critVal << "\n";
cout << "If D-Statistic is less than the critical value then the data fits the Geometric Series model w/ 95% confidence.\n\n";*/
-
+
if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
return data;
}
catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "Standard Error: " << e.what() << " has occurred in the Geom class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
catch(...) {
- cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "An unknown error has occurred in the Geom class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
}
string label;
int numBins = 0;
int count = -1;
- while(in.good())
- {
+ while(in.good()) {
if(count > numBins)
count = 0;
- if(count == 0)
- {
+ if(count == 0) {
cout << label << "\n";
in >> numBins;
}
int numBins = 0;
int count = -1;
int line = 1;
- while(in.good())
- {
+ while(in.good()) {
if(count > numBins)
count = 0;
- if(count == 0)
- {
+ if(count == 0) {
cout << line << "\n";
in >> numBins;
line++;
if (key == "list" ) { listfile = value; inputFileName = value; fileroot = value; format = "list"; }
if (key == "rabund" ) { rabundfile = value; inputFileName = value; fileroot = value; format = "rabund"; }
if (key == "sabund" ) { sabundfile = value; inputFileName = value; fileroot = value; format = "sabund"; }
- if (key == "fasta" ) { fastafile = value; inputFileName = value; fileroot = value; format = "fasta"; }
+ if (key == "fasta" ) { fastafile = value; inputFileName = value; fileroot = value; format = "fasta"; }
+ if (key == "nexus" ) { nexusfile = value; inputFileName = value; fileroot = value; format = "nexus"; }
+ if (key == "clustal" ) { clustalfile = value; inputFileName = value; fileroot = value; format = "clustal"; }
if (key == "tree" ) { treefile = value; inputFileName = value; fileroot = value; format = "tree"; }
if (key == "shared" ) { sharedfile = value; inputFileName = value; fileroot = value; format = "sharedfile"; }
if (key == "name" ) { namefile = value; }
if (key == "step") { step = value; }
if (key == "form") { form = value; }
if (key == "sorted") { sorted = value; }
+ if (key == "vertical") { vertical = value; }
+ if (key == "trump") { trump = value; }
+ if (key == "filter") { filter = value; }
+ if (key == "soft") { soft = value; }
if (key == "scale") { scale = value; }
+
+
if (key == "line") {//stores lines to be used in a set
lines.clear();
labels.clear();
if (key == "rabund" ) { rabundfile = value; inputFileName = value; fileroot = value; format = "rabund"; }
if (key == "sabund" ) { sabundfile = value; inputFileName = value; fileroot = value; format = "sabund"; }
if (key == "fasta" ) { fastafile = value; inputFileName = value; fileroot = value; format = "fasta"; }
+ if (key == "nexus" ) { nexusfile = value; inputFileName = value; fileroot = value; format = "nexus"; }
+ if (key == "clustal" ) { clustalfile = value; inputFileName = value; fileroot = value; format = "clustal"; }
if (key == "tree" ) { treefile = value; inputFileName = value; fileroot = value; format = "tree"; }
if (key == "shared" ) { sharedfile = value; inputFileName = value; fileroot = value; format = "sharedfile"; }
if (key == "name" ) { namefile = value; }
if (key == "step") { step = value; }
if (key == "form") { form = value; }
if (key == "sorted") { sorted = value; }
+ if (key == "vertical") { vertical = value; }
+ if (key == "trump") { trump = value; }
+ if (key == "filter") { filter = value; }
+ if (key == "soft") { soft = value; }
if (key == "scale") { scale = value; }
+
+
+
if (key == "line") {//stores lines to be used in a vector
lines.clear();
//input defaults for calculators
if (commandName == "collect.single") {
+
if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-shannon-npshannon-simpson"; }
Estimators.clear();
splitAtDash(calc, Estimators);
splitAtDash(calc, Estimators);
}
if (commandName == "collect.shared") {
+
if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan"; }
Estimators.clear();
splitAtDash(calc, Estimators);
string GlobalData::getTreeFile() { return treefile; }
string GlobalData::getSharedFile() { return sharedfile; }
string GlobalData::getFastaFile() { return fastafile; }
+string GlobalData::getNexusFile() { return nexusfile; }
+string GlobalData::getClustalFile() { return clustalfile; }
string GlobalData::getCutOff() { return cutoff; }
string GlobalData::getFormat() { return format; }
string GlobalData::getPrecision() { return precision; }
string GlobalData::getStep() { return step; }
string GlobalData::getForm() { return form; }
string GlobalData::getSorted() { return sorted; }
+string GlobalData::getTrump() { return trump; }
+string GlobalData::getSoft() { return soft; }
+string GlobalData::getFilter() { return filter; }
string GlobalData::getScale() { return scale; }
+
void GlobalData::setListFile(string file) { listfile = file; inputFileName = file;}
void GlobalData::setRabundFile(string file) { rabundfile = file; inputFileName = file;}
void GlobalData::setSabundFile(string file) { sabundfile = file; inputFileName = file;}
groupfile = "";
orderfile = "";
fastafile = "";
+ nexusfile = "";
+ clustalfile = "";
treefile = "";
sharedfile = "";
cutoff = "10.00";
step = "0.01";
form = "integral";
sorted = "T"; //F means don't sort, T means sort.
- scale = "log10";
+ vertical = "";
+ trump = "";
+ filter = "";
+ soft = "";
+ scale = "log10";
+
}
//*******************************************************/
GroupMap* gGroupmap;
FullMatrix* gMatrix;
TreeMap* gTreemap;
- string inputFileName, helpRequest, commandName;
+ string inputFileName, helpRequest, commandName, vertical;
bool allLines;
vector<string> Estimators, Groups; //holds estimators to be used
set<int> lines; //hold lines to be used
string getGroupFile();
string getOrderFile();
string getFastaFile();
+ string getNexusFile();
+ string getClustalFile();
string getTreeFile();
string getSharedFile();
string getCutOff();
string getStep();
string getForm();
string getSorted();
+
+ string getTrump();
+ string getSoft();
+ string getFilter();
+
+
string getScale();
+
void setListFile(string);
void setPhylipFile(string);
void setColumnFile(string);
private:
- string phylipfile, columnfile, listfile, rabundfile, sabundfile, namefile, groupfile, orderfile, fastafile, treefile, sharedfile, line, label, randomtree, groups;
- string cutoff, format, precision, method, fileroot, iters, jumble, freq, calc, abund, step, form, sorted, scale;
+
+ string phylipfile, columnfile, listfile, rabundfile, sabundfile, namefile, groupfile, orderfile, fastafile, nexusfile, clustalfile, treefile, sharedfile, line, label, randomtree, groups;
+ string cutoff, format, precision, method, fileroot, iters, jumble, freq, calc, abund, step, form, sorted, trump, soft, filter, scale;
+
static GlobalData* _uniqueInstance;
GlobalData( const GlobalData& ); // Disable copy constructor
--- /dev/null
+/*
+ * goodscoverage.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/8/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "goodscoverage.h"
+#include "calculator.h"
+
+/**********************************************************/
+
+EstOutput GoodsCoverage::getValues(SAbundVector* rank){
+ try {
+ data.resize(1,0);
+
+ double numSingletons = rank->get(1);
+ double totalIndividuals = rank->getNumSeqs();
+
+ data[0] = 1 - numSingletons/totalIndividuals;
+
+ if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
+
+ return data;
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the GoodsCoverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the GoodsCoverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+/***********************************************************************/
+
--- /dev/null
+#ifndef GOODSCOVERAGE_H
+#define GOODSCOVERAGE_H
+
+/*
+ * goodscoverage.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/8/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+#include "calculator.h"
+
+/* This class implements the LogSD estimator on single group.
+It is a child of the calculator class. */
+
+/***********************************************************************/
+
+class GoodsCoverage : public Calculator {
+
+public:
+ GoodsCoverage() : Calculator("goodscoverage", 1, false) {};
+ EstOutput getValues(SAbundVector*);
+ EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
+
+private:
+};
+
+/***********************************************************************/
+
+#endif
double x = .5;
double step = .4999999999;
- while(fabs(snRatio - logS(x)) > .00001) //This uses a binary search to find the value of x.
- {
+ while(fabs(snRatio - logS(x)) > .00001) { //This uses a binary search to find the value of x.
if(logS(x) > snRatio)
x += step;
else
double octSumExp = 0;
double sumExp = 0;
double maxDiff = 0;
- for(int y = 1; y < rank->size(); y++)
- {
- if(y - .5 < pow(2.0, oct))
- {
+ for(int y = 1; y < rank->size(); y++) {
+ if(y - .5 < pow(2.0, oct)) {
octSumObs += rank->get(y);
octSumExp += alpha*pow(x,y)/(y);
}
- else
- {
+ else {
sumObs += octSumObs;
octSumObs = rank->get(y);
oct++;
}
- if(y == rank->size()-1)
- {
+ if(y == rank->size()-1) {
sumObs += octSumObs;
sumExp += octSumExp;
}
cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
cout << "If D Test Statistic is greater than the critical value then the data fits the Log Series Distribution model w/ 95% confidence.\n\n";*/
+
data[0] = (maxDiff + .5)/numSpec;
data[1] = 0.886/sqrt(numSpec);
data[2] = 1.031/sqrt(numSpec);
return data;
}
catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "Standard Error: " << e.what() << " has occurred in the LogSD class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
catch(...) {
- cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "An unknown error has occurred in the LogSD class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
}
int r3Ind = 0;
int sumSpec = 0;
int iqSum = 0;
- for(int i = 1; i < rank->size(); i++)
- {
+ for(int i = 1; i < rank->size(); i++) {
if(r1 != -1 && r3 != -1)
i = rank->size();
sumSpec += rank->get(i);
- if(r1 == -1 && sumSpec >= numSpec*.25)
- {
+ if(r1 == -1 && sumSpec >= numSpec*.25) {
r1 = rank->get(i);
r1Ind = i;
}
- else if(r3 == -1 && sumSpec >= numSpec*.75)
- {
+ else if(r3 == -1 && sumSpec >= numSpec*.75) {
r3 = rank->get(i);
r3Ind = i;
}
return data;
}
catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "Standard Error: " << e.what() << " has occurred in the QStat class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
catch(...) {
- cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "An unknown error has occurred in the QStat class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
}
--- /dev/null
+/*
+ * readclustal.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/24/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "readclustal.h"
+#include <iostream>
+#include <fstream>
+
+/*******************************************************************************/
+ReadClustal::ReadClustal(string file) {
+ try {
+ openInputFile(file, filehandle);
+ clustalFile = file;
+ globaldata = GlobalData::getInstance();
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the ReadTree class function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/*******************************************************************************/
+ReadClustal::~ReadClustal(){
+// for(int i = 0; i < sequencedb.getNumSeqs(); i++)
+// delete sequencedb.get(i);
+}
+/*******************************************************************************/
+void ReadClustal::read() {
+ string temp;
+ string name;
+ string sequence;
+ string firstName = "";
+ for(int i = 0; i < 6; i++)
+ filehandle >> temp;
+
+ int count = 0;
+ int numSeqs = 0;
+ bool firstDone = false;
+
+ while(!filehandle.eof()) {
+ filehandle >> name;
+ if(numSeqs != 0) {
+ if(count == numSeqs)
+ count = 0;
+ }
+ else if(!firstDone && firstName.compare("") == 0)
+ firstName = name;
+ else if(!firstDone && firstName.compare(name) == 0) {
+ numSeqs = count;
+ firstDone = true;
+ count = 0;
+ }
+
+ if(name.find_first_of("*") == -1) {
+ filehandle >> sequence;
+ if(!firstDone) {
+ Sequence newSeq(name, sequence);
+ sequencedb.add(newSeq);
+ }
+ else
+ sequencedb.set(count, sequencedb.get(count).getAligned() + sequence);
+
+ count++;
+ }
+ }
+ filehandle.close();
+}
+
+/*********************************************************************************/
+SequenceDB* ReadClustal::getDB() {
+ return &sequencedb;
+}
--- /dev/null
+#ifndef READCLUSTAL_H
+#define READCLUSTAL_H
+
+/*
+ * readclustal.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/24/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+using namespace std;
+
+#include "globaldata.hpp"
+#include "sequencedb.h"
+#include "mothur.h"
+
+/**********************************************************************************/
+
+class ReadClustal {
+
+ public:
+ ReadClustal(string);
+ ~ReadClustal();
+ void read();
+ SequenceDB* getDB();
+
+ private:
+ GlobalData* globaldata;
+ string clustalFile;
+ ifstream filehandle;
+ SequenceDB sequencedb;
+ int readOk; // readOk = 0 means success, readOk = 1 means error(s).
+
+
+};
+
+#endif
\ No newline at end of file
#include "readdistcommand.h"
#include "readphylip.h"
#include "readcolumn.h"
+#include "readmatrix.hpp"
ReadDistCommand::ReadDistCommand(){
try {
--- /dev/null
+/*
+ * readfasta.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/21/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "readfasta.h"
+#include <iostream>
+#include <fstream>
+
+/*******************************************************************************/
+ReadFasta::ReadFasta(string file) {
+ try {
+ openInputFile(file, filehandle);
+ fastaFile = file;
+ globaldata = GlobalData::getInstance();
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the ReadTree class function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/*******************************************************************************/
+ReadFasta::~ReadFasta(){
+ //for(int i = 0; i < sequencedb.getNumSeqs(); i++)
+ //delete sequencedb.get(i);
+}
+/*******************************************************************************/
+void ReadFasta::read() {
+ string name = "";
+ string sequence = "";
+ string temp;
+ int count = 0;
+ while(!filehandle.eof()){
+ if(count == 0)
+ filehandle >> temp;
+ if(temp.substr(0,1).compare(">") == 0) {
+ if(count != 0) {
+ Sequence newSequence(name, sequence);
+ sequencedb.add(newSequence);
+ sequence = "";
+ }
+ else
+ count++;
+ name = temp.substr(1,temp.length()-1);
+ }
+ else
+ sequence += temp;
+
+ filehandle >> temp;
+ }
+ Sequence newSequence(name, sequence);
+ sequencedb.add(newSequence);
+
+ filehandle.close();
+}
+
+/*********************************************************************************/
+SequenceDB* ReadFasta::getDB() {
+ return &sequencedb;
+}
--- /dev/null
+#ifndef READFASTA_H
+#define READFASTA_H
+
+/*
+ * readfasta.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/21/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+using namespace std;
+
+#include "globaldata.hpp"
+#include "sequencedb.h"
+#include "mothur.h"
+
+/**********************************************************************************/
+
+class ReadFasta {
+
+ public:
+ ReadFasta(string);
+ ~ReadFasta();
+ void read();
+ SequenceDB* getDB();
+
+ private:
+ GlobalData* globaldata;
+ string fastaFile;
+ ifstream filehandle;
+ SequenceDB sequencedb;
+ int readOk; // readOk = 0 means success, readOk = 1 means error(s).
+
+
+};
+
+#endif
\ No newline at end of file
--- /dev/null
+/*
+ * readnexusaln.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/22/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "readnexus.h"
+#include <iostream>
+#include <fstream>
+
+/*******************************************************************************/
+ReadNexus::ReadNexus(string file) {
+ try {
+ openInputFile(file, filehandle);
+ nexusFile = file;
+ globaldata = GlobalData::getInstance();
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the ReadTree class function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/*******************************************************************************/
+ReadNexus::~ReadNexus(){
+// for(int i = 0; i < sequencedb.getNumSeqs(); i++)
+// delete sequencedb.get(i);
+}
+/*******************************************************************************/
+void ReadNexus::read() {
+ string temp;
+ string name;
+ string sequence;
+ for(int i = 0; i < 5; i++)
+ filehandle >> temp;
+
+ int numSeqs = atoi(temp.substr(temp.find_first_of("=")+1, temp.length() - temp.find_first_of("=") - 1).c_str());
+
+ for(int i = 0; i < 9; i++)
+ filehandle >> temp;
+
+ int count = 0;
+ bool firstDone = false;
+ while(!filehandle.eof()){
+ filehandle >> name;
+ filehandle >> sequence;
+ if(name.compare(";") != 0) {
+ if(!firstDone) {
+ Sequence newSeq(name, sequence);
+ sequencedb.add(newSeq);
+ }
+ else
+ sequencedb.set(count, sequencedb.get(count).getAligned() + sequence);
+
+ count++;
+ if(count == numSeqs) {
+ if(!firstDone)
+ firstDone = true;
+ count = 0;
+ }
+ }
+ }
+ filehandle.close();
+}
+
+/*********************************************************************************/
+SequenceDB* ReadNexus::getDB() {
+ return &sequencedb;
+}
--- /dev/null
+#ifndef READNEXUS_H
+#define READNEXUS_H
+
+/*
+ * readnuxus.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/22/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+using namespace std;
+
+#include "globaldata.hpp"
+#include "sequencedb.h"
+#include "mothur.h"
+
+/**********************************************************************************/
+
+class ReadNexus {
+
+ public:
+ ReadNexus(string);
+ ~ReadNexus();
+ void read();
+ SequenceDB* getDB();
+
+ private:
+ GlobalData* globaldata;
+ string nexusFile;
+ ifstream filehandle;
+ SequenceDB sequencedb;
+ int readOk; // readOk = 0 means success, readOk = 1 means error(s).
+
+
+};
+
+#endif
\ No newline at end of file
--- /dev/null
+#ifndef READNEXUSALN_H
+#define READNEXUSALN_H
+
+/*
+ * readnexusaln.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/22/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+using namespace std;
+
+#include "globaldata.hpp"
+#include "sequencedb.h"
+#include "utilities.hpp"
+
+/**********************************************************************************/
+
+class ReadNexus {
+
+ public:
+ ReadNexus(string);
+ ~ReadNexus();
+ void read();
+
+ private:
+ GlobalData* globaldata;
+ string nexusFile;
+ ifstream filehandle;
+ SequenceDB sequencedb;
+ int readOk; // readOk = 0 means success, readOk = 1 means error(s).
+
+
+};
+
+#endif
\ No newline at end of file
/***********************************************************************/
ReadPhylipMatrix::ReadPhylipMatrix(string distFile){
-
- successOpen = openInputFile(distFile, fileHandle);
-
+
+ successOpen = openInputFile(distFile, fileHandle);
+
}
/***********************************************************************/
void ReadPhylipMatrix::read(NameAssignment* nameMap){
- try {
-
- float distance;
- int square, nseqs;
- string name;
- vector<string> matrixNames;
-
- fileHandle >> nseqs >> name;
+ try {
+
+ float distance;
+ int square, nseqs;
+ string name;
+ vector<string> matrixNames;
+
+ fileHandle >> nseqs >> name;
- matrixNames.push_back(name);
+ matrixNames.push_back(name);
- if(nameMap == NULL){
- list = new ListVector(nseqs);
- list->set(0, name);
- }
- else{
- list = new ListVector(nameMap->getListVector());
- if(nameMap->count(name)==0){ cout << "Error: Sequence '" << name << "' was not found in the names file, please correct" << endl; }
- }
-
- char d;
- while((d=fileHandle.get()) != EOF){
-
- if(isalnum(d)){
- square = 1;
- fileHandle.putback(d);
- for(int i=0;i<nseqs;i++){
- fileHandle >> distance;
- }
- break;
- }
- if(d == '\n'){
- square = 0;
- break;
- }
- }
-
- Progress* reading;
-
- if(square == 0){
+ if(nameMap == NULL){
+ list = new ListVector(nseqs);
+ list->set(0, name);
+ }
+ else{
+ list = new ListVector(nameMap->getListVector());
+ if(nameMap->count(name)==0){ cout << "Error: Sequence '" << name << "' was not found in the names file, please correct" << endl; }
+ }
+
+ char d;
+ while((d=fileHandle.get()) != EOF){
+
+ if(isalnum(d)){
+ square = 1;
+ fileHandle.putback(d);
+ for(int i=0;i<nseqs;i++){
+ fileHandle >> distance;
+ }
+ break;
+ }
+ if(d == '\n'){
+ square = 0;
+ break;
+ }
+ }
+
+ Progress* reading;
+
+ if(square == 0){
- reading = new Progress("Reading matrix: ", nseqs * (nseqs - 1) / 2);
-
- int index = 0;
-
- for(int i=1;i<nseqs;i++){
- fileHandle >> name;
- matrixNames.push_back(name);
-
- //there's A LOT of repeated code throughout this method...
- if(nameMap == NULL){
- list->set(i, name);
-
- for(int j=0;j<i;j++){
- fileHandle >> distance;
-
- if (distance == -1) { distance = 1000000; }
-
- if(distance < cutoff){
- PCell value(i, j, distance);
- D->addCell(value);
- }
- index++;
- reading->update(index);
- }
-
- }
- else{
- if(nameMap->count(name)==0){ cout << "Error: Sequence '" << name << "' was not found in the names file, please correct" << endl; }
-
- for(int j=0;j<i;j++){
- fileHandle >> distance;
-
- if (distance == -1) { distance = 1000000; }
-
- if(distance < cutoff){
- PCell value(nameMap->get(matrixNames[i]), nameMap->get(matrixNames[j]), distance);
- D->addCell(value);
- }
- index++;
- reading->update(index);
- }
- }
- }
- }
- else{
+ reading = new Progress("Reading matrix: ", nseqs * (nseqs - 1) / 2);
+
+ int index = 0;
+
+ for(int i=1;i<nseqs;i++){
+ fileHandle >> name;
+ matrixNames.push_back(name);
+
+ //there's A LOT of repeated code throughout this method...
+ if(nameMap == NULL){
+ list->set(i, name);
+
+ for(int j=0;j<i;j++){
+ fileHandle >> distance;
+
+ if (distance == -1) { distance = 1000000; }
+
+ if(distance < cutoff){
+ PCell value(i, j, distance);
+ D->addCell(value);
+ }
+ index++;
+ reading->update(index);
+ }
+
+ }
+ else{
+ if(nameMap->count(name)==0){ cout << "Error: Sequence '" << name << "' was not found in the names file, please correct" << endl; }
+
+ for(int j=0;j<i;j++){
+ fileHandle >> distance;
+
+ if (distance == -1) { distance = 1000000; }
+
+ if(distance < cutoff){
+ PCell value(nameMap->get(matrixNames[i]), nameMap->get(matrixNames[j]), distance);
+ D->addCell(value);
+ }
+ index++;
+ reading->update(index);
+ }
+ }
+ }
+ }
+ else{
- reading = new Progress("Reading matrix: ", nseqs * nseqs);
-
- int index = nseqs;
-
- for(int i=1;i<nseqs;i++){
- fileHandle >> name;
- matrixNames.push_back(name);
-
- if(nameMap == NULL){
- list->set(i, name);
- for(int j=0;j<nseqs;j++){
- fileHandle >> distance;
-
- if (distance == -1) { distance = 1000000; }
-
- if(distance < cutoff && j < i){
- PCell value(i, j, distance);
- D->addCell(value);
- }
- index++;
- reading->update(index);
- }
-
- }
- else{
- if(nameMap->count(name)==0){ cout << "Error: Sequence '" << name << "' was not found in the names file, please correct" << endl; }
-
- for(int j=0;j<nseqs;j++){
- fileHandle >> distance;
-
- if (distance == -1) { distance = 1000000; }
-
- if(distance < cutoff && j < i){
- PCell value(nameMap->get(matrixNames[i]), nameMap->get(matrixNames[j]), distance);
- D->addCell(value);
- }
- index++;
- reading->update(index);
- }
- }
- }
- }
- reading->finish();
- delete reading;
+ reading = new Progress("Reading matrix: ", nseqs * nseqs);
+
+ int index = nseqs;
+
+ for(int i=1;i<nseqs;i++){
+ fileHandle >> name;
+ matrixNames.push_back(name);
+
+ if(nameMap == NULL){
+ list->set(i, name);
+ for(int j=0;j<nseqs;j++){
+ fileHandle >> distance;
+
+ if (distance == -1) { distance = 1000000; }
+
+ if(distance < cutoff && j < i){
+ PCell value(i, j, distance);
+ D->addCell(value);
+ }
+ index++;
+ reading->update(index);
+ }
+
+ }
+ else{
+ if(nameMap->count(name)==0){ cout << "Error: Sequence '" << name << "' was not found in the names file, please correct" << endl; }
+
+ for(int j=0;j<nseqs;j++){
+ fileHandle >> distance;
+
+ if (distance == -1) { distance = 1000000; }
+
+ if(distance < cutoff && j < i){
+ PCell value(nameMap->get(matrixNames[i]), nameMap->get(matrixNames[j]), distance);
+ D->addCell(value);
+ }
+ index++;
+ reading->update(index);
+ }
+ }
+ }
+ }
+ reading->finish();
+ delete reading;
- list->setLabel("0");
- fileHandle.close();
+ list->setLabel("0");
+ fileHandle.close();
- if(nameMap != NULL){
- for(int i=0;i<matrixNames.size();i++){
- nameMap->erase(matrixNames[i]);
- }
- if(nameMap->size() > 0){
- //should probably tell them what is missing if we missed something
- cout << "missed something" << '\t' << nameMap->size() << endl;
- }
- }
+ if(nameMap != NULL){
+ for(int i=0;i<matrixNames.size();i++){
+ nameMap->erase(matrixNames[i]);
+ }
+ if(nameMap->size() > 0){
+ //should probably tell them what is missing if we missed something
+ cout << "missed something" << '\t' << nameMap->size() << endl;
+ }
+ }
- }
- catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the ReadPhylipMatrix class Function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
- exit(1);
- }
- catch(...) {
- cout << "An unknown error has occurred in the ReadPhylipMatrix class function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
- exit(1);
- }
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the ReadPhylipMatrix class Function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the ReadPhylipMatrix class function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
}
/***********************************************************************/
ReadPhylipMatrix::~ReadPhylipMatrix(){
- delete D;
- delete list;
+ delete D;
+ delete list;
}
-
--- /dev/null
+/*
+ * readseqscommand.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/13/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "readseqscommand.h"
+
+//**********************************************************************************************************************
+ReadSeqsCommand::ReadSeqsCommand(){
+ try {
+ globaldata = GlobalData::getInstance();
+
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the ReadOtuCommand class Function ReadOtuCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the ReadOtuCommand class function ReadOtuCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+//**********************************************************************************************************************
+
+ReadSeqsCommand::~ReadSeqsCommand(){
+ //delete readFasta->getDB();
+// delete readNexus->getDB();
+// delete readClustal->getDB();
+// delete readPhylip->getDB();
+}
+
+//**********************************************************************************************************************
+
+int ReadSeqsCommand::execute(){
+ try {
+ filebuf fb;
+
+ //fb.open ("fasta.txt",ios::out);
+// readFasta->read();
+// SequenceDB* db = readFasta->getDB();
+
+ //fb.open("nexus.txt",ios::out);
+// readNexus->read();
+// SequenceDB* db = readNexus->getDB();
+
+ //fb.open("clustal.txt",ios::out);
+// readClustal->read();
+// SequenceDB* db = readClustal->getDB();
+
+ //fb.open("phylip.txt",ios::out);
+// readPhylip->read();
+// SequenceDB* db = readPhylip->getDB();
+
+
+
+ //for(int i = 0; i < db->size(); i++) {
+// cout << db->get(i).getLength() << "\n" << db->get(i).getName() << ": " << db->get(i).getUnaligned() << "\n\n";
+// }
+
+ //ostream os(&fb);
+// db->print(os);
+// fb.close();
+
+ return 0;
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the ReadOtuCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the ReadOtuCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
--- /dev/null
+#ifndef READSEQSCOMMAND_H
+#define READSEQSCOMMAND_H
+/*
+ * readseqscommand.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/13/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "command.hpp"
+#include "readfasta.h"
+#include "readnexus.h"
+#include "readclustal.h"
+#include "readseqsphylip.h"
+#include "inputdata.h"
+#include "groupmap.h"
+#include "sharedcommand.h"
+#include "parselistcommand.h"
+
+/* The read.otu must be run before you execute a collect.single, rarefaction.single, summary.single,
+collect.shared, rarefaction.shared or summary.shared command. Mothur will generate a .list, .rabund and .sabund
+upon completion of the cluster command or you may use your own. The read.otu command parameter options are
+listfile, rabundfile, sabundfile, groupfile and orderfile. The reaad.otu command can be used in two ways.
+The first is to read a listfile, rabundfile or sabundfile and run the collect.single, rarefaction.single or summary.single.
+For this use the read.otu command should be in the following format: read.otu(listfile=yourListFile, orderfile=yourOrderFile).
+The listfile, rabundfile or sabundfile parameter is required, but you may only use one of them.
+The second way to use the read.otu command is to read a listfile and a groupfile so you can use the collect.shared,
+rarefaction.shared or summary.shared commands. In this case the read.otu command should be in the following format:
+read.otu(listfile=yourListFile, groupfile=yourGroupFile). The listfile parameter and groupfile paramaters are required.
+When using the command the second way read.otu command parses the .list file and separates it into groups.
+It outputs a .shared file containing the OTU information for each group. The read.otu command also outputs a .list file for each group. */
+
+class GlobalData;
+
+class ReadSeqsCommand : public Command {
+public:
+ ReadSeqsCommand();
+ ~ReadSeqsCommand();
+ int execute();
+
+private:
+ GlobalData* globaldata;
+ ReadFasta* readFasta;
+ ReadNexus* readNexus;
+ ReadClustal* readClustal;
+ ReadPhylip* readPhylip;
+ InputData* input;
+ Command* shared;
+ Command* parselist;
+ string filename;
+};
+
+#endif
--- /dev/null
+/*
+ * readphylip.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/24/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "readseqsphylip.h"
+#include <iostream>
+#include <fstream>
+
+/*******************************************************************************/
+bool ReadPhylip::isSeq(string seq) {
+ string validChars[] = {"A","G","C","T","U","N","-"};
+
+ for(int i = 0; i < seq.length(); i++) {
+ bool valid = false;
+ string c = seq.substr(i,1);
+ for(int k = 0; k < 7; k++)
+ if(c.compare(validChars[k]) == 0) {
+ valid = true;
+ k = 7;
+ }
+ if(!valid)
+ return false;
+ }
+
+ return true;
+}
+
+/*******************************************************************************/
+ReadPhylip::ReadPhylip(string file) {
+ try {
+ openInputFile(file, filehandle);
+ phylipFile = file;
+ globaldata = GlobalData::getInstance();
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the ReadTree class function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/*******************************************************************************/
+ReadPhylip::~ReadPhylip(){
+// for(int i = 0; i < sequencedb.getNumSeqs(); i++)
+// delete sequencedb.get(i);
+}
+/*******************************************************************************/
+void ReadPhylip::read() {
+ string temp;
+ string name;
+ string sequence;
+
+ int count = 0;
+ int letterCount = 0;
+ int numCols = 0;
+ filehandle >> temp;
+ int numSeqs = atoi(temp.c_str());
+ filehandle >> temp;
+ int numLetters = atoi(temp.c_str());
+
+ bool firstDone = false;
+ bool last = false;
+ filehandle >> name;
+
+ while(!filehandle.eof()) {
+ if(!firstDone) {
+ sequence = "";
+ if(count == 0) {
+ filehandle >> temp;
+ while(isSeq(temp)) {
+ sequence += temp;
+ numCols++;
+ filehandle >> temp;
+ }
+ letterCount += sequence.length();
+ }
+ else {
+ for(int i = 0; i < numCols; i++) {
+ filehandle >> temp;
+ sequence += temp;
+ }
+ if(count < numSeqs-1)
+ filehandle >> temp;
+ }
+ Sequence newSeq(name, sequence);
+ sequencedb.add(newSeq);
+ if(count < numSeqs-1)
+ name = temp;
+ }
+ else {
+ sequence = "";
+ for(int i = 0; i < numCols; i++) {
+ filehandle >> temp;
+ sequence += temp;
+ if(count == 0)
+ letterCount += temp.length();
+ if(letterCount == numLetters && count == 0) {
+ numCols = i + 1;
+ i = numCols;
+ }
+ }
+ if(!(last && count == 0))
+ sequencedb.set(count, sequencedb.get(count).getAligned() + sequence);
+ if(letterCount == numLetters && count == 0)
+ last = true;
+ }
+
+ count++;
+
+ if(count == numSeqs) {
+ firstDone = true;
+ count = 0;
+ }
+ }
+ filehandle.close();
+}
+
+/*********************************************************************************/
+SequenceDB* ReadPhylip::getDB() {
+ return &sequencedb;
+}
--- /dev/null
+#ifndef READPHYLIP_H
+#define READPHYLIP_H
+
+/*
+ * readphylip.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/24/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+using namespace std;
+
+#include "globaldata.hpp"
+#include "sequencedb.h"
+#include "mothur.h"
+
+/**********************************************************************************/
+
+class ReadPhylip {
+
+ public:
+ ReadPhylip(string);
+ ~ReadPhylip();
+ void read();
+ SequenceDB* getDB();
+
+ private:
+ GlobalData* globaldata;
+ string phylipFile;
+ ifstream filehandle;
+ SequenceDB sequencedb;
+ int readOk; // readOk = 0 means success, readOk = 1 means error(s).
+
+ bool isSeq(string);
+};
+
+#endif
\ No newline at end of file
#include "sequence.hpp"
-//********************************************************************************************************************
-
-Sequence::Sequence(){}
+/***********************************************************************/
-//********************************************************************************************************************
-
-Sequence::Sequence(ifstream& fastaFile){
-
- string accession;
- fastaFile >> accession;
- setName(accession);
+Sequence::Sequence() {}
- char letter;
- string sequence;
-
- while(fastaFile && letter != '>'){
-
- letter = fastaFile.get();
-
- if(isalpha(letter)){
-
- sequence += letter;
-
- }
-
- }
- fastaFile.putback(letter);
+/***********************************************************************/
- if(sequence.find_first_of('-') != string::npos){
+Sequence::Sequence(string newName, string sequence) {
+ name = newName;
+ if(sequence.find_first_of('-') != string::npos) {
setAligned(sequence);
}
setUnaligned(sequence);
}
-
//********************************************************************************************************************
-string Sequence::convert2ints(){
+string Sequence::convert2ints() {
if(unaligned == "") { /* need to throw an error */ }
string processed;
- for(int i=0;i<unaligned.length();i++){
+ for(int i=0;i<unaligned.length();i++) {
if(toupper(unaligned[i]) == 'A') { processed += '0'; }
else if(toupper(unaligned[i]) == 'C') { processed += '1'; }
else if(toupper(unaligned[i]) == 'G') { processed += '2'; }
//********************************************************************************************************************
-void Sequence::setName(string seqName){
+void Sequence::setName(string seqName) {
if(seqName[0] == '>') { name = seqName.substr(1); }
else { name = seqName; }
}
void Sequence::setUnaligned(string sequence){
- if(sequence.find_first_of('-') != string::npos){
+ if(sequence.find_first_of('-') != string::npos) {
string temp = "";
- for(int j=0;j<sequence.length();j++){
+ for(int j=0;j<sequence.length();j++) {
if(isalpha(sequence[j])) { temp += sequence[j]; }
}
unaligned = temp;
}
- else{
+ else {
unaligned = sequence;
}
//********************************************************************************************************************
-string Sequence::getSeqName(){
+string Sequence::getName(){
return name;
}
}
//********************************************************************************************************************
+
+int Sequence::getLength(){
+ if(unaligned.length() > aligned.length())
+ return unaligned.length();
+ return aligned.length();
+}
+
+//********************************************************************************************************************
+
+void Sequence::printSequence(ostream& out){
+ string toPrint = unaligned;
+ if(aligned.length() > unaligned.length())
+ toPrint = aligned;
+ out << ">" << name << "\n" << toPrint << "\n";
+}
+
+//********************************************************************************************************************
class Sequence {
public:
Sequence();
- Sequence(ifstream&);
+ Sequence(string, string);
+
void setName(string);
void setUnaligned(string);
void setPairwise(string);
void setAligned(string);
+ void setLength();
+
string convert2ints();
- string getSeqName();
+ string getName();
string getAligned();
string getPairwise();
string getUnaligned();
+ int getLength();
+ void printSequence(ostream&);
+
private:
string name;
string unaligned;
- string pairwise;
string aligned;
+ string pairwise;
+ int length;
+ int lengthAligned;
};
/**************************************************************************************************/
--- /dev/null
+/*
+ * sequencedb.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/13/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "sequencedb.h"
+#include "sequence.hpp"
+#include "mothur.h"
+#include "calculator.h"
+
+
+/***********************************************************************/
+
+SequenceDB::SequenceDB() {}
+
+/***********************************************************************/
+
+SequenceDB::SequenceDB(int newSize) {
+ data.resize(newSize);
+}
+
+/***********************************************************************/
+
+SequenceDB::SequenceDB(ifstream&) {}
+
+/***********************************************************************/
+
+int SequenceDB::getNumSeqs() {
+ return data.size();
+}
+
+/***********************************************************************/
+
+void SequenceDB::set(int index, string newUnaligned) {
+ Sequence newSeq(data[index].getName(), newUnaligned);
+ data[index] = newSeq;
+}
+
+/***********************************************************************/
+
+void SequenceDB::set(int index, Sequence newSeq) {
+ data[index] = newSeq;
+}
+
+/***********************************************************************/
+
+Sequence SequenceDB::get(int index) {
+ return data[index];
+}
+
+/***********************************************************************/
+
+void SequenceDB::changeSize(int newSize) {
+ data.resize(newSize);
+}
+
+/***********************************************************************/
+
+void SequenceDB::clear() {
+ data.clear();
+}
+
+/***********************************************************************/
+
+int SequenceDB::size() {
+ return data.size();
+}
+
+/***********************************************************************/
+
+void SequenceDB::print(ostream& out) {
+ for(int i = 0; i < data.size(); i++)
+ data[i].printSequence(out);
+}
+
+/***********************************************************************/
+
+void SequenceDB::add(Sequence newSequence) {
+ try {
+ data.push_back(newSequence);
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the RAbundVector class Function push_back. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the RAbundVector class function push_back. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
\ No newline at end of file
--- /dev/null
+#ifndef SEQUENCEDB_H
+#define SEQUENCEDB_H
+
+/*
+ * sequencedb.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/13/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+
+using namespace std;
+
+#include "sequence.hpp"
+#include "calculator.h"
+
+
+
+
+class SequenceDB {
+
+public:
+ SequenceDB();
+ SequenceDB(int); //makes data that size
+ SequenceDB(ifstream&); //reads file to fill data
+// ~SequenceDB(); //loops through data and delete each sequence
+
+ int getNumSeqs();
+
+ void set(int, string); //unaligned - should also set length
+ void set(int, Sequence); //unaligned - should also set length
+ Sequence get(int); //returns sequence name at that location
+ void add(Sequence); //adds unaligned sequence
+ void changeSize(int); //resizes data
+ void clear(); //clears data - remeber to loop through and delete the sequences inside or you will have a memory leak
+ int size(); //returns datas size
+ void print(ostream&); //loops through data using sequence class print
+
+private:
+ vector<Sequence> data;
+
+};
+
+#endif
--- /dev/null
+/*
+ * sharedjackknife.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 3/30/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "sharedjackknife.h"
+
+/***************************************************************************************
+
+***************************************************************************************/
+double SharedJackknife::simpson(int abunds[], double numInd, int numBins){
+ double denom = numInd*(numInd-1);
+ double sum = 0;
+ for(int i = 0; i < numBins; i++)
+ sum += (double)abunds[i]*((double)abunds[i]-1)/denom;
+
+ return sum;
+}
+
+/*****************************************************************************************/
+
+double* SharedJackknife::jackknife(){
+ int numBins = groups.at(0)->getNumBins()-1;
+ int cArray[numBins];
+ for(int i = 0; i < numBins; i++)
+ cArray[i] = 0;
+
+ double numInd = 0;
+ for(int i = 0; i < numGroups; i++)
+ for(int j = 0; j < numBins; j++) {
+ int curAbund = groups.at(i)->get(j+1).abundance;
+ cArray[j] += curAbund;
+ numInd += (double)curAbund;
+ }
+
+ double baseD = 1/simpson(cArray, numInd, numBins);
+
+ double pseudoVals[numBins];
+ double jackknifeEstimate = 0;
+ for(int i = 0; i < numGroups; i++) {
+ for(int j = 0; j < numBins-1; j++) {
+ int abundDiff = -groups.at(i)->get(j+1).abundance;
+ if(i > 0)
+ abundDiff += groups.at(i-1)->get(j+1).abundance;
+
+ cArray[j] += abundDiff;
+ numInd += abundDiff;
+ }
+
+ double curD = 1/simpson(cArray, numInd, numBins);
+ pseudoVals[i] = (double)numGroups*(baseD - curD) + curD;
+ jackknifeEstimate += pseudoVals[i];
+ }
+ jackknifeEstimate /= (double)numGroups;
+
+ double variance = 0;
+ for(int i = 0; i < numGroups; i++)
+ variance += pow(pseudoVals[i]-jackknifeEstimate, 2);
+
+ variance /= (double)numGroups*((double)numGroups-1);
+ double stErr = sqrt(variance);
+ TDTable table;
+ double confLimit = 0;
+ if(numGroups <= 30)
+ confLimit = table.getConfLimit(numGroups-1, 1);
+ else
+ confLimit = 1.645;
+
+ confLimit *= stErr;
+
+ double* rdata = new double[3];
+ rdata[0] = baseD;
+ rdata[1] = jackknifeEstimate - confLimit;
+ rdata[2] = jackknifeEstimate + confLimit;
+
+ return rdata;
+}
+
+/************************************************************************************************
+************************************************************************************************/
+
+EstOutput SharedJackknife::getValues(vector<SharedRAbundVector*> vectorShared){ //Fix this for collect, mistake was that it was made with summary in mind.
+ try {
+ SharedRAbundVector* shared1 = vectorShared[0];
+ SharedRAbundVector* shared2 = vectorShared[1];
+ if(numGroups == -1) {
+ globaldata = GlobalData::getInstance();
+ numGroups = globaldata->Groups.size();
+ }
+
+ if(callCount == numGroups*(numGroups-1)/2) {
+ currentCallDone = true;
+ callCount = 0;
+ }
+ callCount++;
+
+ if(currentCallDone) {
+ groups.clear();
+ currentCallDone = false;
+ }
+
+ if(groups.size() != numGroups) {
+ if(groups.size() == 0)
+ groups.push_back(shared1);
+ groups.push_back(shared2);
+ }
+
+ if(groups.size() == numGroups && callCount < numGroups) {
+ data.resize(3,0);
+
+ double* rdata = jackknife();
+ data[0] = rdata[0];
+ data[1] = rdata[1];
+ data[2] = rdata[2];
+
+ //cout << "sT = " << data[0] << " lower confLimit = " << data[1] << " upper confLimit = " << data[2] << "\n";
+ if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
+ if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
+ if (isnan(data[2]) || isinf(data[0])) { data[2] = 0; }
+
+ /*for(int i = 0; i < groups.size(); i++)
+ cout << groups.at(i)->getGroup() << " ";
+ cout << "\n";
+ cout << groups.size() << " " << data[0] << " " << data[1] << " " << data[2] << "\n\n";*/
+
+ return data;
+ }
+
+ data.resize(3,0);
+ data[0] = 0;
+ data[1] = 0;
+ data[2] = 0;
+
+ if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
+ if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
+ if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
+ return data;
+ }
+
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the SharedJackknife class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the SharedJackknife class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+/***********************************************************************/
+
--- /dev/null
+#ifndef SHAREDJACKKNIFE_H
+#define SHAREDJACKKNIFE_H
+
+/*
+ * sharedjackknife.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 3/30/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "calculator.h"
+#include "globaldata.hpp"
+
+/*This class implements the SharedJackknife estimator.
+It is a child of the calculator class.*/
+
+/***********************************************************************/
+
+class SharedJackknife : public Calculator {
+
+public:
+ SharedJackknife() : numGroups(-1), callCount(0), count(0), currentCallDone(true), Calculator("sharedjackknife", 3, false) {};
+ EstOutput getValues(SAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>);
+
+private:
+ GlobalData* globaldata;
+ int numGroups, callCount, count;
+ bool currentCallDone;
+ vector<SharedRAbundVector*> groups;
+ double simpson(int[], double, int);
+ double* jackknife();
+};
+
+/***********************************************************************/
+
+#endif
+
EstOutput KSTest::getValues(vector<SharedRAbundVector*> shared){
try {
- data.resize(2,0);
-
+ data.resize(3,0);
+
//Must return shared1 and shared2 to original order at conclusion of kstest
vector <individual> initData1 = shared[0]->getData();
vector <individual> initData2 = shared[1]->getData();
data[0] = DStatistic;
data[1] = critVal;
+ data[2] = 0;
+
if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
--- /dev/null
+/*
+ * sharedmarczewski.cpp
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/8/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "sharedmarczewski.h"
+
+EstOutput SharedMarczewski::getValues(vector<SharedRAbundVector*> vectorShared){
+ try {
+ SharedRAbundVector* shared1 = vectorShared[0];
+ SharedRAbundVector* shared2 = vectorShared[1];
+
+ data.resize(1,0);
+
+ double a = 0;
+ double b = 0;
+ double c = 0;
+ for(int i = 1; i < shared1->size(); i++)
+ {
+ int abund1 = shared1->get(i).abundance;
+ int abund2 = shared2->get(i).abundance;
+
+ if(abund1 > 0 && abund2 > 0)
+ a++;
+ else if(abund1 > 0 && abund2 == 0)
+ b++;
+ else if(abund1 == 0 && abund2 > 0)
+ c++;
+ }
+ data[0] = (b+c)/(a+b+c);
+
+ if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
+
+ return data;
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the SharedMarczewski class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the SharedMarczewski class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
+/***********************************************************************/
--- /dev/null
+#ifndef SHAREDMARCZEWSKI_H
+#define SHAREDMARCZEWSKI_H
+
+/*
+ * sharedmarczewski.h
+ * Mothur
+ *
+ * Created by Thomas Ryabin on 4/8/09.
+ * Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "calculator.h"
+
+/***********************************************************************/
+
+class SharedMarczewski : public Calculator {
+
+public:
+ SharedMarczewski() : Calculator("sharedmarczewski", 1, false) {};
+ EstOutput getValues(SAbundVector*) {return data;};
+ EstOutput getValues(vector<SharedRAbundVector*>);
+private:
+};
+
+/***********************************************************************/
+
+#endif
\ No newline at end of file
#include "qstat.h"
#include "bergerparker.h"
#include "bstick.h"
+#include "goodscoverage.h"
#include "coverage.h"
//**********************************************************************************************************************
sumCalculators.push_back(new Bootstrap());
}else if (globaldata->Estimators[i] == "nseqs") {
sumCalculators.push_back(new NSeqs());
+ }else if (globaldata->Estimators[i] == "goodscoverage") {
+ sumCalculators.push_back(new GoodsCoverage());
}
}
}
#include "sharedlennon.h"
#include "sharedmorisitahorn.h"
#include "sharedbraycurtis.h"
+#include "sharedjackknife.h"
+#include "sharedwhittaker.h"
//**********************************************************************************************************************
single["logseries"] = "logseries";
single["qstat"] = "qstat";
single["bstick"] = "bstick";
+ single["goodscoverage"] = "goodscoverage";
single["nseqs"] = "nseqs";
single["coverage"] = "coverage";
-
single["default"] = "default";
}
catch(exception& e) {
summary["qstat"] = "qstat";
summary["bstick"] = "bstick";
summary["nseqs"] = "nseqs";
+ summary["goodscoverage"]= "goodscoverage";
summary["coverage"] = "coverage";
summary["default"] = "default";
}
commands["read.dist"] = "read.dist";
commands["read.otu"] = "read.otu";
commands["read.tree"] = "read.tree";
+ commands["read.seqs"] = "read.seqs";
commands["bin.seqs"] = "bin.seqs";
commands["get.oturep"] = "get.oturep";
commands["cluster"] = "cluster";
commands["bootstrap.shared"] = "bootstrap.shared";
commands["concensus"] = "concensus";
commands["help"] = "help";
+ commands["filter.seqs"] = "filter.seqs";
commands["quit"] = "quit";
bool valid = false;
vector<string> cParams = commandParameters[command];
int numParams = cParams.size();
- for(int i = 0; i < numParams; i++)
- {
- if(cParams.at(i).compare(parameter) == 0)
- {
+ for(int i = 0; i < numParams; i++) {
+ if(cParams.at(i).compare(parameter) == 0) {
valid = true;
i = numParams;
}
}
- if(!valid)
- {
+ if(!valid) {
cout << "'" << parameter << "' is not a valid parameter for the " << command << " command.\n";
cout << "The valid paramters for the " << command << " command are: ";
for(int i = 0; i < numParams-1; i++)
vector<string> values;
splitAtDash(value, values);
- for(int i = 0; i < values.size(); i++)
- {
+ for(int i = 0; i < values.size(); i++) {
value = values.at(i);
valid = convertTest(value, pVal);
Special Cases
*********************************************************************************************************/
- if(parameter.compare("precision") == 0)
- {
+ if(parameter.compare("precision") == 0) {
double logNum = log10((double)pVal);
double diff = (double)((int)logNum - logNum);
- if(diff != 0)
- {
+ if(diff != 0) {
cout << "The precision parameter can only take powers of 10 as a value (e.g. 10,1000,1000, etc.)\n";
return false;
}
c = 0;
else if(range.at(4).compare("only") == 0)
c = 1;
- else
- {
+ else {
cout << "The range can only be 'between' or 'only' the bounding numbers.\n";
return false;
}
d = 0;
else if(range.at(0).compare(">=") == 0 || range[3].compare("=>") == 0)
d = 1;
- else
- {
+ else {
cout << "The parameter value can only be '>', '>=', or '=>' the lower bounding number.\n";
return false;
}
e = 0;
else if(range.at(2).compare("<=") == 0 || range[4].compare("=<") == 0)
e = 1;
- else
- {
+ else {
cout << "The parameter value can only be '<', '<=', or '=<' the upper bounding number.\n";
return false;
}
bool b0 = pVal < b;
bool b1 = pVal <= b;
- if(c != 1)
- {
- if(a != piSentinel && b == piSentinel)
- {
+ if(c != 1) {
+ if(a != piSentinel && b == piSentinel) {
if(d == 0)
valid = a0;
else
valid = a1;
}
- else if(a == piSentinel && b != piSentinel)
- {
+ else if(a == piSentinel && b != piSentinel) {
if(e == 0)
valid = b0;
else
valid = b1;
}
- else
- {
+ else {
if(d == 0 && e == 0)
valid = (a0 && b0);
else if(d == 0 && e == 1)
valid = (a1 && b1);
}
}
- else
- {
+ else {
if(a != piSentinel && b == piSentinel)
valid = (pVal == a);
else if(a == piSentinel && b != piSentinel)
}
- if(!valid)
- {
+ if(!valid) {
cout << "The '" << parameter << "' parameter needs to be ";
if(c == 1)
cout << "either '" << a << "' or '" << b << "'.\n";
- else
- {
- if(a != piSentinel)
- {
+ else {
+ if(a != piSentinel) {
cout << ">";
if(d != 0)
cout << "=";
cout << " '" << a << "'";
}
if(b == piSentinel)
- cout << ".\n";
+ cout << "'.\n";
else if(a != piSentinel)
cout << " and ";
- if(b != piSentinel)
- {
+ if(b != piSentinel) {
cout << "<";
if(e != 0)
cout << "=";
- cout << " '" << b << ".\n";
+ cout << " '" << b << "'.\n";
}
}
return false;
string readtreeArray[] = {"tree","group"};
commandParameters["read.tree"] = addParameters(readtreeArray, sizeof(readtreeArray)/sizeof(string));
+ string readseqsArray[] = {"fasta","phylip","clustal","nexus","line"};
+ commandParameters["read.seqs"] = addParameters(readseqsArray, sizeof(readseqsArray)/sizeof(string));
+
string clusterArray[] = {"cutoff","precision","method"};
commandParameters["cluster"] = addParameters(clusterArray, sizeof(clusterArray)/sizeof(string));
string heatmapArray[] = {"groups","line","label","sorted","scale"};
commandParameters["heatmap"] = addParameters(heatmapArray, sizeof(heatmapArray)/sizeof(string));
+
+ string filterseqsArray[] = {"trump", "soft", "filter"};
+ commandParameters["filter.seqs"] = addParameters(filterseqsArray, sizeof(filterseqsArray)/sizeof(string));
string vennArray[] = {"groups","line","label","calc"};
commandParameters["venn"] = addParameters(vennArray, sizeof(vennArray)/sizeof(string));
string jumbleArray[] = {">","0", "<","1", "only"};
parameterRanges["jumble"] = addParameters(jumbleArray, rangeSize);
- string freqArray[] = {">","1", "<","NA", "between"};
+ string freqArray[] = {">=","1", "<","NA", "between"};
parameterRanges["freq"] = addParameters(freqArray, rangeSize);
//string lineArray[] = {">=","1", "<","NA", "between"};
string abundArray[] = {">=","5", "<","NA", "between"};
parameterRanges["abund"] = addParameters(abundArray, rangeSize);
+
+ string softArray[] = {">=","0", "<=","100", "between"};
+ parameterRanges["soft"] = addParameters(softArray, rangeSize);
}
catch(exception& e) {
cout << "Standard Error: " << e.what() << " has occurred in the ValidParameters class Function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
projects / mothur.git / commitdiff