X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=splitmatrix.cpp;h=f6b5c4d81cf5446adc5e5253808380be7ebbf28d;hp=0929cb09167aebb2362bc6f29bf01c710d79ec8d;hb=d1c97b8c04bb75faca1e76ffad60b37a4d789d3d;hpb=2a7d1455e8cfe4f67a7173f3a7249762c5436217 diff --git a/splitmatrix.cpp b/splitmatrix.cpp index 0929cb0..f6b5c4d 100644 --- a/splitmatrix.cpp +++ b/splitmatrix.cpp @@ -9,18 +9,36 @@ #include "splitmatrix.h" #include "phylotree.h" +#include "distancecommand.h" +#include "seqsummarycommand.h" /***********************************************************************/ -SplitMatrix::SplitMatrix(string distfile, string name, string tax, float c, string t, bool l){ +SplitMatrix::SplitMatrix(string distfile, string name, string count, string tax, float c, string t, bool l){ m = MothurOut::getInstance(); distFile = distfile; cutoff = c; namefile = name; method = t; taxFile = tax; + countfile = count; large = l; } +/***********************************************************************/ + +SplitMatrix::SplitMatrix(string ffile, string name, string count, string tax, float c, float cu, string t, int p, bool cl, string output){ + m = MothurOut::getInstance(); + fastafile = ffile; + namefile = name; + countfile = count; + taxFile = tax; + cutoff = c; //tax level cutoff + distCutoff = cu; //for fasta method if you are creating distance matrix you need a cutoff for that + method = t; + processors = p; + classic = cl; + outputDir = output; +} /***********************************************************************/ @@ -29,12 +47,13 @@ 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(); map temp; - temp[distFile] = namefile; + if (namefile != "") { temp[distFile] = namefile; } + else { temp[distFile] = countfile; } dists.push_back(temp); } @@ -51,6 +70,8 @@ int SplitMatrix::splitDistance(){ if (large) { splitDistanceLarge(); } else { splitDistanceRAM(); } + + return 0; } catch(exception& e) { @@ -63,7 +84,7 @@ int SplitMatrix::splitDistance(){ int SplitMatrix::splitClassify(){ try { cutoff = int(cutoff); - + map seqGroup; map::iterator it; map::iterator it2; @@ -73,29 +94,25 @@ int SplitMatrix::splitClassify(){ //build tree from users taxonomy file PhyloTree* phylo = new PhyloTree(); - ifstream in; - openInputFile(taxFile, in); - - //read in users taxonomy file and add sequences to tree - string seqname, tax; - while(!in.eof()){ - in >> seqname >> tax; gobble(in); - - phylo->addSeqToTree(seqname, tax); - } - in.close(); + map temp; + m->readTax(taxFile, temp); + + for (map::iterator itTemp = temp.begin(); itTemp != temp.end();) { + phylo->addSeqToTree(itTemp->first, itTemp->second); + temp.erase(itTemp++); + } phylo->assignHeirarchyIDs(0); //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,29 +122,146 @@ int SplitMatrix::splitClassify(){ } } } - - ifstream dFile; - openInputFile(distFile, dFile); - ofstream outFile; + + 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 copyGroups = seqGroup; + map::iterator it; + set names; + for (int i = 0; i < numGroups; i++) { //remove old temp files, just in case - remove((distFile + "." + toString(i) + ".temp").c_str()); + m->mothurRemove((fastafile + "." + toString(i) + ".temp")); } + + ifstream in; + m->openInputFile(fastafile, in); + + //parse fastafile + ofstream outFile; + while (!in.eof()) { + Sequence query(in); m->gobble(in); + if (query.getName() != "") { + + it = seqGroup.find(query.getName()); + + //save names in case no namefile is given + if ((namefile == "") && (countfile == "")) { names.insert(query.getName()); } + + if (it != seqGroup.end()) { //not singleton + m->openOutputFileAppend((fastafile + "." + toString(it->second) + ".temp"), outFile); + query.printSequence(outFile); + outFile.close(); + + copyGroups.erase(query.getName()); + } + } + } + in.close(); + + //warn about sequence in groups that are not in fasta file + for(it = copyGroups.begin(); it != copyGroups.end(); it++) { + m->mothurOut("ERROR: " + it->first + " is missing from your fastafile. This could happen if your taxonomy file is not unique and your fastafile is, or it could indicate and error."); m->mothurOutEndLine(); + exit(1); + } + + copyGroups.clear(); + + //process each distance file + for (int i = 0; i < numGroups; i++) { + + string options = ""; + if (classic) { options = "fasta=" + (fastafile + "." + toString(i) + ".temp") + ", processors=" + toString(processors) + ", output=lt"; } + else { options = "fasta=" + (fastafile + "." + toString(i) + ".temp") + ", processors=" + toString(processors) + ", cutoff=" + toString(distCutoff); } + if (outputDir != "") { options += ", outputdir=" + outputDir; } + + m->mothurOut("/******************************************/"); m->mothurOutEndLine(); + + Command* command = new DistanceCommand(options); + + m->mothurOut("/******************************************/"); m->mothurOutEndLine(); + + command->execute(); + delete command; + + m->mothurRemove((fastafile + "." + toString(i) + ".temp")); + + //remove old names files just in case + if (namefile != "") { m->mothurRemove((namefile + "." + toString(i) + ".temp")); } + else { m->mothurRemove((countfile + "." + toString(i) + ".temp")); } + } + + //restore old fasta file name since dist.seqs overwrites it with the temp files + m->setFastaFile(fastafile); + + vector tempDistFiles; + for(int i=0;ihasPath(fastafile); } + string tempDistFile = ""; + if (classic) { tempDistFile = outputDir + m->getRootName(m->getSimpleName((fastafile + "." + toString(i) + ".temp"))) + "phylip.dist";} + else { tempDistFile = outputDir + m->getRootName(m->getSimpleName((fastafile + "." + toString(i) + ".temp"))) + "dist"; } + tempDistFiles.push_back(tempDistFile); + } + + splitNames(seqGroup, numGroups, tempDistFiles); + + if (m->control_pressed) { for (int i = 0; i < dists.size(); i++) { m->mothurRemove((dists[i].begin()->first)); m->mothurRemove((dists[i].begin()->second)); } dists.clear(); } + + 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; + + ofstream outFile; + ifstream dFile; + m->openInputFile(distFile, dFile); + for (int i = 0; i < numGroups; i++) { //remove old temp files, just in case + m->mothurRemove((distFile + "." + toString(i) + ".temp")); + } + //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; float dist; - if (m->control_pressed) { dFile.close(); for (int i = 0; i < numGroups; i++) { remove((distFile + "." + toString(i) + ".temp").c_str()); } } + if (m->control_pressed) { dFile.close(); for (int i = 0; i < numGroups; i++) { m->mothurRemove((distFile + "." + toString(i) + ".temp")); } } - dFile >> seqA >> seqB >> dist; gobble(dFile); + dFile >> seqA >> seqB >> dist; m->gobble(dFile); //if both sequences are in the same group then they are within the cutoff it = seqGroup.find(seqA); @@ -135,12 +269,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) { - openOutputFileAppend((distFile + "." + toString(it->second) + ".temp"), outFile); + if (numOutputs[it->second] > 30) { + m->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]++; @@ -149,75 +284,41 @@ int SplitMatrix::splitClassify(){ } } dFile.close(); - + + string inputFile = namefile; + if (countfile != "") { inputFile = countfile; } + + vector tempDistFiles; for (int i = 0; i < numGroups; i++) { //remove old temp files, just in case - remove((namefile + "." + toString(i) + ".temp").c_str()); + string tempDistFile = distFile + "." + toString(i) + ".temp"; + tempDistFiles.push_back(tempDistFile); + m->mothurRemove((inputFile + "." + toString(i) + ".temp")); //write out any remaining buffers - if (numOutputs[it->second] > 0) { - openOutputFileAppend((distFile + "." + toString(i) + ".temp"), outFile); + if (numOutputs[i] > 0) { + m->openOutputFileAppend((distFile + "." + toString(i) + ".temp"), outFile); outFile << outputs[i]; outFile.close(); outputs[i] = ""; numOutputs[i] = 0; + validDistances[i] = true; } } - ifstream bigNameFile; - openInputFile(namefile, bigNameFile); - - singleton = namefile + ".extra.temp"; - ofstream remainingNames; - openOutputFile(singleton, remainingNames); - - bool wroteExtra = false; - - string name, nameList; - while(!bigNameFile.eof()){ - bigNameFile >> name >> nameList; gobble(bigNameFile); - - //did this sequence get assigned a group - it = seqGroup.find(name); - - if (it != seqGroup.end()) { - openOutputFileAppend((namefile + "." + toString(it->second) + ".temp"), outFile); - outFile << name << '\t' << nameList << endl; - outFile.close(); - }else{ - wroteExtra = true; - remainingNames << name << '\t' << nameList << endl; - } - } - bigNameFile.close(); - remainingNames.close(); - - if (!wroteExtra) { - remove(singleton.c_str()); - singleton = "none"; - } - - for(int i=0;i temp; - temp[tempDistFile] = tempNameFile; - dists.push_back(temp); - } - + splitNames(seqGroup, numGroups, tempDistFiles); + if (m->control_pressed) { for (int i = 0; i < dists.size(); i++) { - remove((dists[i].begin()->first).c_str()); - remove((dists[i].begin()->second).c_str()); + m->mothurRemove((dists[i].begin()->first)); + m->mothurRemove((dists[i].begin()->second)); } dists.clear(); } return 0; - } catch(exception& e) { - m->errorOut(e, "SplitMatrix", "splitClassify"); + m->errorOut(e, "SplitMatrix", "splitDistanceFileByTax"); exit(1); } } @@ -234,9 +335,9 @@ int SplitMatrix::splitDistanceLarge(){ int numGroups = 0; - ofstream outFile; + //ofstream outFile; ifstream dFile; - openInputFile(distFile, dFile); + m->openInputFile(distFile, dFile); while(dFile){ string seqA, seqB; @@ -244,7 +345,7 @@ int SplitMatrix::splitDistanceLarge(){ dFile >> seqA >> seqB >> dist; - if (m->control_pressed) { dFile.close(); for(int i=0;i 0){ remove((distFile + "." + toString(i) + ".temp").c_str()); } } return 0; } + if (m->control_pressed) { dFile.close(); for(int i=0;i 0){ m->mothurRemove((distFile + "." + toString(i) + ".temp")); } } return 0; } if(dist < cutoff){ //cout << "in cutoff: " << dist << endl; @@ -310,6 +411,7 @@ int SplitMatrix::splitDistanceLarge(){ //have we reached the max buffer size if (numOutputs[groupID] > 60) { //write out sequence + ofstream outFile; outFile.open(fileName.c_str(), ios::app); outFile << outputs[groupID] << seqA << '\t' << seqB << '\t' << dist << endl; outFile.close(); @@ -336,7 +438,7 @@ int SplitMatrix::splitDistanceLarge(){ //if groupB is written to file it is above buffer size so read and write to new merged file if (wroteOutPut[groupIDB]) { string fileName2 = distFile + "." + toString(groupIDB) + ".temp"; - ifstream fileB(fileName2.c_str(), ios::ate); + /*ifstream fileB(fileName2.c_str(), ios::ate); outFile.open(fileName.c_str(), ios::app); @@ -371,17 +473,22 @@ int SplitMatrix::splitDistanceLarge(){ outFile << temp.substr(0, lastRead); delete memblock; - fileB.close(); - remove(fileName2.c_str()); + fileB.close();*/ + m->appendFiles(fileName2, fileName); + m->mothurRemove(fileName2); + //write out the merged memory if (numOutputs[groupID] > 60) { - outFile << outputs[groupID]; + ofstream tempOut; + m->openOutputFile(fileName, tempOut); + tempOut << outputs[groupID]; outputs[groupID] = ""; numOutputs[groupID] = 0; + tempOut.close(); } - outFile.close(); + //outFile.close(); wroteOutPut[groupID] = true; wroteOutPut[groupIDB] = false; @@ -396,7 +503,7 @@ int SplitMatrix::splitDistanceLarge(){ if (wroteOutPut[groupIDA]) { string fileName2 = distFile + "." + toString(groupIDA) + ".temp"; - ifstream fileB(fileName2.c_str(), ios::ate); + /*ifstream fileB(fileName2.c_str(), ios::ate); outFile.open(fileName.c_str(), ios::app); @@ -431,17 +538,21 @@ int SplitMatrix::splitDistanceLarge(){ delete memblock; - fileB.close(); - remove(fileName2.c_str()); + fileB.close();*/ + m->appendFiles(fileName2, fileName); + m->mothurRemove(fileName2); //write out the merged memory if (numOutputs[groupID] > 60) { - outFile << outputs[groupID]; + ofstream tempOut; + m->openOutputFile(fileName, tempOut); + tempOut << outputs[groupID]; outputs[groupID] = ""; numOutputs[groupID] = 0; + tempOut.close(); } - outFile.close(); + //outFile.close(); wroteOutPut[groupID] = true; wroteOutPut[groupIDA] = false; @@ -450,20 +561,33 @@ int SplitMatrix::splitDistanceLarge(){ } } } - gobble(dFile); + m->gobble(dFile); } dFile.close(); - + + vector tempDistFiles; for (int i = 0; i < numGroups; i++) { + string fileName = distFile + "." + toString(i) + ".temp"; + tempDistFiles.push_back(fileName); + //remove old names files just in case + if (numOutputs[i] > 0) { - string fileName = distFile + "." + toString(i) + ".temp"; + ofstream outFile; outFile.open(fileName.c_str(), ios::app); outFile << outputs[i]; outFile.close(); } } - - splitNames(groups); + + map seqGroup; + for (int i = 0; i < groups.size(); i++) { + for (set::iterator itNames = groups[i].begin(); itNames != groups[i].end();) { + seqGroup[*itNames] = i; + groups[i].erase(itNames++); + } + } + + splitNames(seqGroup, numGroups, tempDistFiles); return 0; } @@ -473,73 +597,104 @@ int SplitMatrix::splitDistanceLarge(){ } } //******************************************************************************************************************** -int SplitMatrix::splitNames(vector >& groups){ +int SplitMatrix::splitNames(map& seqGroup, int numGroups, vector& tempDistFiles){ try { - int numGroups = groups.size(); - - ifstream bigNameFile(namefile.c_str()); - if(!bigNameFile){ - cerr << "Error: We can't open the name file\n"; - exit(1); - } - - map nameMap; - string name, nameList; - while(bigNameFile){ - bigNameFile >> name >> nameList; - nameMap[name] = nameList; - gobble(bigNameFile); - } - bigNameFile.close(); - - for(int i=0;i 0){ - string fileName = namefile + "." + toString(i) + ".temp"; - ofstream smallNameFile(fileName.c_str(), ios::ate); - - for(set::iterator gIt=groups[i].begin();gIt!=groups[i].end();gIt++){ - map::iterator nIt = nameMap.find(*gIt); - if (nIt != nameMap.end()) { - smallNameFile << nIt->first << '\t' << nIt->second << endl; - nameMap.erase(nIt); - }else{ - m->mothurOut((*gIt) + " is in your distance file and not in your namefile. Please correct."); m->mothurOutEndLine(); exit(1); - } - } - smallNameFile.close(); - } - } - - //names of singletons - if (nameMap.size() != 0) { - singleton = namefile + ".extra.temp"; - ofstream remainingNames(singleton.c_str(), ios::ate); - for(map::iterator nIt=nameMap.begin();nIt!=nameMap.end();nIt++){ - remainingNames << nIt->first << '\t' << nIt->second << endl; - } - remainingNames.close(); - }else { singleton = "none"; } - + ofstream outFile; + map::iterator it; + + string inputFile = namefile; + if (countfile != "") { inputFile = countfile; } + + for(int i=0;imothurRemove((inputFile + "." + toString(i) + ".temp")); } + + singleton = inputFile + ".extra.temp"; + ofstream remainingNames; + m->openOutputFile(singleton, remainingNames); + + bool wroteExtra = false; + + ifstream bigNameFile; + m->openInputFile(inputFile, bigNameFile); + + //grab header line + string headers = ""; + if (countfile != "") { headers = m->getline(bigNameFile); m->gobble(bigNameFile); } + + string name, nameList; + while(!bigNameFile.eof()){ + bigNameFile >> name >> nameList; + m->getline(bigNameFile); m->gobble(bigNameFile); //extra getline is for rest of countfile line if groups are given. + + //did this sequence get assigned a group + it = seqGroup.find(name); + + if (it != seqGroup.end()) { + m->openOutputFileAppend((inputFile + "." + toString(it->second) + ".temp"), outFile); + outFile << name << '\t' << nameList << endl; + outFile.close(); + }else{ + wroteExtra = true; + remainingNames << name << '\t' << nameList << endl; + } + } + bigNameFile.close(); + for(int i=0;i 0){ - string tempNameFile = namefile + "." + toString(i) + ".temp"; - string tempDistFile = distFile + "." + toString(i) + ".temp"; - + string tempNameFile = inputFile + "." + toString(i) + ".temp"; + string tempDistFile = tempDistFiles[i]; + + //if there are valid distances + ifstream fileHandle; + fileHandle.open(tempDistFile.c_str()); + if(fileHandle) { + m->gobble(fileHandle); + if (!fileHandle.eof()) { //check map temp; + if (countfile != "") { + //add header + ofstream out; + string newtempNameFile = tempNameFile + "2"; + m->openOutputFile(newtempNameFile, out); + out << "Representative_Sequence\ttotal" << endl; + out.close(); + m->appendFiles(tempNameFile, newtempNameFile); + m->mothurRemove(tempNameFile); + m->renameFile(newtempNameFile, tempNameFile); + } temp[tempDistFile] = tempNameFile; dists.push_back(temp); + }else{ + ifstream in; + m->openInputFile(tempNameFile, in); + + while(!in.eof()) { + in >> name >> nameList; m->gobble(in); + wroteExtra = true; + remainingNames << name << '\t' << nameList << endl; + } + in.close(); + m->mothurRemove(tempNameFile); } + } + fileHandle.close(); } - if (m->control_pressed) { - for (int i = 0; i < dists.size(); i++) { - remove((dists[i].begin()->first).c_str()); - remove((dists[i].begin()->second).c_str()); - } - dists.clear(); - } + remainingNames.close(); + + if (!wroteExtra) { + m->mothurRemove(singleton); + singleton = "none"; + }else if (countfile != "") { + //add header + ofstream out; + string newtempNameFile = singleton + "2"; + m->openOutputFile(newtempNameFile, out); + out << "Representative_Sequence\ttotal" << endl; + out.close(); + m->appendFiles(singleton, newtempNameFile); + m->mothurRemove(singleton); + m->renameFile(newtempNameFile, singleton); + } return 0; } @@ -557,7 +712,7 @@ int SplitMatrix::splitDistanceRAM(){ int numGroups = 0; ifstream dFile; - openInputFile(distFile, dFile); + m->openInputFile(distFile, dFile); while(dFile){ string seqA, seqB; @@ -565,7 +720,7 @@ int SplitMatrix::splitDistanceRAM(){ dFile >> seqA >> seqB >> dist; - if (m->control_pressed) { dFile.close(); for(int i=0;i 0){ remove((distFile + "." + toString(i) + ".temp").c_str()); } } return 0; } + if (m->control_pressed) { dFile.close(); for(int i=0;i 0){ m->mothurRemove((distFile + "." + toString(i) + ".temp")); } } return 0; } if(dist < cutoff){ //cout << "in cutoff: " << dist << endl; @@ -640,21 +795,31 @@ int SplitMatrix::splitDistanceRAM(){ } } } - gobble(dFile); + m->gobble(dFile); } dFile.close(); + vector tempDistFiles; for (int i = 0; i < numGroups; i++) { + string fileName = distFile + "." + toString(i) + ".temp"; + tempDistFiles.push_back(fileName); if (outputs[i] != "") { ofstream outFile; - string fileName = distFile + "." + toString(i) + ".temp"; outFile.open(fileName.c_str(), ios::ate); outFile << outputs[i]; outFile.close(); } } - - splitNames(groups); + + map seqGroup; + for (int i = 0; i < groups.size(); i++) { + for (set::iterator itNames = groups[i].begin(); itNames != groups[i].end();) { + seqGroup[*itNames] = i; + groups[i].erase(itNames++); + } + } + + splitNames(seqGroup, numGroups, tempDistFiles); return 0; }