From e946299389ed49b4afc4262a4e8070ee3d804769 Mon Sep 17 00:00:00 2001 From: SarahsWork Date: Wed, 6 Mar 2013 08:54:53 -0500 Subject: [PATCH] <= change for anova and homova. finished chimera.uchime dereplicate=t issue. --- amovacommand.cpp | 2 +- chimerauchimecommand.cpp | 47 ++++++++++++++++++++-------------------- counttable.cpp | 26 +++++++++++++++++----- homovacommand.cpp | 2 +- 4 files changed, 46 insertions(+), 31 deletions(-) diff --git a/amovacommand.cpp b/amovacommand.cpp index a78aafc..492a096 100644 --- a/amovacommand.cpp +++ b/amovacommand.cpp @@ -307,7 +307,7 @@ double AmovaCommand::runAMOVA(ofstream& AMOVAFile, map > gro for(int i=0;i > randomizedGroup = getRandomizedGroups(groupSampleMap); double ssWithinRand = calcSSWithin(randomizedGroup); - if(ssWithinRand < ssWithinOrig){ counter++; } + if(ssWithinRand <= ssWithinOrig){ counter++; } } double pValue = (double)counter / (double) iters; diff --git a/chimerauchimecommand.cpp b/chimerauchimecommand.cpp index 3ed6170..cb155eb 100644 --- a/chimerauchimecommand.cpp +++ b/chimerauchimecommand.cpp @@ -727,10 +727,8 @@ int ChimeraUchimeCommand::execute(){ if(processors == 1) { totalSeqs = driverGroups(outputFileName, newFasta, accnosFileName, alnsFileName, newCountFile, 0, groups.size(), groups); - //read my own - if (hasCount && !dups) { - CountTable newCount; newCount.readTable(nameFile); - + if (hasCount && dups) { + CountTable c; c.readTable(nameFile); if (!m->isBlank(newCountFile)) { ifstream in2; m->openInputFile(newCountFile, in2); @@ -738,12 +736,12 @@ int ChimeraUchimeCommand::execute(){ string name, group; while (!in2.eof()) { in2 >> name >> group; m->gobble(in2); - newCount.setAbund(name, group, 0); + c.setAbund(name, group, 0); } in2.close(); } m->mothurRemove(newCountFile); - newCount.printTable(newCountFile); + c.printTable(newCountFile); } }else { totalSeqs = createProcessesGroups(outputFileName, newFasta, accnosFileName, alnsFileName, newCountFile, groups, nameFile, groupFile, fastaFileNames[s]); } @@ -757,25 +755,26 @@ int ChimeraUchimeCommand::execute(){ m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(totalSeqs) + " sequences. " + toString(totalChimeras) + " chimeras were found."); m->mothurOutEndLine(); m->mothurOut("The number of sequences checked may be larger than the number of unique sequences because some sequences are found in several samples."); m->mothurOutEndLine(); }else { - /*if (hasCount) { //removed empty seqs - ofstream out2; - m->openOutputFile(accnosFileName, out2); - + + if (hasCount) { + set doNotRemove; CountTable c; c.readTable(newCountFile); - vector nseqs = c.getNamesOfSeqs(); - vector ngroups = c.getNamesOfGroups(); - for (int l = 0; l < nseqs.size(); l++) { - if (c.getNumSeqs(nseqs[l]) == 0) { - c.remove(nseqs[l]); - out2 << nseqs[l] << endl; - } + vector namesInTable = c.getNamesOfSeqs(); + for (int i = 0; i < namesInTable.size(); i++) { + int temp = c.getNumSeqs(namesInTable[i]); + if (temp == 0) { c.remove(namesInTable[i]); } + else { doNotRemove.insert((namesInTable[i])); } } - for (int l = 0; l < ngroups.size(); l++) { - if (c.getGroupCount(ngroups[l]) == 0) { c.removeGroup(ngroups[l]); } + //remove names we want to keep from accnos file. + set accnosNames = m->readAccnos(accnosFileName); + ofstream out2; + m->openOutputFile(accnosFileName, out2); + for (set::iterator it = accnosNames.begin(); it != accnosNames.end(); it++) { + if (doNotRemove.count(*it) == 0) { out2 << (*it) << endl; } } out2.close(); c.printTable(newCountFile); - }*/ + } } if (hasCount) { delete cparser; } @@ -1893,7 +1892,7 @@ int ChimeraUchimeCommand::createProcessesGroups(string outputFName, string filen #endif - + //read my own if (hasCount && dups) { if (!m->isBlank(accnos + ".byCount")) { @@ -1909,7 +1908,7 @@ int ChimeraUchimeCommand::createProcessesGroups(string outputFName, string filen } m->mothurRemove(accnos + ".byCount"); } - + //append output files for(int i=0;iappendFiles((outputFName + toString(processIDS[i]) + ".temp"), outputFName); @@ -1941,8 +1940,8 @@ int ChimeraUchimeCommand::createProcessesGroups(string outputFName, string filen } //print new *.pick.count_table - if (hasCount && dups) { newCount.printTable(newCountFile); } - + if (hasCount && dups) { newCount.printTable(newCountFile); } + return num; } diff --git a/counttable.cpp b/counttable.cpp index ad0b2da..c16bf92 100644 --- a/counttable.cpp +++ b/counttable.cpp @@ -283,7 +283,23 @@ int CountTable::printTable(string file) { for (int i = 0; i < groups.size(); i++) { out << groups[i] << '\t'; } out << endl; - for (map::iterator itNames = indexNameMap.begin(); itNames != indexNameMap.end(); itNames++) { + map reverse; //use this to preserve order + for (map::iterator it = indexNameMap.begin(); it !=indexNameMap.end(); it++) { reverse[it->second] = it->first; } + + for (int i = 0; i < totals.size(); i++) { + map::iterator itR = reverse.find(i); + + if (itR != reverse.end()) { //will equal end if seqs were removed because remove just removes from indexNameMap + out << itR->second << '\t' << totals[i] << '\t'; + if (hasGroups) { + for (int j = 0; j < groups.size(); j++) { + out << counts[i][j] << '\t'; + } + } + out << endl; + } + } + /*for (map::iterator itNames = indexNameMap.begin(); itNames != indexNameMap.end(); itNames++) { out << itNames->first << '\t' << totals[itNames->second] << '\t'; if (hasGroups) { @@ -292,7 +308,7 @@ int CountTable::printTable(string file) { } } out << endl; - } + }*/ out.close(); return 0; } @@ -364,7 +380,7 @@ int CountTable::getGroupCount(string groupName) { if (hasGroups) { map::iterator it = indexGroupMap.find(groupName); if (it == indexGroupMap.end()) { - m->mothurOut("[ERROR]: " + groupName + " is not in your count table. Please correct.\n"); m->control_pressed = true; + m->mothurOut("[ERROR]: group " + groupName + " is not in your count table. Please correct.\n"); m->control_pressed = true; }else { return totalGroups[it->second]; } @@ -384,11 +400,11 @@ int CountTable::getGroupCount(string seqName, string groupName) { if (hasGroups) { map::iterator it = indexGroupMap.find(groupName); if (it == indexGroupMap.end()) { - m->mothurOut("[ERROR]: " + groupName + " is not in your count table. Please correct.\n"); m->control_pressed = true; + m->mothurOut("[ERROR]: group " + groupName + " is not in your count table. Please correct.\n"); m->control_pressed = true; }else { map::iterator it2 = indexNameMap.find(seqName); if (it2 == indexNameMap.end()) { - m->mothurOut("[ERROR]: " + seqName + " is not in your count table. Please correct.\n"); m->control_pressed = true; + m->mothurOut("[ERROR]: seq " + seqName + " is not in your count table. Please correct.\n"); m->control_pressed = true; }else { return counts[it2->second][it->second]; } diff --git a/homovacommand.cpp b/homovacommand.cpp index 49b8b64..488b68e 100644 --- a/homovacommand.cpp +++ b/homovacommand.cpp @@ -308,7 +308,7 @@ double HomovaCommand::runHOMOVA(ofstream& HOMOVAFile, map > vector ssWithinRandVector; map > randomizedGroup = getRandomizedGroups(groupSampleMap); double bValueRand = calcBValue(randomizedGroup, ssWithinRandVector); - if(bValueRand > bValueOrig){ counter++; } + if(bValueRand >= bValueOrig){ counter++; } } double pValue = (double) counter / (double) iters; -- 2.39.2