A7825503116519F70002E2DD /* chimerabellerophoncommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = chimerabellerophoncommand.cpp; sourceTree = "<group>"; };
A78434881162224F00100BE0 /* chimeraccodecommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = chimeraccodecommand.h; sourceTree = "<group>"; };
A78434891162224F00100BE0 /* chimeraccodecommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = chimeraccodecommand.cpp; sourceTree = "<group>"; };
+ A7AE6302121C3408001DE6FC /* sharedrabundfloatvector.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedrabundfloatvector.h; sourceTree = "<group>"; };
+ A7AE6303121C3408001DE6FC /* sharedrabundfloatvector.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedrabundfloatvector.cpp; sourceTree = "<group>"; };
A7BBDA7B11B5694E006E6551 /* classifyotucommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = classifyotucommand.h; sourceTree = "<group>"; };
A7BBDA7C11B5694E006E6551 /* classifyotucommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = classifyotucommand.cpp; sourceTree = "<group>"; };
A7D215C811996C6E00F13F13 /* clearcutcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = clearcutcommand.h; sourceTree = "<group>"; };
A7DA212B113FECD400BF472F /* sharedordervector.h */,
A7DA212C113FECD400BF472F /* sharedrabundvector.cpp */,
A7DA212D113FECD400BF472F /* sharedrabundvector.h */,
- A7DA212E113FECD400BF472F /* sharedsabundvector.cpp */,
+ A7AE6302121C3408001DE6FC /* sharedrabundfloatvector.h */,
+ A7AE6303121C3408001DE6FC /* sharedrabundfloatvector.cpp */,
A7DA212F113FECD400BF472F /* sharedsabundvector.h */,
+ A7DA212E113FECD400BF472F /* sharedsabundvector.cpp */,
A7DA214C113FECD400BF472F /* sparsematrix.cpp */,
A7DA214D113FECD400BF472F /* sparsematrix.hpp */,
A7DA214E113FECD400BF472F /* suffixdb.cpp */,
data.resize(3,0);
double ace, acelci, acehci;
- int nrare = 0;
- int srare = 0;
- int sabund = 0;
+ double nrare = 0;
+ double srare = 0;
+ double sabund = 0;
double Cace, term1, gamace;
- int numsum = 0;
+ double numsum = 0;
double maxRank = (double)rank->getMaxRank();
}
else if(i>abund) {sabund += rank->get(i);}
}
- int sobs = srare + sabund;
+ double sobs = srare + sabund;
if (nrare == 0){ Cace = 0.0000; }
else { Cace = 1.0000 -(double)rank->get(1)/(double)nrare; }
I have also added the forumlae to calculate the 95% confidence intervals.
*/
- int j,D_s=0,nn=0,ww=0,Max_Index=rank->getMaxRank()+1;
+ double j,D_s=0,nn=0,ww=0;
+ int Max_Index=rank->getMaxRank()+1;
double pp, temp1, temp2;
vector<double> Part_N_Part_F(Max_Index+1,0.0);
int ableToOpen;
ifstream in;
-
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
ableToOpen = openInputFile(candidateFileNames[i], in, "noerror");
}
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + candidateFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
int tag = 2001;
vector<unsigned long int> MPIPos;
MPIWroteAccnos = false;
-
+
MPI_Status status;
MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
MPI_Comm_size(MPI_COMM_WORLD, &processors);
//if there is a possibility that this sequence should be reversed
if (candidateSeq->getNumBases() < numBasesNeeded) {
- string wasBetter = "";
+ string wasBetter = "";
//if the user wants you to try the reverse
if (flip) {
//get reverse compliment
delete nast;
nast = nast2;
needToDeleteCopy = true;
+ wasBetter = "\treverse complement produced a better alignment, so mothur used the reverse complement.";
}else{
- wasBetter = "\treverse complement did NOT produce a better alignment, please check sequence.";
+ wasBetter = "\treverse complement did NOT produce a better alignment so it was not used, please check sequence.";
delete nast2;
delete copy;
}
int length = MPIPos[start+i+1] - MPIPos[start+i];
char* buf4 = new char[length];
- memcpy(buf4, outputString.c_str(), length);
+ //memcpy(buf4, outputString.c_str(), length);
MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
delete buf4;
if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
- istringstream iss (tempBuf,istringstream::in);
+ istringstream iss (tempBuf,istringstream::in);
+
Sequence* candidateSeq = new Sequence(iss);
+
int origNumBases = candidateSeq->getNumBases();
string originalUnaligned = candidateSeq->getUnaligned();
int numBasesNeeded = origNumBases * threshold;
Alignment::Alignment(int A) : nCols(A), nRows(A) {
try {
+
m = MothurOut::getInstance();
alignment.resize(nRows); // For the Gotoh and Needleman-Wunsch we initialize the dynamic programming
for(int i=0;i<nRows;i++){ // matrix by initializing a matrix that is A x A. By default we will set A
ifstream kmerFileTest(kmerDBName.c_str());
search->readKmerDB(kmerFileTest);
}
-
+
search->setNumSeqs(numSeqs);
}
}
//vector<double> bootData(3,0);
data.resize(1,0);
double maxRank = (double)rank->getMaxRank();
- int sampled = rank->getNumSeqs();
- int sobs = rank->getNumBins();
+ double sampled = rank->getNumSeqs();
+ double sobs = rank->getNumBins();
double boot = (double)sobs;
sum += 1/(double)i;
return sum;
}
-
+/***********************************************************************/
RAbundVector BStick::getRAbundVector(SAbundVector* rank){
vector <int> rData;
int mr = 1;
int ns = 0;
for(int i = rank->size()-1; i > 0; i--) {
- int cur = rank->get(i);
+ double cur = rank->get(i);
if(mr == 1 && cur > 0)
mr = i;
nb += cur;
int ableToOpen;
ifstream in;
-
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
ableToOpen = openInputFile(fastaFileNames[i], in, "noerror");
}
}
in.close();
-
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
int ableToOpen;
ifstream in;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ableToOpen = openInputFile(fastaFileNames[i], in, "noerror");
//if you can't open it, try default location
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
//erase from file list
int ableToOpen;
ifstream in;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ableToOpen = openInputFile(fastaFileNames[i], in, "noerror");
//if you can't open it, try default location
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + fastaFileNames[i] +". It will be disregarded."); m->mothurOutEndLine();
//erase from file list
int ableToOpen;
ifstream in;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ableToOpen = openInputFile(nameFileNames[i], in, "noerror");
//if you can't open it, try default location
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + nameFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
//erase from file list
int ableToOpen;
ifstream in;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ableToOpen = openInputFile(fastaFileNames[i], in, "noerror");
//if you can't open it, try default location
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
//erase from file list
int ableToOpen;
ifstream in;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ableToOpen = openInputFile(fastaFileNames[i], in, "noerror");
//if you can't open it, try default location
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
//erase from file list
int ableToOpen;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ifstream in;
ableToOpen = openInputFile(fastaFileNames[i], in, "noerror");
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
//erase from file list
}
int ableToOpen;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ifstream in;
ableToOpen = openInputFile(namefileNames[i], in, "noerror");
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + namefileNames[i] + ". It will be disregarded."); m->mothurOutEndLine(); abort = true;
//erase from file list
}
int ableToOpen;
- #ifdef USE_MPI
- int pid;
- MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
- MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
- if (pid == 0) {
- #endif
-
ifstream in;
ableToOpen = openInputFile(groupfileNames[i], in, "noerror");
}
in.close();
- #ifdef USE_MPI
- for (int j = 1; j < processors; j++) {
- MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
- }
- }else{
- MPI_Status status;
- MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
- }
-
- #endif
-
if (ableToOpen == 1) {
m->mothurOut("Unable to open " + groupfileNames[i] + ". It will be disregarded."); m->mothurOutEndLine(); groupfileNames[i] = "";
//erase from file list
temp = validParameter.validFile(parameters, "shuffle", false); if (temp == "not found"){ temp = "F"; }
shuffle = isTrue(temp);
- temp = validParameter.validFile(parameters, "neighbor", false); if (temp == "not found"){ temp = "F"; }
+ temp = validParameter.validFile(parameters, "neighbor", false); if (temp == "not found"){ temp = "T"; }
neighbor = isTrue(temp);
temp = validParameter.validFile(parameters, "DNA", false); if (temp == "not found"){ temp = "F"; }
m->mothurOut("The seed parameter allows you to explicitly set the PRNG seed to a specific value. \n");
m->mothurOut("The norandom parameter allows you to attempt joins deterministically, default=F. \n");
m->mothurOut("The shuffle parameter allows you to randomly shuffle the distance matrix, default=F. \n");
- m->mothurOut("The neighbor parameter allows you to use traditional Neighbor-Joining algorithm, default=F. \n");
+ m->mothurOut("The neighbor parameter allows you to use traditional Neighbor-Joining algorithm, default=T. \n");
m->mothurOut("The DNA parameter allows you to indicate your fasta file contains DNA sequences, default=F. \n");
m->mothurOut("The protein parameter allows you to indicate your fasta file contains protein sequences, default=F. \n");
rabundFile.close();
listFile.close();
- if (saveCutoff != cutoff) { m->mothurOut("changed cutoff to " + toString(cutoff)); m->mothurOutEndLine(); }
+ if (saveCutoff != cutoff) {
+ if (hard) { saveCutoff = ceilDist(saveCutoff, precision); }
+ else { saveCutoff = roundDist(saveCutoff, precision); }
+
+ m->mothurOut("changed cutoff to " + toString(cutoff)); m->mothurOutEndLine();
+ }
m->mothurOutEndLine();
m->mothurOut("Output File Names: "); m->mothurOutEndLine();
remove(thisDistFile.c_str());
remove(thisNamefile.c_str());
- if (saveCutoff != cutoff) { m->mothurOut("Cutoff was " + toString(cutoff) + " changed cutoff to " + toString(saveCutoff)); m->mothurOutEndLine(); }
+ if (saveCutoff != cutoff) {
+ if (hard) { saveCutoff = ceilDist(saveCutoff, precision); }
+ else { saveCutoff = roundDist(saveCutoff, precision); }
+
+ m->mothurOut("Cutoff was " + toString(cutoff) + " changed cutoff to " + toString(saveCutoff)); m->mothurOutEndLine();
+ }
if (saveCutoff < smallestCutoff) { smallestCutoff = saveCutoff; }
}
MPI_Status status;
MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
- cout << pid << "is in create filter " << endl;
+
MPI_File outMPI;
MPI_File tempMPI;
MPI_File inMPI;
out << '>' << seq.getName() << endl << filterSeq << endl;
}
gobble(in);
+
+ //report progress
+ if((i+1) % 100 == 0){ m->mothurOut(toString(i+1)); m->mothurOutEndLine(); }
}
+
+ //report progress
+ if((line->num) % 100 != 0){ m->mothurOut(toString(line->num)); m->mothurOutEndLine(); }
+
out.close();
in.close();
//loop through and create all the processes you want
while (process != processors) {
- int pid = vfork();
+ int pid = fork();
if (pid > 0) {
processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
process++;
}else if (pid == 0){
driverCreateFilter(F, filename, lines[process]);
+
+ //write out filter counts to file
+ filename += toString(getpid()) + "filterValues.temp";
+ ofstream out;
+ openOutputFile(filename, out);
+
+ for (int k = 0; k < alignmentLength; k++) { out << F.a[k] << '\t'; } out << endl;
+ for (int k = 0; k < alignmentLength; k++) { out << F.t[k] << '\t'; } out << endl;
+ for (int k = 0; k < alignmentLength; k++) { out << F.g[k] << '\t'; } out << endl;
+ for (int k = 0; k < alignmentLength; k++) { out << F.c[k] << '\t'; } out << endl;
+ for (int k = 0; k < alignmentLength; k++) { out << F.gap[k] << '\t'; } out << endl;
+
+ out.close();
+
exit(0);
}else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
}
wait(&temp);
}
+ //parent reads in and combine Filter info
+ for (int i = 0; i < processIDS.size(); i++) {
+ string tempFilename = filename + toString(processIDS[i]) + "filterValues.temp";
+ ifstream in;
+ openInputFile(tempFilename, in);
+
+ int temp;
+ for (int k = 0; k < alignmentLength; k++) { in >> temp; F.a[k] += temp; } gobble(in);
+ for (int k = 0; k < alignmentLength; k++) { in >> temp; F.t[k] += temp; } gobble(in);
+ for (int k = 0; k < alignmentLength; k++) { in >> temp; F.g[k] += temp; } gobble(in);
+ for (int k = 0; k < alignmentLength; k++) { in >> temp; F.c[k] += temp; } gobble(in);
+ for (int k = 0; k < alignmentLength; k++) { in >> temp; F.gap[k] += temp; } gobble(in);
+
+ in.close();
+ remove(tempFilename.c_str());
+ }
+
return exitCommand;
#endif
}
data.resize(3,0);
rdata = getRAbundVector(rank);
- int numInd = rdata.getNumSeqs();
- int numSpec = rdata.getNumBins();
- int min = rdata.get(rdata.size()-1);
+ double numInd = rdata.getNumSeqs();
+ double numSpec = rdata.getNumBins();
+ double min = rdata.get(rdata.size()-1);
double k = .5;
double step = .49999;
#include "heatmap.h"
//**********************************************************************************************************************
-HeatMap::HeatMap(string sort, string scale, string dir){
+HeatMap::HeatMap(string sort, string scale, int num, int fsize, string dir){
try {
globaldata = GlobalData::getInstance();
m = MothurOut::getInstance();
sorted = sort;
scaler = scale;
outputDir = dir;
+ numOTU = num;
+ fontSize = fsize;
}
catch(exception& e) {
m->errorOut(e, "HeatMap", "HeatMap");
//white backround
outsvg << "<rect fill=\"white\" stroke=\"white\" x=\"0\" y=\"0\" width=\"" + toString(300) + "\" height=\"" + toString((rabund->getNumBins()*5 + 120)) + "\"/>";
- outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((150) - 40) + "\" y=\"25\">Heatmap at distance " + rabund->getLabel() + "</text>\n";
+ outsvg << "<text fill=\"black\" class=\"seri\" font-size=\"" + toString(fontSize) + "\" x=\"" + toString((150) - 40) + "\" y=\"25\">Heatmap at distance " + rabund->getLabel() + "</text>\n";
//output legend and color labels
string color;
string HeatMap::getPic(vector<SharedRAbundVector*> lookup) {
try {
+
+ int numBinsToDisplay = lookup[0]->size();
+
+ if (numOTU != 0) { //user want to display a portion of the otus
+ if (numOTU < numBinsToDisplay) { numBinsToDisplay = numOTU; }
+ }
+
//sort lookup so shared bins are on top
- if (isTrue(sorted) == true) { sortSharedVectors(lookup); }
+ if (sorted != "none") { sortSharedVectors(lookup); }
vector<vector<string> > scaleRelAbund;
vector<float> maxRelAbund(lookup.size(), 0.0);
//white backround
outsvg << "<rect fill=\"white\" stroke=\"white\" x=\"0\" y=\"0\" width=\"" + toString(lookup.size() * 300) + "\" height=\"" + toString((lookup[0]->getNumBins()*5 + 120)) + "\"/>";
- outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((lookup.size() * 150) - 40) + "\" y=\"25\">Heatmap at distance " + lookup[0]->getLabel() + "</text>\n";
+ outsvg << "<text fill=\"black\" class=\"seri\" font-size=\"" + toString(fontSize) + "\" text-anchor=\"middle\" x=\"" + toString((lookup.size() * 150) - 40) + "\" y=\"25\">Heatmap at distance " + lookup[0]->getLabel() + "</text>\n";
//column labels
for (int h = 0; h < lookup.size(); h++) {
- outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString(((300 * (h+1)) - 150) - ((int)lookup[h]->getGroup().length() / 2)) + "\" y=\"50\">" + lookup[h]->getGroup() + "</text>\n";
+ outsvg << "<text fill=\"black\" class=\"seri\" font-size=\"" + toString(fontSize) + "\" x=\"" + toString(((300 * (h+1)) - 150) - ((int)lookup[h]->getGroup().length() / 2)) + "\" y=\"50\">" + lookup[h]->getGroup() + "</text>\n";
}
-
+
//output legend and color labels
string color;
int x = 0;
printLegend(y, superMaxRelAbund);
y = 70;
- for (int i = 0; i < scaleRelAbund[0].size(); i++) {
+ for (int i = 0; i < numBinsToDisplay; i++) {
for (int j = 0; j < scaleRelAbund.size(); j++) {
if (m->control_pressed) { outsvg.close(); return "control"; }
}
//**********************************************************************************************************************
-void HeatMap::sortSharedVectors(vector<SharedRAbundVector*>& lookup){
+int HeatMap::sortSharedVectors(vector<SharedRAbundVector*>& lookup){
try {
- //copy lookup and then clear it to refill with sorted.
- //loop though lookup and determine if they are shared
- //if they are then insert in the front
- //if not push to back
-
+
vector<SharedRAbundVector*> looktemp;
map<int, int> place; //spot in lookup where you insert shared by, ie, 3 -> 2 if they are shared by 3 inset into location 2.
map<int, int>::iterator it;
- int count;
+ /****************** find order of otus **********************/
+ if (sorted == "shared") {
+ place = orderShared(lookup);
+ }else if (sorted == "topotu") {
+ place = orderTopOtu(lookup);
+ }else if (sorted == "topgroup") {
+ place = orderTopGroup(lookup);
+ }else { m->mothurOut("Error: invalid sort option."); m->mothurOutEndLine(); return 1; }
+
+
+ /******************* create copy of lookup *********************/
//create and initialize looktemp as a copy of lookup
for (int i = 0; i < lookup.size(); i++) {
SharedRAbundVector* temp = new SharedRAbundVector(lookup[i]->getNumBins());
}
looktemp.push_back(temp);
}
+
+ /************************ fill lookup in order given by place *********************/
+ //for each bin
+ for (int i = 0; i < looktemp[0]->size(); i++) { //place
+ //fill lookup // 2 -> 1
+ for (int j = 0; j < looktemp.size(); j++) { // 3 -> 2
+ int newAbund = looktemp[j]->getAbundance(i); // 1 -> 3
+ lookup[j]->set(place[i], newAbund, looktemp[j]->getGroup()); //binNumber, abundance, group
+ }
+ }
- //clear out lookup to create sorted lookup -- Sarah look at - this is causing segmentation faults
- for (int j = 0; j < lookup.size(); j++) {
-// delete lookup[j];
+ //delete looktemp -- Sarah look at - this is causing segmentation faults
+ for (int j = 0; j < looktemp.size(); j++) {
+// delete looktemp[j];
}
- lookup.clear(); //doesn't this do the job?
-
- //create and initialize lookup to empty vectors
- for (int i = 0; i < looktemp.size(); i++) {
- SharedRAbundVector* temp = new SharedRAbundVector();
- temp->setLabel(looktemp[i]->getLabel());
- temp->setGroup(looktemp[i]->getGroup());
- lookup.push_back(temp);
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "HeatMap", "sortSharedVectors");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+map<int, int> HeatMap::orderShared(vector<SharedRAbundVector*>& lookup){
+ try {
+
+ map<int, int> place; //spot in lookup where you insert shared by, ie, 3 -> 2 if they are shared by 3 inset into location 2.
+ map<int, int>::iterator it;
+
+ vector<int> sharedBins;
+ vector<int> uniqueBins;
+
+ //for each bin
+ for (int i = 0; i < lookup[0]->size(); i++) {
+ int count = 0;
+
+ //is this bin shared
+ for (int j = 0; j < lookup.size(); j++) { if (lookup[j]->getAbundance(i) != 0) { count++; } }
- //initialize place map
- place[i] = 0;
+ if (count < 2) { uniqueBins.push_back(i); }
+ else { sharedBins.push_back(i); }
}
+ //fill place
+ for (int i = 0; i < sharedBins.size(); i++) { place[sharedBins[i]] = i; }
+ for (int i = 0; i < uniqueBins.size(); i++) { place[uniqueBins[i]] = (sharedBins.size() + i); }
+
+ return place;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "HeatMap", "orderShared");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+map<int, int> HeatMap::orderTopOtu(vector<SharedRAbundVector*>& lookup){
+ try {
+
+ map<int, int> place; //spot in lookup where you insert shared by, ie, 3 -> 2 if they are shared by 3 inset into location 2.
+ map<int, int>::iterator it;
+
+ vector<binCount> totals;
//for each bin
- for (int i = 0; i < looktemp[0]->size(); i++) {
- count = 0;
- bool updatePlace = false;
- //for each group
- for (int j = 0; j < looktemp.size(); j++) {
- if (looktemp[j]->getAbundance(i) != 0) { count++; }
- }
+ for (int i = 0; i < lookup[0]->size(); i++) {
+ int total = 0;
- //fill lookup
- for (int j = 0; j < looktemp.size(); j++) {
- //if they are not shared then push to back, if they are not insert in front
- if (count < 2) { lookup[j]->push_back(looktemp[j]->getAbundance(i), looktemp[j]->getGroup()); }
- //they are shared by some
- else { lookup[j]->insert(looktemp[j]->getAbundance(i), place[count], looktemp[j]->getGroup()); updatePlace = true; }
- }
+ for (int j = 0; j < lookup.size(); j++) { total += lookup[j]->getAbundance(i); }
- if (updatePlace == true) {
- //move place holders below where you entered up to "make space" for you entry
- for (it = place.begin(); it!= place.end(); it++) {
- if (it->first < count) { it->second++; }
- }
+ binCount temp(i, total);
+
+ totals.push_back(temp);
+ }
+
+ sort(totals.begin(), totals.end(), comparebinCounts);
+
+ //fill place
+ for (int i = 0; i < totals.size(); i++) { place[totals[i].bin] = i; }
+
+ return place;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "HeatMap", "orderTopOtu");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+map<int, int> HeatMap::orderTopGroup(vector<SharedRAbundVector*>& lookup){
+ try {
+
+ map<int, int> place; //spot in lookup where you insert shared by, ie, 3 -> 2 if they are shared by 3 inset into location 2.
+ map<int, int>::iterator it;
+
+ vector < vector<binCount> > totals; //totals[0] = bin totals for group 0, totals[1] = bin totals for group 1, ...
+ totals.resize(lookup.size());
+
+ //for each bin
+ for (int i = 0; i < lookup[0]->size(); i++) {
+ for (int j = 0; j < lookup.size(); j++) {
+ binCount temp(i, (lookup[j]->getAbundance(i)));
+ totals[j].push_back(temp);
}
}
- //delete looktemp -- Sarah look at - this is causing segmentation faults
- for (int j = 0; j < looktemp.size(); j++) {
-// delete looktemp[j];
+ for (int i = 0; i < totals.size(); i++) { sort(totals[i].begin(), totals[i].end(), comparebinCounts); }
+
+ //fill place
+ //grab the top otu for each group adding it if its not already added
+ int count = 0;
+ for (int i = 0; i < totals[0].size(); i++) {
+
+ for (int j = 0; j < totals.size(); j++) {
+ it = place.find(totals[j][i].bin);
+
+ if (it == place.end()) { //not added yet
+ place[totals[j][i].bin] = count;
+ count++;
+ }
+ }
}
+
+ return place;
}
catch(exception& e) {
- m->errorOut(e, "HeatMap", "sortSharedVectors");
+ m->errorOut(e, "HeatMap", "orderTopGroup");
exit(1);
}
}
-
//**********************************************************************************************************************
void HeatMap::printLegend(int y, float maxbin) {
exit(1);
}
}
+//**********************************************************************************************************************
+
#include "datavector.hpp"
#include "globaldata.hpp"
+/***********************************************************************/
+struct binCount {
+ int bin;
+ int abund;
+ binCount(int i, int j) : bin(i), abund(j) {}
+};
+/***********************************************************************/
+//sorts highest abund to lowest
+inline bool comparebinCounts(binCount left, binCount right){
+ return (left.abund > right.abund);
+}
/***********************************************************************/
class HeatMap {
public:
- HeatMap(string, string, string);
+ HeatMap(string, string, int, int, string);
~HeatMap(){};
string getPic(RAbundVector*);
string getPic(vector<SharedRAbundVector*>);
private:
- void sortSharedVectors(vector<SharedRAbundVector*>& );
+ int sortSharedVectors(vector<SharedRAbundVector*>& );
void printLegend(int, float);
GlobalData* globaldata;
string format, sorted, groupComb, scaler, outputDir;
ofstream outsvg;
MothurOut* m;
+ int numOTU, fontSize;
+
+ map<int, int> orderTopGroup(vector<SharedRAbundVector*>&);
+ map<int, int> orderTopOtu(vector<SharedRAbundVector*>&);
+ map<int, int> orderShared(vector<SharedRAbundVector*>&);
};
else {
//valid paramters for this command
- string AlignArray[] = {"groups","label","sorted","scale","outputdir","inputdir"};
+ string AlignArray[] = {"groups","label","sorted","scale","fontsize","numotu","outputdir","inputdir"};
vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
OptionParser parser(option);
globaldata->Groups = Groups;
}
- sorted = validParameter.validFile(parameters, "sorted", false); if (sorted == "not found") { sorted = "T"; }
+ string temp = validParameter.validFile(parameters, "numotu", false); if (temp == "not found") { temp = "0"; }
+ convert(temp, numOTU);
+
+ temp = validParameter.validFile(parameters, "fontsize", false); if (temp == "not found") { temp = "24"; }
+ convert(temp, fontSize);
+
+ sorted = validParameter.validFile(parameters, "sorted", false);
+ if (sorted == "not found") {
+ //if numOTU is used change default
+ if (numOTU != 0) { sorted = "topotu"; }
+ else { sorted = "shared"; }
+ }
scale = validParameter.validFile(parameters, "scale", false); if (scale == "not found") { scale = "log10"; }
+ if ((sorted != "none") && (sorted != "shared") && (sorted != "topotu") && (sorted != "topgroup")) { m->mothurOut(sorted + " is not a valid sorting option. Sorted options are: none, shared, topotu, topgroup"); m->mothurOutEndLine(); abort=true; }
}
}
void HeatMapCommand::help(){
try {
m->mothurOut("The heatmap.bin command can only be executed after a successful read.otu command.\n");
- m->mothurOut("The heatmap.bin command parameters are groups, sorted, scale label. No parameters are required.\n");
+ m->mothurOut("The heatmap.bin command parameters are groups, sorted, scale, numotu, fontsize and label. No parameters are required.\n");
m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included in your heatmap.\n");
- m->mothurOut("The sorted parameter allows you to choose to see the file with the shared otus at the top or the otus in the order they appear in your input file. \n");
+ m->mothurOut("The sorted parameter allows you to order the otus displayed, default=shared, meaning display the shared otus first. Other options for sorted are none, meaning the exact representation of your otus, \n");
+ m->mothurOut("topotu, meaning the otus with the greatest abundance when totaled across groups, topgroup, meaning the top otus for each group. \n");
m->mothurOut("The scale parameter allows you to choose the range of color your bin information will be displayed with.\n");
+ m->mothurOut("The numotu parameter allows you to display only the top N otus, by default all the otus are displayed. You could choose to look at the top 10, by setting numotu=10. The default for sorted is topotu when numotu is used.\n");
m->mothurOut("The group names are separated by dashes. The label parameter allows you to select what distance levels you would like a heatmap created for, and are also separated by dashes.\n");
+ m->mothurOut("The fontsize parameter allows you to adjust the font size of the picture created, default=24.\n");
m->mothurOut("The heatmap.bin command should be in the following format: heatmap.bin(groups=yourGroups, sorted=yourSorted, label=yourLabels).\n");
m->mothurOut("Example heatmap.bin(groups=A-B-C, sorted=F, scale=log10).\n");
m->mothurOut("The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n");
- m->mothurOut("The default value for sorted is T meaning you want the shared otus on top, you may change it to F meaning the exact representation of your input file.\n");
m->mothurOut("The default value for scale is log10; your other options are log2 and linear.\n");
m->mothurOut("The heatmap.bin command outputs a .svg file for each label you specify.\n");
m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
if (abort == true) { return 0; }
- heatmap = new HeatMap(sorted, scale, outputDir);
+ heatmap = new HeatMap(sorted, scale, numOTU, fontSize, outputDir);
format = globaldata->getFormat();
string lastLabel;
set<string> labels; //holds labels to be used
string format, groups, sorted, scale, label, outputDir;
vector<string> Groups, outputNames;
+ int numOTU, fontSize;
};
double jack, jacklci, jackhci;
- double maxRank = (double)rank->getMaxRank();
+ int maxRank = (double)rank->getMaxRank();
int S = rank->getNumBins();
double N[maxOrder+1];
double LogSD::logS(double x){
return -(1-x)*log(1-x)/x;
}
+/***********************************************************************/
EstOutput LogSD::getValues(SAbundVector* rank){
try {
SAbundVector *rank = &rankw;*/
data.resize(3,0);
- int numInd = rank->getNumSeqs();
- int numSpec = rank->getNumBins();
+ double numInd = rank->getNumSeqs();
+ double numSpec = rank->getNumBins();
double snRatio = (double)numSpec/numInd;
double x = .5;
double step = .4999999999;
}
double alpha = numInd*(1-x)/x;
- int oct = 1;
+ double oct = 1;
double octSumObs = 0;
double sumObs = 0;
double octSumExp = 0;
}
/***********************************************************************/
+inline void gobble(istringstream& f){
+
+ char d;
+ while(isspace(d=f.get())) {;}
+ f.putback(d);
+
+}
+
+/***********************************************************************/
+
inline string getline(istringstream& fileHandle) {
try {
vector<unsigned long int> positions;
ifstream inFASTA;
openInputFile(filename, inFASTA);
-
+
string input;
while(!inFASTA.eof()){
- input = getline(inFASTA); gobble(inFASTA);
+ input = getline(inFASTA);
if (input.length() != 0) {
if(input[0] == '>'){ unsigned long int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1); }
}
+ gobble(inFASTA); //has to be here since windows line endings are 2 characters and mess up the positions
}
inFASTA.close();
string input;
while(!in.eof()){
unsigned long int lastpos = in.tellg();
- input = getline(in); gobble(in);
+ input = getline(in);
if (input.length() != 0) {
unsigned long int pos = in.tellg();
if (pos != -1) { positions.push_back(pos - input.length() - 1); }
else { positions.push_back(lastpos); }
}
+ gobble(in); //has to be here since windows line endings are 2 characters and mess up the positions
}
in.close();
float npShannon = 0.0000;
double maxRank = (double)rank->getMaxRank();
- int sampled = rank->getNumSeqs();
+ double sampled = rank->getNumSeqs();
double Chat = 1.0000 - (double)rank->get(1)/(double)sampled;
int r3 = -1;
int r1Ind = 0;
int r3Ind = 0;
- int sumSpec = 0;
- int iqSum = 0;
+ double sumSpec = 0;
+ double iqSum = 0;
for(int i = 1; i < rank->size(); i++) {
if(r1 != -1 && r3 != -1)
i = rank->size();
string label;
label = shared[0]->getLabel();
- int fARare, fBRare, S12Rare, S12Abund, S12, f11, tempA, tempB, t10, t01, t11, t21, t12, t22, C12Numerator;
+ double fARare, fBRare, S12Rare, S12Abund, S12, f11, tempA, tempB, t10, t01, t11, t21, t12, t22, C12Numerator;
fARare = 0; fBRare = 0; S12Rare = 0; S12Abund = 0; S12 = 0; f11 = 0; t10 = 0; t01 = 0; t11= 0; t21= 0; t12= 0; t22= 0; C12Numerator = 0;
- float Sharedace, C12, part1, part2, part3, part4, part5, Gamma1, Gamma2, Gamma3;
+ double Sharedace, C12, part1, part2, part3, part4, part5, Gamma1, Gamma2, Gamma3;
/*fARare = number of OTUs with one individual found in A and less than or equal to 10 in B.
fBRare = number of OTUs with one individual found in B and less than or equal to 10 in A.
if (isnan(Gamma3) || isinf(Gamma3)) { Gamma3 = 0; }
if (isnan(part1) || isinf(part1)) { part1 = 0; }
if (isnan(part2) || isinf(part2)) { part2 = 0; }
-
+
part3 = fARare * Gamma1;
part4 = fBRare * Gamma2;
part5 = f11 * Gamma3;
EstOutput Anderberg::getValues(vector<SharedRAbundVector*> shared) {
try {
- int S1, S2, S12, tempA, tempB;
+ double S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
/*S1, S2 = number of OTUs observed or estimated in A and B
try {
data.resize(1,0);
- int sumSharedA, sumSharedB, sumSharedAB, tempA, tempB;
+ double sumSharedA, sumSharedB, sumSharedAB, tempA, tempB;
sumSharedA = 0; sumSharedB = 0; sumSharedAB = 0;
/*Xi, Yi = abundance of the ith shared OTU in A and B
EstOutput Jclass::getValues(vector<SharedRAbundVector*> shared) {
try {
- int S1, S2, S12, tempA, tempB;
+ double S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
/*S1, S2 = number of OTUs observed or estimated in A and B
* sharedkulczynski.cpp
* Mothur
*
- * Created by John Westcott on 3/24/09.
+ * Created by Sarah Westcott on 3/24/09.
* Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
*
*/
EstOutput Kulczynski::getValues(vector<SharedRAbundVector*> shared) {
try {
- int S1, S2, S12, tempA, tempB;
+ double S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
/*S1, S2 = number of OTUs observed or estimated in A and B
EstOutput KulczynskiCody::getValues(vector<SharedRAbundVector*> shared) {
try {
- int S1, S2, S12, tempA, tempB;
+ double S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
/*S1, S2 = number of OTUs observed or estimated in A and B
EstOutput Lennon::getValues(vector<SharedRAbundVector*> shared) {
try {
- int S1, S2, S12, tempA, tempB, min;
+ double S1, S2, S12, tempA, tempB, min;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0; min = 0;
/*S1, S2 = number of OTUs observed or estimated in A and B
try {
data.resize(1,0);
- float Atotal, Btotal, tempA, tempB;
+ double Atotal, Btotal, tempA, tempB;
Atotal = 0; Btotal = 0;
- float morhorn, sumSharedA, sumSharedB, a, b, d;
+ double morhorn, sumSharedA, sumSharedB, a, b, d;
morhorn = 0.0; sumSharedA = 0.0; sumSharedB = 0.0; a = 0.0; b = 0.0; d = 0.0;
//get the total values we need to calculate the theta denominator sums
EstOutput Ochiai::getValues(vector<SharedRAbundVector*> shared) {
try {
- int S1, S2, S12, tempA, tempB;
+ double S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
/*S1, S2 = number of OTUs observed or estimated in A and B
}
};
+struct individualFloat {
+ string group;
+ int bin;
+ float abundance;
+ bool operator()(const individual& i1, const individual& i2) {
+ return (i1.abundance > i2.abundance);
+ }
+};
+
+
#include "sabundvector.hpp"
#include "rabundvector.hpp"
#include "sharedrabundvector.h"
EstOutput SharedSobs::getValues(vector<SharedRAbundVector*> shared){
try {
data.resize(1,0);
- int observed = 0;
+ double observed = 0;
//loop through the species in each group
for (int k = 0; k < shared[0]->size(); k++) {
EstOutput SharedSobsCS::getValues(vector<SharedRAbundVector*> shared){
try {
data.resize(1,0);
- int observed = 0;
+ double observed = 0;
int numGroups = shared.size();
for (int i = 0; i < shared[0]->size(); i++) {
EstOutput SorClass::getValues(vector<SharedRAbundVector*> shared) {
try {
- int S1, S2, S12, tempA, tempB;
+ double S1, S2, S12, tempA, tempB;
S1 = 0; S2 = 0; S12 = 0; tempA = 0; tempB = 0;
/*S1, S2 = number of OTUs observed or estimated in A and B
try {
data.resize(1,0);
- int Atotal, Btotal, tempA, tempB;
+ double Atotal, Btotal, tempA, tempB;
Atotal = 0; Btotal = 0;
- float numerator, denominator, thetaN, sumSharedA, sumSharedB, a, b, d;
+ double numerator, denominator, thetaN, sumSharedA, sumSharedB, a, b, d;
numerator = 0.0; denominator = 0.0; thetaN = 0.0; sumSharedA = 0.0; sumSharedB = 0.0; a = 0.0; b = 0.0; d = 0.0;
//get the total values we need to calculate the theta denominator sums
try {
data.resize(3,0.0000);
- float Atotal = 0;
- float Btotal = 0;
- float thetaYC = 0;
- float pi = 0;
- float qi = 0;
- float a = 0;
- float b = 0;
- float d = 0;
+ double Atotal = 0;
+ double Btotal = 0;
+ double thetaYC = 0;
+ double pi = 0;
+ double qi = 0;
+ double a = 0;
+ double b = 0;
+ double d = 0;
- float sumPcubed = 0;
- float sumQcubed = 0;
- float sumPQsq = 0;
- float sumPsqQ = 0;
+ double sumPcubed = 0;
+ double sumQcubed = 0;
+ double sumPQsq = 0;
+ double sumPsqQ = 0;
//get the total values we need to calculate the theta denominator sums
for (int i = 0; i < shared[0]->size(); i++) {
if (isnan(thetaYC) || isinf(thetaYC)) { thetaYC = 0; }
- float varA = 4 / Atotal * (sumPcubed - a * a);
- float varB = 4 / Btotal * (sumQcubed - b * b);
- float varD = sumPQsq / Atotal + sumPsqQ / Btotal - d * d * (1/Atotal + 1/Btotal);
- float covAD = 2 / Atotal * (sumPsqQ - a * d);
- float covBD = 2 / Btotal * (sumPQsq - b* d);
+ double varA = 4 / Atotal * (sumPcubed - a * a);
+ double varB = 4 / Btotal * (sumQcubed - b * b);
+ double varD = sumPQsq / Atotal + sumPsqQ / Btotal - d * d * (1/Atotal + 1/Btotal);
+ double covAD = 2 / Atotal * (sumPsqQ - a * d);
+ double covBD = 2 / Btotal * (sumPQsq - b* d);
- float varT = d * d * (varA + varB) / pow(a + b - d, (float)4.0) + pow(a+b, (float)2.0) * varD / pow(a+b-d, (float)4.0)
- - 2.0 * (a + b) * d / pow(a + b - d, (float)4.0) * (covAD + covBD);
+ double varT = d * d * (varA + varB) / pow(a + b - d, (double)4.0) + pow(a+b, (double)2.0) * varD / pow(a+b-d, (double)4.0)
+ - 2.0 * (a + b) * d / pow(a + b - d, (double)4.0) * (covAD + covBD);
- float ci = 1.95 * sqrt(varT);
+ double ci = 1.95 * sqrt(varT);
data[0] = thetaYC;
data[1] = thetaYC - ci;
remove((fileroot + "rare.list").c_str());
remove((fileroot + "abund.list").c_str());
- wroteListFile["rare"] = false;
- wroteListFile["abund"] = false;
+ outputNames.push_back((fileroot + "rare.list"));
+ outputNames.push_back((fileroot + "abund.list"));
}else{
for (int i=0; i<Groups.size(); i++) {
remove((fileroot + Groups[i] + ".rare.list").c_str());
remove((fileroot + Groups[i] + ".abund.list").c_str());
-
- wroteListFile[(Groups[i] + ".rare")] = false;
- wroteListFile[(Groups[i] + ".abund")] = false;
+
+ outputNames.push_back((fileroot + Groups[i] + ".rare.list"));
+ outputNames.push_back((fileroot + Groups[i] + ".abund.list"));
}
}
delete input;
- for (map<string, bool>::iterator itBool = wroteListFile.begin(); itBool != wroteListFile.end(); itBool++) {
- string filename = fileroot + itBool->first;
- if ((itBool->first == "rare") || (itBool->first == "abund")) {
- filename = fileroot + itBool->first + ".list";
- }
- if (itBool->second) { //we wrote to this file
- outputNames.push_back(filename);
- }else{
- remove(filename.c_str());
- }
- }
-
if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
-
}else { //you are using the namefile to determine abundance
if (outputDir == "") { outputDir = hasPath(namefile); }
ofstream aout;
ofstream rout;
- if (rareNames.size() != 0) {
- string rare = outputDir + getRootName(getSimpleName(listfile)) + "rare.list";
- wroteListFile["rare"] = true;
- openOutputFileAppend(rare, rout);
- rout << thisList->getLabel() << '\t' << numRareBins << '\t';
- }
+ string rare = outputDir + getRootName(getSimpleName(listfile)) + "rare.list";
+ openOutputFileAppend(rare, rout);
+ outputNames.push_back(rare);
- if (abundNames.size() != 0) {
- string abund = outputDir + getRootName(getSimpleName(listfile)) + "abund.list";
- wroteListFile["abund"] = true;
- openOutputFileAppend(abund, aout);
- aout << thisList->getLabel() << '\t' << numAbundBins << '\t';
- }
+ string abund = outputDir + getRootName(getSimpleName(listfile)) + "abund.list";
+ openOutputFileAppend(abund, aout);
+ outputNames.push_back(abund);
+
+ if (rareNames.size() != 0) { rout << thisList->getLabel() << '\t' << numRareBins << '\t'; }
+ if (abundNames.size() != 0) { aout << thisList->getLabel() << '\t' << numAbundBins << '\t'; }
for (int i = 0; i < thisList->getNumBins(); i++) {
if (m->control_pressed) { break; }
else { aout << bin << '\t'; }
}
- if (rareNames.size() != 0) { rout << endl; rout.close(); }
- if (abundNames.size() != 0) { aout << endl; aout.close(); }
-
+ if (rareNames.size() != 0) { rout << endl; }
+ if (abundNames.size() != 0) { aout << endl; }
+
+ rout.close();
+ aout.close();
+
}else{ //parse names by abundance and group
string fileroot = outputDir + getRootName(getSimpleName(listfile));
ofstream* temp;
ofstream* temp2;
- map<string, bool> wroteFile;
+ //map<string, bool> wroteFile;
map<string, ofstream*> filehandles;
map<string, ofstream*>::iterator it3;
//end list vector
for (it3 = filehandles.begin(); it3 != filehandles.end(); it3++) {
(*(filehandles[it3->first])) << thisList->getLabel() << '\t' << groupNumBins[it3->first] << '\t' << groupVector[it3->first] << endl; // label numBins listvector for that group
- wroteListFile[it3->first] = true;
(*(filehandles[it3->first])).close();
delete it3->second;
}
ofstream aout;
ofstream rout;
+ string rare = outputDir + getRootName(getSimpleName(namefile)) + "rare.names";
+ openOutputFile(rare, rout);
+ outputNames.push_back(rare);
+
+ string abund = outputDir + getRootName(getSimpleName(namefile)) + "abund.names";
+ openOutputFile(abund, aout);
+ outputNames.push_back(abund);
+
if (rareNames.size() != 0) {
- string rare = outputDir + getRootName(getSimpleName(namefile)) + "rare.names";
- openOutputFile(rare, rout);
- outputNames.push_back(rare);
-
for (set<string>::iterator itRare = rareNames.begin(); itRare != rareNames.end(); itRare++) {
rout << (*itRare) << '\t' << nameMap[(*itRare)] << endl;
}
- rout.close();
}
+ rout.close();
if (abundNames.size() != 0) {
- string abund = outputDir + getRootName(getSimpleName(namefile)) + "abund.names";
- openOutputFile(abund, aout);
- outputNames.push_back(abund);
-
for (set<string>::iterator itAbund = abundNames.begin(); itAbund != abundNames.end(); itAbund++) {
aout << (*itAbund) << '\t' << nameMap[(*itAbund)] << endl;
}
- aout.close();
}
-
+ aout.close();
+
}else{ //parse names by abundance and group
string fileroot = outputDir + getRootName(getSimpleName(namefile));
ofstream* temp;
ofstream* temp2;
- map<string, bool> wroteFile;
map<string, ofstream*> filehandles;
map<string, ofstream*>::iterator it3;
openOutputFile(fileroot + Groups[i] + ".rare.names", *(filehandles[Groups[i]+".rare"]));
openOutputFile(fileroot + Groups[i] + ".abund.names", *(filehandles[Groups[i]+".abund"]));
-
- wroteFile[Groups[i] + ".rare"] = false;
- wroteFile[Groups[i] + ".abund"] = false;
}
for (map<string, string>::iterator itName = nameMap.begin(); itName != nameMap.end(); itName++) {
}
}
- for (itout = outputStrings.begin(); itout != outputStrings.end(); itout++) {
- *(filehandles[itout->first]) << itout->second << endl;
- wroteFile[itout->first] = true;
- }
+ for (itout = outputStrings.begin(); itout != outputStrings.end(); itout++) { *(filehandles[itout->first]) << itout->second << endl; }
}
for (it3 = filehandles.begin(); it3 != filehandles.end(); it3++) {
(*(filehandles[it3->first])).close();
- if (wroteFile[it3->first] == true) { outputNames.push_back(fileroot + it3->first + ".names"); }
- else { remove((it3->first).c_str()); }
+ outputNames.push_back(fileroot + it3->first + ".names");
delete it3->second;
}
}
ofstream aout;
ofstream rout;
- if (rareNames.size() != 0) {
- string rare = outputDir + getRootName(getSimpleName(inputFile)) + tag + "rare.accnos";
- openOutputFile(rare, rout);
- outputNames.push_back(rare);
-
- for (set<string>::iterator itRare = rareNames.begin(); itRare != rareNames.end(); itRare++) {
- rout << (*itRare) << endl;
- }
- rout.close();
+
+ string rare = outputDir + getRootName(getSimpleName(inputFile)) + tag + "rare.accnos";
+ openOutputFile(rare, rout);
+ outputNames.push_back(rare);
+
+ for (set<string>::iterator itRare = rareNames.begin(); itRare != rareNames.end(); itRare++) {
+ rout << (*itRare) << endl;
}
+ rout.close();
+
+ string abund = outputDir + getRootName(getSimpleName(inputFile)) + tag + "abund.accnos";
+ openOutputFile(abund, aout);
+ outputNames.push_back(abund);
- if (abundNames.size() != 0) {
- string abund = outputDir + getRootName(getSimpleName(inputFile)) + tag + "abund.accnos";
- openOutputFile(abund, aout);
- outputNames.push_back(abund);
-
- for (set<string>::iterator itAbund = abundNames.begin(); itAbund != abundNames.end(); itAbund++) {
- aout << (*itAbund) << endl;
- }
- aout.close();
+ for (set<string>::iterator itAbund = abundNames.begin(); itAbund != abundNames.end(); itAbund++) {
+ aout << (*itAbund) << endl;
}
+ aout.close();
+
}else{ //parse names by abundance and group
string fileroot = outputDir + getRootName(getSimpleName(inputFile));
ofstream* temp;
ofstream* temp2;
- map<string, bool> wroteFile;
map<string, ofstream*> filehandles;
map<string, ofstream*>::iterator it3;
openOutputFile(fileroot + tag + Groups[i] + ".rare.accnos", *(filehandles[Groups[i]+".rare"]));
openOutputFile(fileroot + tag + Groups[i] + ".abund.accnos", *(filehandles[Groups[i]+".abund"]));
-
- wroteFile[Groups[i] + ".rare"] = false;
- wroteFile[Groups[i] + ".abund"] = false;
}
//write rare
if (inUsersGroups(group, Groups)) { //only add if this is in a group we want
*(filehandles[group+".rare"]) << *itRare << endl;
- wroteFile[group+".rare"] = true;
}
}
if (inUsersGroups(group, Groups)) { //only add if this is in a group we want
*(filehandles[group+".abund"]) << *itAbund << endl;
- wroteFile[group+".abund"] = true;
}
}
//close files
for (it3 = filehandles.begin(); it3 != filehandles.end(); it3++) {
(*(filehandles[it3->first])).close();
- if (wroteFile[it3->first] == true) { outputNames.push_back(fileroot + tag + it3->first + ".accnos"); }
- else { remove((fileroot + tag + it3->first + ".accnos").c_str()); }
+ outputNames.push_back(fileroot + tag + it3->first + ".accnos");
delete it3->second;
}
}
ofstream aout;
ofstream rout;
- if (rareNames.size() != 0) {
- string rare = outputDir + getRootName(getSimpleName(groupfile)) + tag + "rare.groups";
- openOutputFile(rare, rout);
- outputNames.push_back(rare);
- }
+ string rare = outputDir + getRootName(getSimpleName(groupfile)) + tag + "rare.groups";
+ openOutputFile(rare, rout);
+ outputNames.push_back(rare);
+
+ string abund = outputDir + getRootName(getSimpleName(groupfile)) + tag + "abund.groups";
+ openOutputFile(abund, aout);
+ outputNames.push_back(abund);
- if (abundNames.size() != 0) {
- string abund = outputDir + getRootName(getSimpleName(groupfile)) + tag + "abund.groups";
- openOutputFile(abund, aout);
- outputNames.push_back(abund);
- }
-
-
for (map<string, string>::iterator itName = nameMap.begin(); itName != nameMap.end(); itName++) {
vector<string> names;
splitAtComma(itName->second, names); //parses bin into individual sequence names
}
}
- if (rareNames.size() != 0) { rout.close(); }
- if (abundNames.size() != 0) { aout.close(); }
+ rout.close();
+ aout.close();
}else{ //parse names by abundance and group
string fileroot = outputDir + getRootName(getSimpleName(groupfile));
ofstream* temp;
ofstream* temp2;
- map<string, bool> wroteFile;
map<string, ofstream*> filehandles;
map<string, ofstream*>::iterator it3;
openOutputFile(fileroot + tag + Groups[i] + ".rare.groups", *(filehandles[Groups[i]+".rare"]));
openOutputFile(fileroot + tag + Groups[i] + ".abund.groups", *(filehandles[Groups[i]+".abund"]));
-
- wroteFile[Groups[i] + ".rare"] = false;
- wroteFile[Groups[i] + ".abund"] = false;
}
for (map<string, string>::iterator itName = nameMap.begin(); itName != nameMap.end(); itName++) {
if (inUsersGroups(group, Groups)) { //only add if this is in a group we want
*(filehandles[group+rareAbund]) << names[i] << '\t' << group << endl;
- wroteFile[group+rareAbund] = true;
}
}
}
for (it3 = filehandles.begin(); it3 != filehandles.end(); it3++) {
(*(filehandles[it3->first])).close();
- if (wroteFile[it3->first] == true) { outputNames.push_back(fileroot + tag + it3->first + ".groups"); }
- else { remove((fileroot + tag + it3->first + ".groups").c_str()); }
+ outputNames.push_back(fileroot + tag + it3->first + ".groups");
delete it3->second;
}
}
ofstream aout;
ofstream rout;
- if (rareNames.size() != 0) {
- string rare = outputDir + getRootName(getSimpleName(fastafile)) + tag + "rare.fasta";
- openOutputFile(rare, rout);
- outputNames.push_back(rare);
- }
-
- if (abundNames.size() != 0) {
- string abund = outputDir + getRootName(getSimpleName(fastafile)) + tag + "abund.fasta";
- openOutputFile(abund, aout);
- outputNames.push_back(abund);
- }
-
-
+ string rare = outputDir + getRootName(getSimpleName(fastafile)) + tag + "rare.fasta";
+ openOutputFile(rare, rout);
+ outputNames.push_back(rare);
+
+ string abund = outputDir + getRootName(getSimpleName(fastafile)) + tag + "abund.fasta";
+ openOutputFile(abund, aout);
+ outputNames.push_back(abund);
+
//open input file
ifstream in;
openInputFile(fastafile, in);
}
}
in.close();
- if (rareNames.size() != 0) { rout.close(); }
- if (abundNames.size() != 0) { aout.close(); }
+ rout.close();
+ aout.close();
}else{ //parse names by abundance and group
string fileroot = outputDir + getRootName(getSimpleName(fastafile));
ofstream* temp;
ofstream* temp2;
- map<string, bool> wroteFile;
map<string, ofstream*> filehandles;
map<string, ofstream*>::iterator it3;
openOutputFile(fileroot + tag + Groups[i] + ".rare.fasta", *(filehandles[Groups[i]+".rare"]));
openOutputFile(fileroot + tag + Groups[i] + ".abund.fasta", *(filehandles[Groups[i]+".abund"]));
-
- wroteFile[Groups[i] + ".rare"] = false;
- wroteFile[Groups[i] + ".abund"] = false;
}
//open input file
if (inUsersGroups(group, Groups)) { //only add if this is in a group we want
seq.printSequence(*(filehandles[group+rareAbund]));
- wroteFile[group+rareAbund] = true;
}else if(group == "not found") {
m->mothurOut(seq.getName() + " is not in your groupfile. Ignoring."); m->mothurOutEndLine();
}
for (it3 = filehandles.begin(); it3 != filehandles.end(); it3++) {
(*(filehandles[it3->first])).close();
- if (wroteFile[it3->first] == true) { outputNames.push_back(fileroot + tag + it3->first + ".fasta"); }
- else { remove((fileroot + tag + it3->first + ".fasta").c_str()); }
+ outputNames.push_back(fileroot + tag + it3->first + ".fasta");
delete it3->second;
}
}
vector<string> Groups;
bool abort, allLines, accnos;
int cutoff;
- map<string, bool> wroteListFile;
+ //map<string, bool> wroteListFile;
map<string, string> nameMap;
projects / mothur.git / commitdiff