From: pschloss Date: Thu, 23 Apr 2009 20:01:31 +0000 (+0000) Subject: more calculator edits and added coverage and whittaker X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=477e76a8a79b60f6cd4253324dd830bdea25e3e9 more calculator edits and added coverage and whittaker --- diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 1c8df76..6fd4a06 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -115,6 +115,8 @@ 7E412F490F8D21B600381DD0 /* slibshuff.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7E412F480F8D21B600381DD0 /* slibshuff.cpp */; }; 7E412FEA0F8D3E2C00381DD0 /* libshuff.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7E412FE90F8D3E2C00381DD0 /* libshuff.cpp */; }; 7E4130F80F8E58FA00381DD0 /* dlibshuff.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7E4130F60F8E58FA00381DD0 /* dlibshuff.cpp */; }; + 7EC3D4550FA0FFF900338DA5 /* coverage.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7EC3D4500FA0FFF900338DA5 /* coverage.cpp */; }; + 7EC3D4560FA0FFF900338DA5 /* whittaker.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7EC3D4530FA0FFF900338DA5 /* whittaker.cpp */; }; 8DD76F6A0486A84900D96B5E /* Mothur.1 in CopyFiles */ = {isa = PBXBuildFile; fileRef = C6859E8B029090EE04C91782 /* Mothur.1 */; }; A70B53AA0F4CD7AD0064797E /* getgroupcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A70B53A40F4CD7AD0064797E /* getgroupcommand.cpp */; }; A70B53AB0F4CD7AD0064797E /* getlabelcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A70B53A60F4CD7AD0064797E /* getlabelcommand.cpp */; }; @@ -123,7 +125,6 @@ EB1216E50F61ACFB004A865F /* bstick.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB1216E40F61ACFB004A865F /* bstick.cpp */; }; EB1217230F61C9AC004A865F /* sharedkstest.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB1217220F61C9AC004A865F /* sharedkstest.cpp */; }; EB6E68190F5F1C780044E17B /* qstat.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB6E68180F5F1C780044E17B /* qstat.cpp */; }; - EB6F015B0F6AC1670048081A /* sharedbdiversity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB6F015A0F6AC1670048081A /* sharedbdiversity.cpp */; }; EB9303EB0F534D9400E8EF26 /* logsd.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB9303EA0F534D9400E8EF26 /* logsd.cpp */; }; EB9303F90F53517300E8EF26 /* geom.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB9303F80F53517300E8EF26 /* geom.cpp */; }; /* End PBXBuildFile section */ @@ -369,6 +370,11 @@ 7E412FE90F8D3E2C00381DD0 /* libshuff.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = libshuff.cpp; sourceTree = ""; }; 7E4130F60F8E58FA00381DD0 /* dlibshuff.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = dlibshuff.cpp; sourceTree = ""; }; 7E4130F70F8E58FA00381DD0 /* dlibshuff.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = dlibshuff.h; sourceTree = ""; }; + 7EC3D4500FA0FFF900338DA5 /* coverage.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = coverage.cpp; sourceTree = ""; }; + 7EC3D4510FA0FFF900338DA5 /* coverage.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = coverage.h; sourceTree = ""; }; + 7EC3D4520FA0FFF900338DA5 /* sharedanderberg.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedanderberg.h; sourceTree = ""; }; + 7EC3D4530FA0FFF900338DA5 /* whittaker.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = whittaker.cpp; sourceTree = ""; }; + 7EC3D4540FA0FFF900338DA5 /* whittaker.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = whittaker.h; sourceTree = ""; }; 8DD76F6C0486A84900D96B5E /* mothur */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = mothur; sourceTree = BUILT_PRODUCTS_DIR; }; A70B53A40F4CD7AD0064797E /* getgroupcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = getgroupcommand.cpp; sourceTree = ""; }; A70B53A50F4CD7AD0064797E /* getgroupcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = getgroupcommand.h; sourceTree = ""; }; @@ -385,8 +391,6 @@ EB1217220F61C9AC004A865F /* sharedkstest.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedkstest.cpp; sourceTree = ""; }; EB6E68170F5F1C780044E17B /* qstat.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = qstat.h; sourceTree = ""; }; EB6E68180F5F1C780044E17B /* qstat.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = qstat.cpp; sourceTree = ""; }; - EB6F01590F6AC1670048081A /* sharedbdiversity.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedbdiversity.h; sourceTree = ""; }; - EB6F015A0F6AC1670048081A /* sharedbdiversity.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedbdiversity.cpp; sourceTree = ""; }; EB9303E90F534D9400E8EF26 /* logsd.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = logsd.h; sourceTree = ""; }; EB9303EA0F534D9400E8EF26 /* logsd.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = logsd.cpp; sourceTree = ""; }; EB9303F70F53517300E8EF26 /* geom.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = geom.h; sourceTree = ""; }; @@ -537,10 +541,13 @@ 37D928220F21331F001D4494 /* sharedace.cpp */, 37519A690F80E6EB00FED5E8 /* sharedanderbergs.h */, 37519A6A0F80E6EB00FED5E8 /* sharedanderbergs.cpp */, + 7EC3D4500FA0FFF900338DA5 /* coverage.cpp */, + 7EC3D4510FA0FFF900338DA5 /* coverage.h */, + 7EC3D4520FA0FFF900338DA5 /* sharedanderberg.h */, + 7EC3D4530FA0FFF900338DA5 /* whittaker.cpp */, + 7EC3D4540FA0FFF900338DA5 /* whittaker.h */, 375873FB0F7D64DA0040F377 /* sharedbraycurtis.h */, 375873FC0F7D64DA0040F377 /* sharedbraycurtis.cpp */, - EB6F01590F6AC1670048081A /* sharedbdiversity.h */, - EB6F015A0F6AC1670048081A /* sharedbdiversity.cpp */, 37D928250F21331F001D4494 /* sharedchao1.h */, 37D928240F21331F001D4494 /* sharedchao1.cpp */, 37D928290F21331F001D4494 /* sharedjabund.h */, @@ -851,7 +858,6 @@ EB1216880F619B83004A865F /* bergerparker.cpp in Sources */, EB1216E50F61ACFB004A865F /* bstick.cpp in Sources */, EB1217230F61C9AC004A865F /* sharedkstest.cpp in Sources */, - EB6F015B0F6AC1670048081A /* sharedbdiversity.cpp in Sources */, 375873EC0F7D64520040F377 /* fullmatrix.cpp in Sources */, 375873EF0F7D646F0040F377 /* heatmap.cpp in Sources */, 375873F20F7D64800040F377 /* heatmapcommand.cpp in Sources */, @@ -877,6 +883,8 @@ 375AA1390F9E433D008EF9B8 /* readcolumn.cpp in Sources */, 375AA13A0F9E433D008EF9B8 /* readotu.cpp in Sources */, 375AA13B0F9E433D008EF9B8 /* readphylip.cpp in Sources */, + 7EC3D4550FA0FFF900338DA5 /* coverage.cpp in Sources */, + 7EC3D4560FA0FFF900338DA5 /* whittaker.cpp in Sources */, ); runOnlyForDeploymentPostprocessing = 0; }; diff --git a/bootstrapsharedcommand.cpp b/bootstrapsharedcommand.cpp index 56611a8..5c2e171 100644 --- a/bootstrapsharedcommand.cpp +++ b/bootstrapsharedcommand.cpp @@ -35,7 +35,7 @@ BootSharedCommand::BootSharedCommand(){ if (validCalculator->isValidCalculator("boot", globaldata->Estimators[i]) == true) { if (globaldata->Estimators[i] == "jabund") { treeCalculators.push_back(new JAbund()); - }else if (globaldata->Estimators[i] == "sorensonabund") { + }else if (globaldata->Estimators[i] == "sorabund") { treeCalculators.push_back(new SorAbund()); }else if (globaldata->Estimators[i] == "jclass") { treeCalculators.push_back(new Jclass()); diff --git a/bstick.cpp b/bstick.cpp index 7f12d1d..eabf087 100644 --- a/bstick.cpp +++ b/bstick.cpp @@ -46,7 +46,7 @@ RAbundVector BStick::getRAbundVector(SAbundVector* rank){ /***************************************************************************/ EstOutput BStick::getValues(SAbundVector* rank){ try { - data.resize(2,0); + data.resize(3,0); rdata = getRAbundVector(rank); double numInd = (double)rdata.getNumSeqs(); double numSpec = (double)rdata.getNumBins(); @@ -64,27 +64,16 @@ EstOutput BStick::getValues(SAbundVector* rank){ maxDiff = diff; } - double DStatistic = maxDiff/numInd; - double critVal = 0; - /*cout << "BStick:\n"; - cout << "D-Statistic = " << DStatistic << "\n"; - cout << "Critical value for 95% confidence interval = ";*/ - if(rdata.size() > 20) - { - critVal = .886/sqrt(rdata.size()); - } - else - { - KOSTable table; - critVal = table.getConfLimit(numSpec); - } + 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";*/ - data[0] = DStatistic; - data[1] = critVal; 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; } diff --git a/collectsharedcommand.cpp b/collectsharedcommand.cpp index ad392e7..bc0609e 100644 --- a/collectsharedcommand.cpp +++ b/collectsharedcommand.cpp @@ -20,7 +20,7 @@ #include "sharedthetayc.h" #include "sharedthetan.h" #include "sharedkstest.h" -#include "sharedbdiversity.h" +#include "whittaker.h" #include "sharednseqs.h" #include "sharedochiai.h" #include "sharedanderbergs.h" @@ -54,7 +54,7 @@ CollectSharedCommand::CollectSharedCommand(){ cDisplays.push_back(new CollectDisplay(new SharedAce(), new SharedOneColumnFile(fileNameRoot+"shared.ace"))); }else if (globaldata->Estimators[i] == "jabund") { cDisplays.push_back(new CollectDisplay(new JAbund(), new SharedOneColumnFile(fileNameRoot+"jabund"))); - }else if (globaldata->Estimators[i] == "sorensonabund") { + }else if (globaldata->Estimators[i] == "sorabund") { cDisplays.push_back(new CollectDisplay(new SorAbund(), new SharedOneColumnFile(fileNameRoot+"sorabund"))); }else if (globaldata->Estimators[i] == "jclass") { cDisplays.push_back(new CollectDisplay(new Jclass(), new SharedOneColumnFile(fileNameRoot+"jclass"))); @@ -70,8 +70,8 @@ CollectSharedCommand::CollectSharedCommand(){ cDisplays.push_back(new CollectDisplay(new ThetaN(), new SharedOneColumnFile(fileNameRoot+"thetan"))); }else if (globaldata->Estimators[i] == "kstest") { cDisplays.push_back(new CollectDisplay(new KSTest(), new SharedOneColumnFile(fileNameRoot+"kstest"))); - }else if (globaldata->Estimators[i] == "bdiversity") { - cDisplays.push_back(new CollectDisplay(new BDiversity(), new SharedOneColumnFile(fileNameRoot+"bdiversity"))); + }else if (globaldata->Estimators[i] == "whittaker") { + 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") { diff --git a/coverage.cpp b/coverage.cpp new file mode 100644 index 0000000..a93cbf3 --- /dev/null +++ b/coverage.cpp @@ -0,0 +1,33 @@ +/* + * coverage.cpp + * Mothur + * + * Created by Pat Schloss on 4/22/09. + * Copyright 2009 Patrick D. Schloss. All rights reserved. + * + */ + +#include "coverage.h" + +/***********************************************************************/ +EstOutput Coverage::getValues(SAbundVector* rank){ + + try { + data.resize(1,0); + + data[0] = 1. - rank->get(1) / (double)rank->getNumSeqs(); + + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +}; + + +/***********************************************************************/ diff --git a/coverage.h b/coverage.h new file mode 100644 index 0000000..c17f8e7 --- /dev/null +++ b/coverage.h @@ -0,0 +1,31 @@ +#ifndef COVERAGE_H +#define COVERAGE_H + +/* + * coverage.h + * Mothur + * + * Created by Pat Schloss on 4/22/09. + * Copyright 2009 Patrick D. Schloss. All rights reserved. + * + */ + +#include "calculator.h" + +/* This class implements the coverage estimator on single group. + It is a child of the calculator class. */ + +/***********************************************************************/ + +class Coverage : public Calculator { + +public: + Coverage() : Calculator("coverage", 1) {}; + EstOutput getValues(SAbundVector*); + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;}; +}; + + +/***********************************************************************/ + +#endif diff --git a/geom.cpp b/geom.cpp index 31580db..755c1ec 100644 --- a/geom.cpp +++ b/geom.cpp @@ -75,26 +75,18 @@ EstOutput Geom::getValues(SAbundVector* rank){ } - double DStatistic = maxDiff/numInd; - double critVal = 0; + /*cout << "Geom:\n"; cout << "D-Statistic = " << DStatistic << "\n"; cout << "Critical value for 95% confidence interval = ";*/ - if(rdata.size() > 20) - { - data[1] = 1.031/sqrt(rdata.size()); - data[2] = 0.886/sqrt(rdata.size()); - } - else - { - KOSTable table; - critVal = table.getConfLimit(numSpec); - } + + data[0] = maxDiff/numInd; + data[1] = 0.886/sqrt(numSpec); + data[2] = 1.031/sqrt(numSpec); + /*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";*/ - - data[0] = DStatistic; - + 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; } diff --git a/geom.h b/geom.h index 05fc0d7..f88cb4a 100644 --- a/geom.h +++ b/geom.h @@ -11,7 +11,7 @@ */ #include "calculator.h" -/* This class implements the LogSD estimator on single group. +/* This class implements the geometric estimator on single group. It is a child of the calculator class. */ /***********************************************************************/ diff --git a/globaldata.cpp b/globaldata.cpp index 714b5ad..312b314 100644 --- a/globaldata.cpp +++ b/globaldata.cpp @@ -173,7 +173,7 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ splitAtDash(calc, Estimators); } if (commandName == "collect.shared") { - if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan"; } + if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan"; } Estimators.clear(); splitAtDash(calc, Estimators); } @@ -183,7 +183,7 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ splitAtDash(calc, Estimators); } if (commandName == "summary.shared") { - if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan"; } + if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan"; } Estimators.clear(); splitAtDash(calc, Estimators); } diff --git a/helpcommand.cpp b/helpcommand.cpp index 6cf2134..ac01330 100644 --- a/helpcommand.cpp +++ b/helpcommand.cpp @@ -86,8 +86,8 @@ int HelpCommand::execute(){ cout << "The collect.shared command parameters are label, line, freq, jumble, calc and groups. No parameters are required, but you may not use " << "\n"; cout << "both the line and label parameters at the same time. The collect.shared command should be in the following format: " << "\n"; cout << "collect.shared(label=yourLabel, line=yourLines, freq=yourFreq, jumble=yourJumble, calc=yourEstimators, groups=yourGroups)." << "\n"; - cout << "Example collect.shared(label=unique-.01-.03, line=0-5-10, freq=10, jumble=1, groups=B-C, calc=sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n"; - cout << "The default values for jumble is 1 (meaning jumble, if it’s set to 0 then it will not jumble), freq is 100 and calc are sharedsobs-sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan." << "\n"; + cout << "Example collect.shared(label=unique-.01-.03, line=0-5-10, freq=10, jumble=1, groups=B-C, calc=sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n"; + cout << "The default values for jumble is 1 (meaning jumble, if it’s set to 0 then it will not jumble), freq is 100 and calc are sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan." << "\n"; cout << "The default value for groups is all the groups in your groupfile." << "\n"; validCalcs->printCalc("shared", cout); cout << "The label and line parameters are used to analyze specific lines in your input." << "\n"; @@ -150,9 +150,9 @@ int HelpCommand::execute(){ cout << "The summary.shared command parameters are label, line, jumble and calc. No parameters are required, but you may not use " << "\n"; cout << "both the line and label parameters at the same time. The summary.shared command should be in the following format: " << "\n"; cout << "summary.shared(label=yourLabel, line=yourLines, jumble=yourJumble, calc=yourEstimators, groups=yourGroups)." << "\n"; - cout << "Example summary.shared(label=unique-.01-.03, line=0,5,10, jumble=1, groups=B-C, calc=sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n"; + cout << "Example summary.shared(label=unique-.01-.03, line=0,5,10, jumble=1, groups=B-C, calc=sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan)." << "\n"; validCalcs->printCalc("sharedsummary", cout); - cout << "The default value for jumble is 1 (meaning jumble, if it’s set to 0 then it will not jumble) and calc is sharedsobs-sharedchao-sharedace-jabund-sorensonabund-jclass-sorclass-jest-sorest-thetayc-thetan" << "\n"; + cout << "The default value for jumble is 1 (meaning jumble, if it’s set to 0 then it will not jumble) and calc is sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan" << "\n"; cout << "The default value for groups is all the groups in your groupfile." << "\n"; cout << "The label and line parameters are used to analyze specific lines in your input." << "\n"; cout << "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed. You must enter at least 2 valid groups." << "\n"; @@ -232,7 +232,7 @@ int HelpCommand::execute(){ cout << "The groups parameter allows you to specify which of the groups in your groupfile you would like included used." << "\n"; cout << "The group names are separated by dashes. The line and label allow you to select what distance levels you would like trees created for, and are also separated by dashes." << "\n"; cout << "The tree.shared command should be in the following format: tree.shared(groups=yourGroups, calc=yourCalcs, line=yourLines, label=yourLabels)." << "\n"; - cout << "Example tree.shared(groups=A-B-C, line=1-3-5, calc=jabund-sorensonabund)." << "\n"; + cout << "Example tree.shared(groups=A-B-C, line=1-3-5, calc=jabund-sorabund)." << "\n"; cout << "The default value for groups is all the groups in your groupfile." << "\n"; cout << "There is no default value for calc." << "\n"; validCalcs->printCalc("treegroup", cout); @@ -244,7 +244,7 @@ int HelpCommand::execute(){ cout << "The groups parameter allows you to specify which of the groups in your groupfile you would like included used." << "\n"; cout << "The group names are separated by dashes. The line and label allow you to select what distance levels you would like trees created for, and are also separated by dashes." << "\n"; cout << "The bootstrap.shared command should be in the following format: bootstrap.shared(groups=yourGroups, calc=yourCalcs, line=yourLines, label=yourLabels, iters=yourIters)." << "\n"; - cout << "Example bootstrap.shared(groups=A-B-C, line=1-3-5, calc=jabund-sorensonabund, iters=100)." << "\n"; + cout << "Example bootstrap.shared(groups=A-B-C, line=1-3-5, calc=jabund-sorabund, iters=100)." << "\n"; cout << "The default value for groups is all the groups in your groupfile." << "\n"; cout << "There is no default value for calc. The default for iters is 1000." << "\n"; validCalcs->printCalc("boot", cout); diff --git a/logsd.cpp b/logsd.cpp index 14d2e44..e959ab8 100644 --- a/logsd.cpp +++ b/logsd.cpp @@ -78,17 +78,21 @@ EstOutput LogSD::getValues(SAbundVector* rank){ if(diff > maxDiff) maxDiff = diff; } - - double DStatistic = (maxDiff + .5)/numSpec; + + /*cout << "LogSD:\n"; cout << "D Test Statistic = " << DStatistic << "\n"; 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] = DStatistic; - data[1] = .89196/sqrt(numSpec); + data[0] = (maxDiff + .5)/numSpec; + data[1] = 0.886/sqrt(numSpec); + data[2] = 1.031/sqrt(numSpec); + 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) { diff --git a/sharedace.cpp b/sharedace.cpp index 791b526..8730596 100644 --- a/sharedace.cpp +++ b/sharedace.cpp @@ -77,13 +77,13 @@ EstOutput SharedAce::getValues(SharedRAbundVector* shared1, SharedRAbundVector* C12 = 1 - (C12Numerator /(float) t11); part1 = S12Rare / (float)C12; part2 = 1 / (float)C12; - + //calculate gammas Gamma1 = ((S12Rare * t21) / (float)((C12 * t10 * t11)) - 1); Gamma2 = ((S12Rare * t12) / (float)((C12 * t01 * t11)) - 1); Gamma3 = ((S12Rare / C12) * (S12Rare / C12)) * ( t22 / (float)(t10 * t01 * t11)); Gamma3 = Gamma3 - ((S12Rare * t11) / (float)(C12 * t01 * t10)) - Gamma1 - Gamma2; - + if (isnan(Gamma1) || isinf(Gamma1)) { Gamma1 = 0; } if (isnan(Gamma2) || isinf(Gamma2)) { Gamma2 = 0; } if (isnan(Gamma3) || isinf(Gamma3)) { Gamma3 = 0; } diff --git a/sharedbdiversity.cpp b/sharedbdiversity.cpp deleted file mode 100644 index ae2753d..0000000 --- a/sharedbdiversity.cpp +++ /dev/null @@ -1,125 +0,0 @@ -/* - * bdiversity.cpp - * Mothur - * - * Created by Thomas Ryabin on 3/13/09. - * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. - * - */ - -#include "sharedbdiversity.h" - -/***********************************************************************/ - -double BDiversity::whitt(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ - double nz1 = (double)shared1->numNZ(); - double nz2 = (double)shared2->numNZ(); - double sum = nz1; - for(int i = 1; i < shared1->size(); i++) - { - if(shared1->get(i).abundance == 0 && shared2->get(i).abundance != 0) - sum++; - } - return 2*sum/(nz1+nz2)-1; -} - -/***********************************************************************/ - -double BDiversity::ms(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ - 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++; - } - return (b+c)/(a+b+c); -} - -/***********************************************************************/ - -double BDiversity::sor(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ - double sum = 0; - double asum = 0; - double bsum = 0; - for(int i = 1; i < shared1->size(); i++) - { - int abund1 = shared1->get(i).abundance; - int abund2 = shared2->get(i).abundance; - - asum += abund1; - bsum += abund2; - if(abund1 >= abund2) - sum += abund2; - else - sum += abund1; - } - return 2*sum/(asum+bsum); -} - -/***********************************************************************/ - -double BDiversity::mor(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ - double multSum = 0; - double powSum1 = 0; - double powSum2 = 0; - double ind1 = 0; - double ind2 = 0; - for(int i = 1; i < shared1->size(); i++) - { - double abund1 = (double)shared1->get(i).abundance; - double abund2 = (double)shared2->get(i).abundance; - multSum += abund1*abund2; - powSum1 += pow(abund1, 2); - powSum2 += pow(abund2, 2); - ind1 += abund1; - ind2 += abund2; - } - return 2*multSum / ((powSum1/pow(ind1, 2) + powSum2/pow(ind2, 2)) * (ind1*ind2)); -} -/***********************************************************************/ - -EstOutput BDiversity::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ - try { - data.resize(4,0); - - double whitt1 = whitt(shared1, shared2); - double jac1 = 1-ms(shared1, shared2); - double sor1 = sor(shared1, shared2); - double mor1 = mor(shared1, shared2); - /*cout << "Whittaker's Measure: " << whitt1 << "\n"; - cout << "Marczewski-Steinhaus distance: " << jac1 << "\n"; - cout << "Sorensen index: " << sor1 << "\n"; - cout << "Morisita-Horn index: " << mor1 << "\n";*/ - - data[0] = whitt1; - data[1] = jac1; - data[2] = sor1; - data[3] = mor1; - 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; } - if (isnan(data[3]) || isinf(data[1])) { data[3] = 0; } - - return data; - } - - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BDiversity class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} - -/***********************************************************************/ \ No newline at end of file diff --git a/sharedbdiversity.h b/sharedbdiversity.h deleted file mode 100644 index 07b4fe0..0000000 --- a/sharedbdiversity.h +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef BDIVERSITY_H -#define BDIVERSITY_H -/* - * bdiversity.h - * Mothur - * - * Created by Thomas Ryabin on 3/13/09. - * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. - * - */ -#include "calculator.h" - -/*This class implements the BDiversity estimator on 2 groups. -It is a child of the calculator class.*/ - -/***********************************************************************/ - -class BDiversity : public Calculator { - -public: - BDiversity() : Calculator("bdiversity", 3) {}; - EstOutput getValues(SAbundVector*) {return data;}; - EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*); -private: - double whitt(SharedRAbundVector*, SharedRAbundVector*); - double ms(SharedRAbundVector*, SharedRAbundVector*); - double sor(SharedRAbundVector*, SharedRAbundVector*); - double mor(SharedRAbundVector*, SharedRAbundVector*); -}; - -/***********************************************************************/ - -#endif \ No newline at end of file diff --git a/sharedchao1.cpp b/sharedchao1.cpp index 824115f..02e2411 100644 --- a/sharedchao1.cpp +++ b/sharedchao1.cpp @@ -39,6 +39,7 @@ EstOutput SharedChao1::getValues(SharedRAbundVector* sharedA, SharedRAbundVector //does B have one or two if (tempB == 1) { f1B++; } else if (tempB == 2) { f2B++; } + } } diff --git a/summarysharedcommand.cpp b/summarysharedcommand.cpp index eea7b25..e6a7d9a 100644 --- a/summarysharedcommand.cpp +++ b/summarysharedcommand.cpp @@ -21,7 +21,7 @@ #include "sharedthetayc.h" #include "sharedthetan.h" #include "sharedkstest.h" -#include "sharedbdiversity.h" +#include "whittaker.h" #include "sharedochiai.h" #include "sharedanderbergs.h" #include "sharedkulczynski.h" @@ -53,7 +53,7 @@ SummarySharedCommand::SummarySharedCommand(){ sumCalculators.push_back(new SharedAce()); }else if (globaldata->Estimators[i] == "jabund") { sumCalculators.push_back(new JAbund()); - }else if (globaldata->Estimators[i] == "sorensonabund") { + }else if (globaldata->Estimators[i] == "sorabund") { sumCalculators.push_back(new SorAbund()); }else if (globaldata->Estimators[i] == "jclass") { sumCalculators.push_back(new Jclass()); @@ -86,8 +86,8 @@ SummarySharedCommand::SummarySharedCommand(){ }else if (globaldata->Estimators[i] == "braycurtis") { sumCalculators.push_back(new BrayCurtis()); } - else if (globaldata->Estimators[i] == "bdiversity") { - sumCalculators.push_back(new BDiversity()); + else if (globaldata->Estimators[i] == "whittaker") { + sumCalculators.push_back(new Whittaker()); } } diff --git a/treegroupscommand.cpp b/treegroupscommand.cpp index 6393684..7791d5e 100644 --- a/treegroupscommand.cpp +++ b/treegroupscommand.cpp @@ -33,7 +33,7 @@ TreeGroupCommand::TreeGroupCommand(){ if (validCalculator->isValidCalculator("treegroup", globaldata->Estimators[i]) == true) { if (globaldata->Estimators[i] == "jabund") { treeCalculators.push_back(new JAbund()); - }else if (globaldata->Estimators[i] == "sorensonabund") { + }else if (globaldata->Estimators[i] == "sorabund") { treeCalculators.push_back(new SorAbund()); }else if (globaldata->Estimators[i] == "jclass") { treeCalculators.push_back(new Jclass()); diff --git a/validcalculator.cpp b/validcalculator.cpp index 25cc5c5..72d4041 100644 --- a/validcalculator.cpp +++ b/validcalculator.cpp @@ -208,7 +208,7 @@ void ValidCalculators::initialShared() { shared["sharedchao"] = "sharedchao"; shared["sharedace"] = "sharedace"; shared["jabund"] = "jabund"; - shared["sorensonabund"] = "sorensonabund"; + shared["sorabund"] = "sorabund"; shared["jclass"] = "jclass"; shared["sorclass"] = "sorclass"; shared["jest"] = "jest"; @@ -216,7 +216,7 @@ void ValidCalculators::initialShared() { shared["thetayc"] = "thetayc"; shared["thetan"] = "thetan"; shared["kstest"] = "kstest"; - shared["bdiversity"] = "bdiversity"; + shared["whittaker"] = "whittaker"; shared["sharednseqs"] = "sharednseqs"; shared["ochiai"] = "ochiai"; shared["anderberg"] = "anderberg"; @@ -276,7 +276,7 @@ void ValidCalculators::initialSummary() { summary["bergerparker"] = "bergerparker"; summary["geometric"] = "geometric"; summary["bootstrap"] = "bootstrap"; - summary["logsd"] = "logsd"; + summary["logseries"] = "logseries"; summary["qstat"] = "qstat"; summary["bstick"] = "bstick"; summary["nseqs"] = "nseqs"; @@ -300,7 +300,7 @@ void ValidCalculators::initialSharedSummary() { sharedsummary["sharedchao"] = "sharedchao"; sharedsummary["sharedace"] = "sharedace"; sharedsummary["jabund"] = "jabund"; - sharedsummary["sorensonabund"] = "sorensonabund"; + sharedsummary["sorabund"] = "sorabund"; sharedsummary["jclass"] = "jclass"; sharedsummary["sorclass"] = "sorclass"; sharedsummary["jest"] = "jest"; @@ -308,7 +308,7 @@ void ValidCalculators::initialSharedSummary() { sharedsummary["thetayc"] = "thetayc"; sharedsummary["thetan"] = "thetan"; sharedsummary["kstest"] = "kstest"; - sharedsummary["bdiversity"] = "bdiversity"; + sharedsummary["whittaker"] = "whittaker"; sharedsummary["sharednseqs"] = "sharednseqs"; sharedsummary["ochiai"] = "ochiai"; sharedsummary["anderberg"] = "anderberg"; @@ -390,7 +390,7 @@ void ValidCalculators::initialVennShared() { void ValidCalculators::initialTreeGroups() { try { treegroup["jabund"] = "jabund"; - treegroup["sorensonabund"] = "sorensonabund"; + treegroup["sorabund"] = "sorabund"; treegroup["jclass"] = "jclass"; treegroup["sorclass"] = "sorclass"; treegroup["jest"] = "jest"; @@ -412,7 +412,7 @@ void ValidCalculators::initialTreeGroups() { void ValidCalculators::initialBoot() { try { boot["jabund"] = "jabund"; - boot["sorensonabund"] = "sorensonabund"; + boot["sorabund"] = "sorabund"; boot["jclass"] = "jclass"; boot["sorclass"] = "orclass"; boot["jest"] = "jest"; diff --git a/whittaker.cpp b/whittaker.cpp new file mode 100644 index 0000000..95f1c9c --- /dev/null +++ b/whittaker.cpp @@ -0,0 +1,39 @@ +/* + * whittaker.cpp + * Mothur + * + * Created by Pat Schloss on 4/23/09. + * Copyright 2009 Patrick D. Schloss. All rights reserved. + * + */ + +#include "whittaker.h" + +/***********************************************************************/ + +EstOutput Whittaker::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ + try{ + data.resize(1); + + int countA = 0; + int countB = 0; + int sTotal = shared1->getNumBins(); + for(int i=0;igetAbundance(i) != 0){ countA++; } + if(shared2->getAbundance(i) != 0){ countB++; } + } + + data[0] = 2*sTotal/(float)(countA+countB)-1; + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Whittaker class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Whittaker class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ \ No newline at end of file diff --git a/whittaker.h b/whittaker.h new file mode 100644 index 0000000..22bf788 --- /dev/null +++ b/whittaker.h @@ -0,0 +1,29 @@ +#ifndef WHITTAKER_H +#define WHITTAKER_H +/* + * whittaker.h + * Mothur + * + * Created by Thomas Ryabin on 3/13/09. + * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. + * + */ +#include "calculator.h" + +/*This class implements the Whittaker estimator on 2 groups. + It is a child of the calculator class.*/ + +/***********************************************************************/ + +class Whittaker : public Calculator { + +public: + Whittaker() : Calculator("whittaker", 3) {}; + EstOutput getValues(SAbundVector*) {return data;}; + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*); + +}; + +/***********************************************************************/ + +#endif \ No newline at end of file