From 0e40e23448c2ee274268d85e0d0e65cb9eaeee6f Mon Sep 17 00:00:00 2001 From: westcott Date: Thu, 3 Feb 2011 16:44:19 +0000 Subject: [PATCH] changes for 1.16.0 --- Mothur.xcodeproj/project.pbxproj | 20 ++++++------- catchallcommand.cpp | 5 ++++ corraxescommand.cpp | 48 ++++++++++++++++++++++++-------- splitabundcommand.cpp | 46 ++++++++++-------------------- splitabundcommand.h | 2 +- 5 files changed, 67 insertions(+), 54 deletions(-) diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 112cf3c..c1a1dc9 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -958,6 +958,12 @@ A7E9B75C12D37EC400DA6239 /* mothur.h */, A7E9B75D12D37EC400DA6239 /* mothurout.cpp */, A7E9B75E12D37EC400DA6239 /* mothurout.h */, + A7E9B76112D37EC400DA6239 /* nast.cpp */, + A7E9B76212D37EC400DA6239 /* nast.hpp */, + A7E9B76312D37EC400DA6239 /* nastreport.cpp */, + A7E9B76412D37EC400DA6239 /* nastreport.hpp */, + A7E9B76712D37EC400DA6239 /* noalign.cpp */, + A7E9B76812D37EC400DA6239 /* noalign.hpp */, A7E9B76512D37EC400DA6239 /* needlemanoverlap.cpp */, A7E9B76612D37EC400DA6239 /* needlemanoverlap.hpp */, A7E9B77012D37EC400DA6239 /* observable.h */, @@ -974,13 +980,9 @@ A7E9B7A912D37EC400DA6239 /* rarefact.cpp */, A7E9B7AA12D37EC400DA6239 /* rarefact.h */, A7E9B7AD12D37EC400DA6239 /* rarefactioncurvedata.h */, + 7E6BE10812F710D8007ADDBE /* refchimeratest.h */, + 7E6BE10912F710D8007ADDBE /* refchimeratest.cpp */, A7E9BA5312D39A5E00DA6239 /* read */, - A7E9B76112D37EC400DA6239 /* nast.cpp */, - A7E9B76212D37EC400DA6239 /* nast.hpp */, - A7E9B76312D37EC400DA6239 /* nastreport.cpp */, - A7E9B76412D37EC400DA6239 /* nastreport.hpp */, - A7E9B76712D37EC400DA6239 /* noalign.cpp */, - A7E9B76812D37EC400DA6239 /* noalign.hpp */, A7E9B82312D37EC400DA6239 /* sharedutilities.cpp */, A7E9B82412D37EC400DA6239 /* sharedutilities.h */, A7E9B82D12D37EC400DA6239 /* singlelinkage.cpp */, @@ -1014,8 +1016,6 @@ A7E9BA3812D3956100DA6239 /* commands */ = { isa = PBXGroup; children = ( - 7E6BE10812F710D8007ADDBE /* refchimeratest.h */, - 7E6BE10912F710D8007ADDBE /* refchimeratest.cpp */, A7E9B6AE12D37EC400DA6239 /* command.hpp */, A7E9B65112D37EC300DA6239 /* aligncommand.cpp */, A7E9B65212D37EC300DA6239 /* aligncommand.h */, @@ -1598,7 +1598,7 @@ attributes = { ORGANIZATIONNAME = "Schloss Lab"; }; - buildConfigurationList = 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "Mothur" */; + buildConfigurationList = 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "mothur" */; compatibilityVersion = "Xcode 3.1"; developmentRegion = English; hasScannedForEncodings = 1; @@ -2016,7 +2016,7 @@ defaultConfigurationIsVisible = 0; defaultConfigurationName = Release; }; - 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "Mothur" */ = { + 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "mothur" */ = { isa = XCConfigurationList; buildConfigurations = ( 1DEB928A08733DD80010E9CD /* Debug */, diff --git a/catchallcommand.cpp b/catchallcommand.cpp index 8e29d62..69b9f04 100644 --- a/catchallcommand.cpp +++ b/catchallcommand.cpp @@ -121,6 +121,11 @@ CatchAllCommand::CatchAllCommand(string option) { if (sharedfile == "not open") { sharedfile = ""; abort = true; } else if (sharedfile == "not found") { sharedfile = ""; } + //check for shared file loaded during read.otu + if (sharedfile == "") { + if (globaldata->getSharedFile() != "") { sharedfile = globaldata->getSharedFile(); } + } + string label = validParameter.validFile(parameters, "label", false); if (label == "not found") { label = ""; } else { diff --git a/corraxescommand.cpp b/corraxescommand.cpp index db2e989..eac1f35 100644 --- a/corraxescommand.cpp +++ b/corraxescommand.cpp @@ -382,6 +382,8 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou int CorrAxesCommand::calcSpearman(map >& axes, ofstream& out) { try { + bool hasTies = false; + //format data vector< vector > scores; scores.resize(numaxes); for (map >::iterator it = axes.begin(); it != axes.end(); it++) { @@ -396,6 +398,7 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o for (int i = 0; i < numaxes; i++) { sort(scores[i].begin(), scores[i].end(), compareSpearman); } //find ranks of xi in each axis + vector averageX; averageX.resize(numaxes, 0.0); map > rankAxes; for (int i = 0; i < numaxes; i++) { @@ -407,6 +410,8 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o if (j != scores.size()-1) { // you are not the last so you can look ahead if (scores[i][j].score != scores[i][j+1].score) { // you are done with ties, rank them and continue + if (ties.size() > 1) { hasTies = true; } + averageX[i] += rankTotal; for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankAxes[ties[k].name].push_back(thisrank); @@ -415,12 +420,16 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o rankTotal = 0; } }else { // you are the last one + if (ties.size() > 1) { hasTies = true; } + averageX[i] += rankTotal; for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankAxes[ties[k].name].push_back(thisrank); } } } + + averageX[i] /= (float) scores[i].size(); } @@ -442,6 +451,7 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o map rankOtus; vector ties; + float averageY = 0.0; int rankTotal = 0; for (int j = 0; j < otuScores.size(); j++) { rankTotal += j; @@ -449,6 +459,8 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o if (j != scores.size()-1) { // you are not the last so you can look ahead if (otuScores[j].score != otuScores[j+1].score) { // you are done with ties, rank them and continue + if (ties.size() > 1) { hasTies = true; } + averageY += rankTotal; for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankOtus[ties[k].name] = thisrank; @@ -457,6 +469,8 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o rankTotal = 0; } }else { // you are the last one + if (ties.size() > 1) { hasTies = true; } + averageY += rankTotal; for (int k = 0; k < ties.size(); k++) { float thisrank = rankTotal / (float) ties.size(); rankOtus[ties[k].name] = thisrank; @@ -464,32 +478,44 @@ int CorrAxesCommand::calcSpearman(map >& axes, ofstream& o } } - vector pValues(numaxes); + + averageY /= (float) otuScores.size(); + //hasTies = false; + //calc spearman ranks for each axis for this otu for (int j = 0; j < numaxes; j++) { double di = 0.0; + double denom1 = 0.0; + double denom2 = 0.0; for (int k = 0; k < lookupFloat.size(); k++) { float xi = rankAxes[lookupFloat[k]->getGroup()][j]; float yi = rankOtus[lookupFloat[k]->getGroup()]; - di += ((xi - yi) * (xi - yi)); + if (hasTies) { + di += ((xi - averageX[j]) * (yi - averageY)); + denom1 += ((xi - averageX[j]) * (xi - averageX[j])); + denom2 += ((yi - averageY) * (yi - averageY)); + }else { + di += ((xi - yi) * (xi - yi)); + } } - int n = lookupFloat.size(); - double p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1))); + double p = 0.0; + + if (hasTies) { + double denom = sqrt((denom1 * denom2)); + p = di / denom; + }else { + int n = lookupFloat.size(); + p = 1.0 - ((6 * di) / (float) (n * ((n*n) - 1))); + } out << '\t' << p; - pValues[j] = p; - } - double sum = 0; - for(int k=0;khasPath(listfile); } - - //remove old files so you can append later.... - string fileroot = outputDir + m->getRootName(m->getSimpleName(listfile)); - if (Groups.size() == 0) { - remove((fileroot + "rare.list").c_str()); - remove((fileroot + "abund.list").c_str()); - - outputNames.push_back((fileroot + "rare.list")); - outputNames.push_back((fileroot + "abund.list")); - outputTypes["list"].push_back((fileroot + "rare.list")); - outputTypes["list"].push_back((fileroot + "abund.list")); - }else{ - for (int i=0; i processedLabels; @@ -390,9 +368,11 @@ int SplitAbundCommand::splitList(ListVector* thisList) { } }//end for - writeList(thisList); string tag = thisList->getLabel() + "."; + + writeList(thisList, tag); + if (groupfile != "") { parseGroup(tag); } if (accnos) { writeAccnos(tag); } if (fastafile != "") { parseFasta(tag); } @@ -406,7 +386,7 @@ int SplitAbundCommand::splitList(ListVector* thisList) { } } /**********************************************************************************************************************/ -int SplitAbundCommand::writeList(ListVector* thisList) { +int SplitAbundCommand::writeList(ListVector* thisList, string tag) { try { map filehandles; @@ -428,13 +408,13 @@ int SplitAbundCommand::writeList(ListVector* thisList) { ofstream aout; ofstream rout; - string rare = outputDir + m->getRootName(m->getSimpleName(listfile)) + "rare.list"; - m->openOutputFileAppend(rare, rout); - //outputNames.push_back(rare); + string rare = outputDir + m->getRootName(m->getSimpleName(listfile)) + tag + "rare.list"; + m->openOutputFile(rare, rout); + outputNames.push_back(rare); outputTypes["list"].push_back(rare); - string abund = outputDir + m->getRootName(m->getSimpleName(listfile)) + "abund.list"; - m->openOutputFileAppend(abund, aout); - //outputNames.push_back(abund); + string abund = outputDir + m->getRootName(m->getSimpleName(listfile)) + tag + "abund.list"; + m->openOutputFile(abund, aout); + outputNames.push_back(abund); outputTypes["list"].push_back(abund); if (rareNames.size() != 0) { rout << thisList->getLabel() << '\t' << numRareBins << '\t'; } if (abundNames.size() != 0) { aout << thisList->getLabel() << '\t' << numAbundBins << '\t'; } @@ -470,8 +450,10 @@ int SplitAbundCommand::writeList(ListVector* thisList) { temp2 = new ofstream; filehandles[Groups[i]+".abund"] = temp2; - m->openOutputFileAppend(fileroot + Groups[i] + ".rare.list", *(filehandles[Groups[i]+".rare"])); - m->openOutputFileAppend(fileroot + Groups[i] + ".abund.list", *(filehandles[Groups[i]+".abund"])); + m->openOutputFile(fileroot + Groups[i] + tag + ".rare.list", *(filehandles[Groups[i]+".rare"])); + m->openOutputFile(fileroot + Groups[i] + tag + ".abund.list", *(filehandles[Groups[i]+".abund"])); + outputNames.push_back(fileroot + Groups[i] + tag + ".rare.list"); outputTypes["list"].push_back(fileroot + Groups[i] + tag + ".rare.list"); + outputNames.push_back(fileroot + Groups[i] + tag + ".abund.list"); outputTypes["list"].push_back(fileroot + Groups[i] + tag + ".abund.list"); } map groupVector; diff --git a/splitabundcommand.h b/splitabundcommand.h index 02d61f3..2e409b2 100644 --- a/splitabundcommand.h +++ b/splitabundcommand.h @@ -42,7 +42,7 @@ private: int splitList(ListVector*); int splitNames(); //namefile int writeNames(); - int writeList(ListVector*); + int writeList(ListVector*, string); int writeAccnos(string); int parseGroup(string); int parseFasta(string); -- 2.39.2