From cad05a21b084833b07808c1586e755be48fe7e1a Mon Sep 17 00:00:00 2001 From: westcott Date: Fri, 11 Jun 2010 15:37:05 +0000 Subject: [PATCH] fixed bug in cluster.split with classify method --- clustersplitcommand.cpp | 67 +++++++++++++++++------ clustersplitcommand.h | 2 +- decalc.cpp | 2 +- readtree.cpp | 85 ++++++++++++++++++++++-------- splitmatrix.cpp | 114 ++++++++++++++++++++++++++++++++++------ splitmatrix.h | 8 ++- 6 files changed, 219 insertions(+), 59 deletions(-) diff --git a/clustersplitcommand.cpp b/clustersplitcommand.cpp index c3faab7..363b695 100644 --- a/clustersplitcommand.cpp +++ b/clustersplitcommand.cpp @@ -21,13 +21,14 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { try{ globaldata = GlobalData::getInstance(); abort = false; + format = ""; //allow user to run help if(option == "help") { help(); abort = true; } else { //valid paramters for this command - string Array[] = {"phylip","column","name","cutoff","precision","method","splitmethod","taxonomy","taxlevel","large","showabund","timing","hard","processors","outputdir","inputdir"}; + string Array[] = {"fasta","phylip","column","name","cutoff","precision","method","splitmethod","taxonomy","taxlevel","large","showabund","timing","hard","processors","outputdir","inputdir"}; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); OptionParser parser(option); @@ -84,6 +85,14 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { //if the user has not given a path then, add inputdir. else leave path alone. if (path == "") { parameters["taxonomy"] = inputDir + it->second; } } + + it = parameters.find("fasta"); + //user has given a template file + if(it != parameters.end()){ + path = hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["fasta"] = inputDir + it->second; } + } } //check for required parameters @@ -101,16 +110,25 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { if (namefile == "not open") { abort = true; } else if (namefile == "not found") { namefile = ""; } + fastafile = validParameter.validFile(parameters, "fasta", true); + if (fastafile == "not open") { abort = true; } + else if (fastafile == "not found") { fastafile = ""; } + else { splitmethod = "fasta"; } + taxFile = validParameter.validFile(parameters, "taxonomy", true); if (taxFile == "not open") { abort = true; } else if (taxFile == "not found") { taxFile = ""; } - if ((phylipfile == "") && (columnfile == "")) { m->mothurOut("When executing a cluster.split command you must enter a phylip or a column."); m->mothurOutEndLine(); abort = true; } - else if ((phylipfile != "") && (columnfile != "")) { m->mothurOut("When executing a cluster.split command you must enter ONLY ONE of the following: phylip or column."); m->mothurOutEndLine(); abort = true; } + if ((phylipfile == "") && (columnfile == "") && (fastafile == "")) { m->mothurOut("When executing a cluster.split command you must enter a phylip or a column or fastafile."); m->mothurOutEndLine(); abort = true; } + else if ((phylipfile != "") && (columnfile != "") && (fastafile != "")) { m->mothurOut("When executing a cluster.split command you must enter ONLY ONE of the following: fasta, phylip or column."); m->mothurOutEndLine(); abort = true; } if (columnfile != "") { if (namefile == "") { m->mothurOut("You need to provide a namefile if you are going to use the column format."); m->mothurOutEndLine(); abort = true; } } + + if (fastafile != "") { + if (taxFile == "") { m->mothurOut("You need to provide a taxonomy file if you are using a fasta file to generate the split."); m->mothurOutEndLine(); abort = true; } + } //check for optional parameter and set defaults // ...at some point should added some additional type checking... @@ -131,7 +149,11 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; } convert(temp, processors); - splitmethod = validParameter.validFile(parameters, "splitmethod", false); if (splitmethod == "not found") { splitmethod = "distance"; } + temp = validParameter.validFile(parameters, "splitmethod", false); + if (splitmethod != "fasta") { + if (temp == "not found") { splitmethod = "distance"; } + else { splitmethod = temp; } + } temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found") { temp = "10"; } convert(temp, cutoff); @@ -145,8 +167,8 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { if ((method == "furthest") || (method == "nearest") || (method == "average")) { } else { m->mothurOut("Not a valid clustering method. Valid clustering algorithms are furthest, nearest or average."); m->mothurOutEndLine(); abort = true; } - if ((splitmethod == "distance") || (splitmethod == "classify")) { } - else { m->mothurOut("Not a valid splitting method. Valid splitting algorithms are distance or classify."); m->mothurOutEndLine(); abort = true; } + if ((splitmethod == "distance") || (splitmethod == "classify") || (splitmethod == "fasta")) { } + else { m->mothurOut(splitmethod + " is not a valid splitting method. Valid splitting algorithms are distance, classify or fasta."); m->mothurOutEndLine(); abort = true; } if ((splitmethod == "classify") && (taxFile == "")) { m->mothurOut("You need to provide a taxonomy file if you are going to use the classify splitmethod."); m->mothurOutEndLine(); abort = true; } @@ -168,13 +190,20 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { void ClusterSplitCommand::help(){ try { - m->mothurOut("The cluster.split command parameter options are phylip, column, name, cutoff, precision, method, splitmethod, taxonomy, taxlevel, showabund, timing, hard, large, processors. Phylip or column and name are required.\n"); + m->mothurOut("The cluster.split command parameter options are fasta, phylip, column, name, cutoff, precision, method, splitmethod, taxonomy, taxlevel, showabund, timing, hard, large, processors. Fasta or Phylip or column and name are required.\n"); + m->mothurOut("The cluster.split command can split your files in 3 ways. Splitting by distance file, by classification, or by classification also using a fasta file. \n"); + m->mothurOut("For the distance file method, you need only provide your distance file and mothur will split the file into distinct groups. \n"); + m->mothurOut("For the classification method, you need to provide your distance file and taxonomy file, and set the splitmethod to classify. \n"); + m->mothurOut("You will also need to set the taxlevel you want to split by. mothur will split the sequence into distinct taxonomy groups, and split the distance file based on those groups. \n"); + m->mothurOut("For the classification method using a fasta file, you need to provide your fasta file and taxonomy file. \n"); + m->mothurOut("You will also need to set the taxlevel you want to split by. mothur will split the sequence into distinct taxonomy groups, and create distance files for each grouping. \n"); m->mothurOut("The phylip and column parameter allow you to enter your distance file. \n"); + m->mothurOut("The fasta parameter allows you to enter your aligned fasta file. \n"); m->mothurOut("The name parameter allows you to enter your name file and is required if your distance file is in column format. \n"); m->mothurOut("The cutoff parameter allow you to set the distance you want to cluster to, default is 10.0. \n"); m->mothurOut("The precision parameter allows you specify the precision of the precision of the distances outputted, default=100, meaning 2 decimal places. \n"); m->mothurOut("The method allows you to specify what clustering algorythm you want to use, default=furthest, option furthest, nearest, or average. \n"); - m->mothurOut("The splitmethod parameter allows you to specify how you want to split your distance file before you cluster, default=distance, options distance or classify. \n"); + m->mothurOut("The splitmethod parameter allows you to specify how you want to split your distance file before you cluster, default=distance, options distance, classify or fasta. \n"); m->mothurOut("The taxonomy parameter allows you to enter the taxonomy file for your sequences, this is only valid if you are using splitmethod=classify. Be sure your taxonomy file does not include the probability scores. \n"); m->mothurOut("The taxlevel parameter allows you to specify the taxonomy level you want to use to split the distance file, default=1. \n"); m->mothurOut("The large parameter allows you to indicate that your distance matrix is too large to fit in RAM. The default value is false.\n"); @@ -199,12 +228,14 @@ int ClusterSplitCommand::execute(){ try { if (abort == true) { return 0; } - + time_t estart; //****************** file prep work ******************************// //if user gave a phylip file convert to column file if (format == "phylip") { - + estart = time(NULL); + m->mothurOut("Converting to column format..."); m->mothurOutEndLine(); + ReadCluster* convert = new ReadCluster(distfile, cutoff, outputDir, false); NameAssignment* nameMap = NULL; @@ -229,16 +260,20 @@ int ClusterSplitCommand::execute(){ } delete listToMakeNameFile; delete convert; + + m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to convert the distance file."); m->mothurOutEndLine(); } if (m->control_pressed) { return 0; } - time_t estart = time(NULL); - m->mothurOut("Splitting the distance file..."); m->mothurOutEndLine(); + estart = time(NULL); + m->mothurOut("Splitting the file..."); m->mothurOutEndLine(); //split matrix into non-overlapping groups SplitMatrix* split; - if (splitmethod == "distance") { split = new SplitMatrix(distfile, namefile, taxFile, cutoff, splitmethod, large); } - else { split = new SplitMatrix(distfile, namefile, taxFile, taxLevelCutoff, splitmethod, large); } + if (splitmethod == "distance") { split = new SplitMatrix(distfile, namefile, taxFile, cutoff, splitmethod, large); } + else if (splitmethod == "classify") { split = new SplitMatrix(distfile, namefile, taxFile, taxLevelCutoff, splitmethod, large); } + else if (splitmethod == "fasta") { split = new SplitMatrix(fastafile, taxFile, taxLevelCutoff, splitmethod); } + else { m->mothurOut("Not a valid splitting method. Valid splitting algorithms are distance, classify or fasta."); m->mothurOutEndLine(); return 0; } split->split(); @@ -740,8 +775,8 @@ vector ClusterSplitCommand::cluster(vector< map > distNa listFileNames.clear(); return listFileNames; } - remove(thisDistFile.c_str()); - remove(thisNamefile.c_str()); + //remove(thisDistFile.c_str()); + //remove(thisNamefile.c_str()); } diff --git a/clustersplitcommand.h b/clustersplitcommand.h index c696329..dc46f99 100644 --- a/clustersplitcommand.h +++ b/clustersplitcommand.h @@ -32,7 +32,7 @@ private: vector processIDS; //processid vector outputNames; - string method, fileroot, tag, outputDir, phylipfile, columnfile, namefile, distfile, format, showabund, timing, splitmethod, taxFile; + string method, fileroot, tag, outputDir, phylipfile, columnfile, namefile, distfile, format, showabund, timing, splitmethod, taxFile, fastafile; double cutoff, splitcutoff; int precision, length, processors, taxLevelCutoff; bool print_start, abort, hard, large; diff --git a/decalc.cpp b/decalc.cpp index f5aa850..74cd865 100644 --- a/decalc.cpp +++ b/decalc.cpp @@ -428,7 +428,7 @@ vector< vector > DeCalculator::getQuantiles(vector seqs, float de = calcDE(obsi, exp); float dist = calcDist(query, subject, front, back); - //o << i << '\t' << j << '\t' << dist << '\t' << de << endl; + //cout << i << '\t' << j << '\t' << dist << '\t' << de << endl; dist = ceil(dist); quanMember newScore(de, i, j); diff --git a/readtree.cpp b/readtree.cpp index 1df49b5..c805edc 100644 --- a/readtree.cpp +++ b/readtree.cpp @@ -287,15 +287,55 @@ int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) { if (m->control_pressed) { return -1; } int c = readNodeChar(f); - + string k; +k = c; + cout << "at beginning = " << k < childrenNodes; + while(f.peek() != ')'){ + int child = readNewickInt(f, n, T); + if (child == -1) { return -1; } //reports an error in reading + + childrenNodes.push_back(child); + + //after a child you either have , or ), check for both + if(f.peek()==')'){ break; } + else if (f.peek()==',') { readSpecialChar(f,',',"comma"); } + else { string k; + k = f.peek(); + cout << "in here k = " << k << '\t' << f.tellg() <mothurOut("Error in tree, please correct."); m->mothurOutEndLine(); return -1; } + + //then force into 2 node structure + for (int i = 1; i < childrenNodes.size(); i++) { + + int lc, rc; + if (i == 1) { lc = childrenNodes[i-1]; rc = childrenNodes[i]; } + else { lc = n; rc = childrenNodes[i]; } + cout << i << '\t' << lc << '\t' << rc << endl; + T->tree[n].setChildren(lc,rc); + T->tree[lc].setParent(n); + T->tree[rc].setParent(n); + + T->printTree(); cout << endl; + n++; + } - readSpecialChar(f,',',"comma"); + //to account for extra ++ in looping + n--; + //int lc = readNewickInt(f, n, T); + //if (lc == -1) { return -1; } //reports an error in reading + + //readSpecialChar(f,',',"comma"); - int rc = readNewickInt(f, n, T); - if (rc == -1) { return -1; } //reports an error in reading + //int rc = readNewickInt(f, n, T); + //if (rc == -1) { return -1; } //reports an error in reading if(f.peek()==')'){ readSpecialChar(f,')',"right parenthesis"); @@ -303,7 +343,6 @@ int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) { c=filehandle.get(); while((c!=',') && (c != -1) && (c!= ':') && (c!=';')){ c=filehandle.get(); } filehandle.putback(c); - } if(f.peek() == ':'){ @@ -317,8 +356,10 @@ int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) { } //to account for multifurcating trees generated by fasttree, we are forcing them to be bifurcating - while(f.peek() == ','){ - + /*while(f.peek() == ','){ + string k; + k = f.peek(); + cout << "in here k = " << k << '\t' << f.tellg() <tree[n].setChildren(lc,rc); T->tree[lc].setParent(n); @@ -330,14 +371,11 @@ int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) { readSpecialChar(f,',',"comma"); - int rc = readNewickInt(f, n, T); - cout << n << '\t' << lc << '\t' << rc << endl; + rc = readNewickInt(f, n, T); + if (rc == -1) { return -1; } //reports an error in reading if(f.peek()==')'){ - string k; - k = f.peek(); - cout << "in here " << k << endl; readSpecialChar(f,')',"right parenthesis"); //to pass over labels in trees c=filehandle.get(); @@ -345,25 +383,26 @@ int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) { filehandle.putback(c); if(f.peek() == ':'){ - readSpecialChar(f,':',"colon"); + readSpecialChar(f,':',"colon"); - if(n >= numNodes){ m->mothurOut("Error: Too many nodes in input tree\n"); readOk = -1; return -1; } + if(n >= numNodes){ m->mothurOut("Error: Too many nodes in input tree\n"); readOk = -1; return -1; } - T->tree[n].setBranchLength(readBranchLength(f)); + T->tree[n].setBranchLength(readBranchLength(f)); }else{ T->tree[n].setBranchLength(0.0); } break; } - } - cout << n << '\t' << lc << '\t' << rc << endl; - T->tree[n].setChildren(lc,rc); - T->tree[lc].setParent(n); - T->tree[rc].setParent(n); - T->printTree(); cout << endl; + }*/ + + //T->tree[n].setChildren(lc,rc); + //T->tree[lc].setParent(n); + //T->tree[rc].setParent(n); + //T->printTree(); cout << endl; return n++; + }else{ f.putback(c); string name = ""; diff --git a/splitmatrix.cpp b/splitmatrix.cpp index 0929cb0..cd05cc4 100644 --- a/splitmatrix.cpp +++ b/splitmatrix.cpp @@ -9,6 +9,9 @@ #include "splitmatrix.h" #include "phylotree.h" +#include "sequencedb.h" +#include "onegapdist.h" +#include "dist.h" /***********************************************************************/ @@ -21,6 +24,15 @@ SplitMatrix::SplitMatrix(string distfile, string name, string tax, float c, stri taxFile = tax; large = l; } +/***********************************************************************/ + +SplitMatrix::SplitMatrix(string ffile, string tax, float c, string t){ + m = MothurOut::getInstance(); + fastafile = ffile; + taxFile = tax; + cutoff = c; + method = t; +} /***********************************************************************/ @@ -29,7 +41,7 @@ int SplitMatrix::split(){ if (method == "distance") { splitDistance(); - }else if (method == "classify") { + }else if ((method == "classify") || (method == "fasta")) { splitClassify(); }else { m->mothurOut("Unknown splitting method, aborting split."); m->mothurOutEndLine(); @@ -80,7 +92,6 @@ int SplitMatrix::splitClassify(){ string seqname, tax; while(!in.eof()){ in >> seqname >> tax; gobble(in); - phylo->addSeqToTree(seqname, tax); } in.close(); @@ -89,13 +100,13 @@ int SplitMatrix::splitClassify(){ //make sure the cutoff is not greater than maxlevel if (cutoff > phylo->getMaxLevel()) { m->mothurOut("splitcutoff is greater than the longest taxonomy, using " + toString(phylo->getMaxLevel())); m->mothurOutEndLine(); cutoff = phylo->getMaxLevel(); } - + //for each node in tree for (int i = 0; i < phylo->getNumNodes(); i++) { //is this node within the cutoff TaxNode taxon = phylo->get(i); - + if (taxon.level == cutoff) {//if yes, then create group containing this nodes sequences if (taxon.accessions.size() > 1) { //if this taxon just has one seq its a singleton for (int j = 0; j < taxon.accessions.size(); j++) { @@ -105,7 +116,71 @@ int SplitMatrix::splitClassify(){ } } } - + + delete phylo; + + if (method == "classify") { + splitDistanceFileByTax(seqGroup, numGroups); + }else { + createDistanceFilesFromTax(seqGroup, numGroups); + } + + + return 0; + + } + catch(exception& e) { + m->errorOut(e, "SplitMatrix", "splitClassify"); + exit(1); + } +} +/***********************************************************************/ +int SplitMatrix::createDistanceFilesFromTax(map& seqGroup, int numGroups){ + try { + map::iterator it; + map::iterator it2; + map seqIndexInFasta; + + //read fastafile + SequenceDB alignDB; + + ifstream filehandle; + openInputFile(fastafile, filehandle); + int numSeqs = 0; + while (!filehandle.eof()) { + //input sequence info into sequencedb + Sequence newSequence(filehandle); + + if (newSequence.getName() != "") { + alignDB.push_back(newSequence); + seqIndexInFasta[newSequence.getName()] = numSeqs; + numSeqs++; + } + + //takes care of white space + gobble(filehandle); + } + filehandle.close(); + + Dist* distCalculator = new oneGapDist(); + + +//still not done.... + + + return 0; + } + catch(exception& e) { + m->errorOut(e, "SplitMatrix", "createDistanceFilesFromTax"); + exit(1); + } +} +/***********************************************************************/ +int SplitMatrix::splitDistanceFileByTax(map& seqGroup, int numGroups){ + try { + map::iterator it; + map::iterator it2; + ifstream dFile; openInputFile(distFile, dFile); ofstream outFile; @@ -114,12 +189,15 @@ int SplitMatrix::splitClassify(){ remove((distFile + "." + toString(i) + ".temp").c_str()); } - //for buffering the io to improve speed //allow for 10 dists to be stored, then output. vector outputs; outputs.resize(numGroups, ""); vector numOutputs; numOutputs.resize(numGroups, 0); + //you can have a group made, but their may be no distances in the file for this group if the taxonomy file and distance file don't match + //this can occur if we have converted the phylip to column, since we reduce the size at that step by using the cutoff value + vector validDistances; validDistances.resize(numGroups, false); + //for each distance while(dFile){ string seqA, seqB; @@ -135,12 +213,13 @@ int SplitMatrix::splitClassify(){ if ((it != seqGroup.end()) && (it2 != seqGroup.end())) { //they are both not singletons if (it->second == it2->second) { //they are from the same group so add the distance - if (numOutputs[it->second] > 10) { + if (numOutputs[it->second] > 30) { openOutputFileAppend((distFile + "." + toString(it->second) + ".temp"), outFile); outFile << outputs[it->second] << seqA << '\t' << seqB << '\t' << dist << endl; outFile.close(); outputs[it->second] = ""; numOutputs[it->second] = 0; + validDistances[it->second] = true; }else{ outputs[it->second] += seqA + '\t' + seqB + '\t' + toString(dist) + '\n'; numOutputs[it->second]++; @@ -154,12 +233,13 @@ int SplitMatrix::splitClassify(){ remove((namefile + "." + toString(i) + ".temp").c_str()); //write out any remaining buffers - if (numOutputs[it->second] > 0) { + if (numOutputs[i] > 0) { openOutputFileAppend((distFile + "." + toString(i) + ".temp"), outFile); outFile << outputs[i]; outFile.close(); outputs[i] = ""; numOutputs[i] = 0; + validDistances[i] = true; } } @@ -195,14 +275,17 @@ int SplitMatrix::splitClassify(){ remove(singleton.c_str()); singleton = "none"; } - + for(int i=0;i temp; - temp[tempDistFile] = tempNameFile; - dists.push_back(temp); + map temp; + temp[tempDistFile] = tempNameFile; + dists.push_back(temp); + } } if (m->control_pressed) { @@ -214,10 +297,9 @@ int SplitMatrix::splitClassify(){ } return 0; - } catch(exception& e) { - m->errorOut(e, "SplitMatrix", "splitClassify"); + m->errorOut(e, "SplitMatrix", "splitDistanceFileByTax"); exit(1); } } diff --git a/splitmatrix.h b/splitmatrix.h index 95fe0dc..9467dc1 100644 --- a/splitmatrix.h +++ b/splitmatrix.h @@ -19,7 +19,9 @@ class SplitMatrix { public: - SplitMatrix(string, string, string, float, string, bool); //column formatted distance file, namesfile, cutoff, method + SplitMatrix(string, string, string, float, string, bool); //column formatted distance file, namesfile, cutoff, method, large + SplitMatrix(string, string, float, string); //fastafile, taxFile, cutoff, method + ~SplitMatrix(); int split(); vector< map > getDistanceFiles(); //returns map of distance files -> namefile sorted by distance file size @@ -28,7 +30,7 @@ class SplitMatrix { private: MothurOut* m; - string distFile, namefile, singleton, method, taxFile; + string distFile, namefile, singleton, method, taxFile, fastafile; vector< map< string, string> > dists; float cutoff; bool large; @@ -38,6 +40,8 @@ class SplitMatrix { int splitDistanceLarge(); int splitDistanceRAM(); int splitNames(vector >& groups); + int splitDistanceFileByTax(map&, int); + int createDistanceFilesFromTax(map&, int); }; /******************************************************/ -- 2.39.2