#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;
+}
/***********************************************************************/
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<string, string> temp;
- temp[distFile] = namefile;
+ if (namefile != "") { temp[distFile] = namefile; }
+ else { temp[distFile] = countfile; }
dists.push_back(temp);
}
if (large) { splitDistanceLarge(); }
else { splitDistanceRAM(); }
+
+ return 0;
}
catch(exception& e) {
int SplitMatrix::splitClassify(){
try {
cutoff = int(cutoff);
-
+
map<string, int> seqGroup;
map<string, int>::iterator it;
map<string, int>::iterator it2;
//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<string, string> temp;
+ m->readTax(taxFile, temp);
+
+ for (map<string, string>::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++) {
}
}
}
-
- 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<string, int>& seqGroup, int numGroups){
+ try {
+ map<string, int> copyGroups = seqGroup;
+ map<string, int>::iterator it;
+ set<string> 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<string> tempDistFiles;
+ for(int i=0;i<numGroups;i++){
+ if (outputDir == "") { outputDir = m->hasPath(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<string, int>& seqGroup, int numGroups){
+ try {
+ map<string, int>::iterator it;
+ map<string, int>::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<string> outputs; outputs.resize(numGroups, "");
vector<int> 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<bool> 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);
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]++;
}
}
dFile.close();
-
+
+ string inputFile = namefile;
+ if (countfile != "") { inputFile = countfile; }
+
+ vector<string> 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<numGroups;i++){
- string tempNameFile = namefile + "." + toString(i) + ".temp";
- string tempDistFile = distFile + "." + toString(i) + ".temp";
-
- map<string, string> 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);
}
}
int numGroups = 0;
- ofstream outFile;
+ //ofstream outFile;
ifstream dFile;
- openInputFile(distFile, dFile);
+ m->openInputFile(distFile, dFile);
while(dFile){
string seqA, seqB;
dFile >> seqA >> seqB >> dist;
- if (m->control_pressed) { dFile.close(); for(int i=0;i<numGroups;i++){ if(groups[i].size() > 0){ remove((distFile + "." + toString(i) + ".temp").c_str()); } } return 0; }
+ if (m->control_pressed) { dFile.close(); for(int i=0;i<numGroups;i++){ if(groups[i].size() > 0){ m->mothurRemove((distFile + "." + toString(i) + ".temp")); } } return 0; }
if(dist < cutoff){
//cout << "in cutoff: " << dist << endl;
//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();
//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);
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;
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);
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;
}
}
}
- gobble(dFile);
+ m->gobble(dFile);
}
dFile.close();
-
+
+ vector<string> 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<string, int> seqGroup;
+ for (int i = 0; i < groups.size(); i++) {
+ for (set<string>::iterator itNames = groups[i].begin(); itNames != groups[i].end();) {
+ seqGroup[*itNames] = i;
+ groups[i].erase(itNames++);
+ }
+ }
+
+ splitNames(seqGroup, numGroups, tempDistFiles);
return 0;
}
}
}
//********************************************************************************************************************
-int SplitMatrix::splitNames(vector<set<string> >& groups){
+int SplitMatrix::splitNames(map<string, int>& seqGroup, int numGroups, vector<string>& 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<string, string> nameMap;
- string name, nameList;
- while(bigNameFile){
- bigNameFile >> name >> nameList;
- nameMap[name] = nameList;
- gobble(bigNameFile);
- }
- bigNameFile.close();
-
- for(int i=0;i<numGroups;i++){ //parse names file to match distance files
- int numSeqsInGroup = groups[i].size();
-
- if(numSeqsInGroup > 0){
- string fileName = namefile + "." + toString(i) + ".temp";
- ofstream smallNameFile(fileName.c_str(), ios::ate);
-
- for(set<string>::iterator gIt=groups[i].begin();gIt!=groups[i].end();gIt++){
- map<string,string>::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<string,string>::iterator nIt=nameMap.begin();nIt!=nameMap.end();nIt++){
- remainingNames << nIt->first << '\t' << nIt->second << endl;
- }
- remainingNames.close();
- }else { singleton = "none"; }
-
+ ofstream outFile;
+ map<string, int>::iterator it;
+
+ string inputFile = namefile;
+ if (countfile != "") { inputFile = countfile; }
+
+ for(int i=0;i<numGroups;i++){ m->mothurRemove((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<numGroups;i++){
- if(groups[i].size() > 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<string, string> 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;
}
int numGroups = 0;
ifstream dFile;
- openInputFile(distFile, dFile);
+ m->openInputFile(distFile, dFile);
while(dFile){
string seqA, seqB;
dFile >> seqA >> seqB >> dist;
- if (m->control_pressed) { dFile.close(); for(int i=0;i<numGroups;i++){ if(groups[i].size() > 0){ remove((distFile + "." + toString(i) + ".temp").c_str()); } } return 0; }
+ if (m->control_pressed) { dFile.close(); for(int i=0;i<numGroups;i++){ if(groups[i].size() > 0){ m->mothurRemove((distFile + "." + toString(i) + ".temp")); } } return 0; }
if(dist < cutoff){
//cout << "in cutoff: " << dist << endl;
}
}
}
- gobble(dFile);
+ m->gobble(dFile);
}
dFile.close();
+ vector<string> 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<string, int> seqGroup;
+ for (int i = 0; i < groups.size(); i++) {
+ for (set<string>::iterator itNames = groups[i].begin(); itNames != groups[i].end();) {
+ seqGroup[*itNames] = i;
+ groups[i].erase(itNames++);
+ }
+ }
+
+ splitNames(seqGroup, numGroups, tempDistFiles);
return 0;
}