From 01f6af90c907264686304a5c684c702e94ff40ae Mon Sep 17 00:00:00 2001 From: Sarah Westcott Date: Tue, 22 May 2012 16:12:50 -0400 Subject: [PATCH] added debug output classify.seqs. worked on make.contigs command. added large parameter to shhh.flows that divides the flow file. --- classify.cpp | 10 +- makecontigscommand.cpp | 39 +++- makefile | 4 +- phylotree.cpp | 4 +- sequence.cpp | 4 +- shhhercommand.cpp | 407 ++++++++++++++++++++++++----------------- shhhercommand.h | 7 +- 7 files changed, 298 insertions(+), 177 deletions(-) diff --git a/classify.cpp b/classify.cpp index 2d01183..7726b3e 100644 --- a/classify.cpp +++ b/classify.cpp @@ -249,7 +249,8 @@ int Classify::readTaxonomy(string file) { m->mothurOutEndLine(); m->mothurOut("Reading in the " + file + " taxonomy...\t"); cout.flush(); - + if (m->debug) { m->mothurOut("[DEBUG]: Taxonomies read in...\n"); } + #ifdef USE_MPI int pid, num, processors; vector positions; @@ -298,6 +299,7 @@ int Classify::readTaxonomy(string file) { istringstream iss (tempBuf,istringstream::in); iss >> name; m->gobble(iss); iss >> taxInfo; + if (m->debug) { m->mothurOut("[DEBUG]: name = " + name + " tax = " + taxInfo + "\n"); } taxonomy[name] = taxInfo; phyloTree->addSeqToTree(name, taxInfo); } @@ -312,7 +314,9 @@ int Classify::readTaxonomy(string file) { while (!inTax.eof()) { inTax >> name; m->gobble(inTax); inTax >> taxInfo; - + + if (m->debug) { m->mothurOut("[DEBUG]: name = '" + name + "' tax = '" + taxInfo + "'\n"); } + taxonomy[name] = taxInfo; phyloTree->addSeqToTree(name, taxInfo); @@ -321,6 +325,8 @@ int Classify::readTaxonomy(string file) { } inTax.close(); #endif + + phyloTree->assignHeirarchyIDs(0); diff --git a/makecontigscommand.cpp b/makecontigscommand.cpp index f672840..523edd5 100644 --- a/makecontigscommand.cpp +++ b/makecontigscommand.cpp @@ -392,10 +392,32 @@ int MakeContigsCommand::driver(vector files, string outputFasta, string int numMismatches = 0; string seq1 = fSeq.getAligned(); string seq2 = rSeq.getAligned(); - vector scores1 = fQual.getQualityScores(); vector scores2 = rQual.getQualityScores(); - for (int i = 0; i < length; i++) { + + // if (num < 5) { cout << fSeq.getStartPos() << '\t' << fSeq.getEndPos() << '\t' << rSeq.getStartPos() << '\t' << rSeq.getEndPos() << endl; } + int overlapStart = fSeq.getStartPos(); + int seq2Start = rSeq.getStartPos(); + //bigger of the 2 starting positions is the location of the overlapping start + if (overlapStart < seq2Start) { //seq2 starts later so take from 0 to seq2Start from seq1 + overlapStart = seq2Start; + for (int i = 0; i < overlapStart; i++) { + contig += seq1[i]; + contigScores.push_back(scores1[ABaseMap[i]]); + } + }else { //seq1 starts later so take from 0 to overlapStart from seq2 + for (int i = 0; i < overlapStart; i++) { + contig += seq2[i]; + contigScores.push_back(scores2[BBaseMap[i]]); + } + } + + int seq1End = fSeq.getEndPos(); + int seq2End = rSeq.getEndPos(); + int overlapEnd = seq1End; + if (seq2End < overlapEnd) { overlapEnd = seq2End; } //smallest end position is where overlapping ends + + for (int i = overlapStart; i < overlapEnd; i++) { if (seq1[i] == seq2[i]) { //match, add base and choose highest score contig += seq1[i]; contigScores.push_back(scores1[ABaseMap[i]]); @@ -423,6 +445,19 @@ int MakeContigsCommand::driver(vector files, string outputFasta, string } } + if (seq1End < seq2End) { //seq1 ends before seq2 so take from overlap to length from seq2 + for (int i = overlapEnd; i < length; i++) { + contig += seq2[i]; + contigScores.push_back(scores2[BBaseMap[i]]); + } + }else { //seq2 ends before seq1 so take from overlap to length from seq1 + for (int i = overlapEnd; i < length; i++) { + contig += seq1[i]; + contigScores.push_back(scores1[ABaseMap[i]]); + } + + } + //if (num < 5) { cout << overlapStart << '\t' << overlapEnd << endl << seq1 << endl << seq2 << endl<< contig << endl; } //output outFasta << ">" << fSeq.getName() << endl << contig << endl; outQual << ">" << fSeq.getName() << endl; diff --git a/makefile b/makefile index e06c334..6f8a423 100644 --- a/makefile +++ b/makefile @@ -15,8 +15,8 @@ USEREADLINE ?= yes CYGWIN_BUILD ?= no USECOMPRESSION ?= no MOTHUR_FILES="\"Enter_your_default_path_here\"" -RELEASE_DATE = "\"4/30/2012\"" -VERSION = "\"1.25.0\"" +RELEASE_DATE = "\"5/14/2012\"" +VERSION = "\"1.25.1\"" FORTAN_COMPILER = gfortran FORTRAN_FLAGS = diff --git a/phylotree.cpp b/phylotree.cpp index a9ef6cb..4ed3d8c 100644 --- a/phylotree.cpp +++ b/phylotree.cpp @@ -619,13 +619,15 @@ bool PhyloTree::ErrorCheck(vector templateFileNames){ map::iterator itFind; map taxonomyFileNames = name2Taxonomy; + if (m->debug) { m->mothurOut("[DEBUG]: in error check. Numseqs in template = " + toString(templateFileNames.size()) + ". Numseqs in taxonomy = " + toString(taxonomyFileNames.size()) + ".\n"); } + for (int i = 0; i < templateFileNames.size(); i++) { itFind = taxonomyFileNames.find(templateFileNames[i]); if (itFind != taxonomyFileNames.end()) { //found it so erase it taxonomyFileNames.erase(itFind); }else { - m->mothurOut(templateFileNames[i] + " is in your template file and is not in your taxonomy file. Please correct."); m->mothurOutEndLine(); + m->mothurOut("'" +templateFileNames[i] + "' is in your template file and is not in your taxonomy file. Please correct."); m->mothurOutEndLine(); okay = false; } diff --git a/sequence.cpp b/sequence.cpp index d877f4b..a1e787b 100644 --- a/sequence.cpp +++ b/sequence.cpp @@ -603,7 +603,7 @@ int Sequence::getLongHomoPolymer(){ int Sequence::getStartPos(){ if(startPos == -1){ for(int j = 0; j < alignmentLength; j++) { - if(aligned[j] != '.'){ + if((aligned[j] != '.')&&(aligned[j] != '-')){ startPos = j + 1; break; } @@ -668,7 +668,7 @@ int Sequence::filterFromPos(int end){ int Sequence::getEndPos(){ if(endPos == -1){ for(int j=alignmentLength-1;j>=0;j--){ - if(aligned[j] != '.'){ + if((aligned[j] != '.')&&(aligned[j] != '-')){ endPos = j + 1; break; } diff --git a/shhhercommand.cpp b/shhhercommand.cpp index 508e2a9..a09825e 100644 --- a/shhhercommand.cpp +++ b/shhhercommand.cpp @@ -18,6 +18,7 @@ vector ShhherCommand::setParameters(){ CommandParameter pcutoff("cutoff", "Number", "", "0.01", "", "", "",false,false); parameters.push_back(pcutoff); CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); CommandParameter pmaxiter("maxiter", "Number", "", "1000", "", "", "",false,false); parameters.push_back(pmaxiter); + CommandParameter plarge("large", "Number", "", "-1", "", "", "",false,false); parameters.push_back(plarge); CommandParameter psigma("sigma", "Number", "", "60", "", "", "",false,false); parameters.push_back(psigma); CommandParameter pmindelta("mindelta", "Number", "", "0.000001", "", "", "",false,false); parameters.push_back(pmindelta); CommandParameter porder("order", "String", "", "", "", "", "",false,false); parameters.push_back(porder); @@ -307,6 +308,17 @@ ShhherCommand::ShhherCommand(string option) { temp = validParameter.validFile(parameters, "maxiter", false); if (temp == "not found"){ temp = "1000"; } m->mothurConvert(temp, maxIters); + + temp = validParameter.validFile(parameters, "large", false); if (temp == "not found"){ temp = "0"; } + m->mothurConvert(temp, largeSize); + if (largeSize != 0) { large = true; } + else { large = false; } + if (largeSize < 0) { m->mothurOut("The value of the large cannot be negative.\n"); } + +#ifdef USE_MPI + if (large) { m->mothurOut("The large parameter is not available with the MPI-Enabled version.\n"); large=false; } +#endif + temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found") { temp = "60"; } m->mothurConvert(temp, sigma); @@ -2014,170 +2026,248 @@ int ShhherCommand::createProcesses(vector filenames){ } /**************************************************************************************************/ -int ShhherCommand::driver(vector filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){ +vector ShhherCommand::parseFlowFiles(string filename){ try { + vector files; + int count = 0; - for(int i=start;icontrol_pressed) { break; } - - string flowFileName = filenames[i]; - - m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n"); - m->mothurOut("Reading flowgrams...\n"); - - vector seqNameVector; - vector lengths; - vector flowDataIntI; - vector flowDataPrI; - map nameMap; - vector uniqueFlowgrams; - vector uniqueCount; - vector mapSeqToUnique; - vector mapUniqueToSeq; - vector uniqueLengths; - int numFlowCells; - - int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells); - - if (m->control_pressed) { break; } - - m->mothurOut("Identifying unique flowgrams...\n"); - int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI); - - if (m->control_pressed) { break; } - - m->mothurOut("Calculating distances between flowgrams...\n"); - string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist"; - unsigned long long begTime = time(NULL); - double begClock = clock(); + ifstream in; + m->openInputFile(filename, in); - flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI); - - m->mothurOutEndLine(); - m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n'); - + int thisNumFLows = 0; + in >> thisNumFLows; m->gobble(in); + + while (!in.eof()) { + if (m->control_pressed) { break; } - string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names"; - createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq); - - if (m->control_pressed) { break; } - - m->mothurOut("\nClustering flowgrams...\n"); - string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list"; - cluster(listFileName, distFileName, namesFileName); - - if (m->control_pressed) { break; } + ofstream out; + string outputFileName = filename + toString(count) + ".temp"; + m->openOutputFile(outputFileName, out); + out << thisNumFLows << endl; + files.push_back(outputFileName); - vector otuData; - vector cumNumSeqs; - vector nSeqsPerOTU; - vector > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices - vector > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI - vector seqNumber; //tMaster->anP: the sequence id number sorted by OTU - vector seqIndex; //tMaster->anI; the index that corresponds to seqNumber + for (int i = 0; i < largeSize; i++) { + if (in.eof()) { break; } + string line = m->getline(in); + out << line << endl; + } + out.close(); + count++; + } + in.close(); + + if (m->control_pressed) { for (int i = 0; i < files.size(); i++) { m->mothurRemove(files[i]); } files.clear(); } + + return files; + } + catch(exception& e) { + m->errorOut(e, "ShhherCommand", "parseFlowFiles"); + exit(1); + } +} +/**************************************************************************************************/ - - int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap); +int ShhherCommand::driver(vector filenames, string thisCompositeFASTAFileName, string thisCompositeNamesFileName, int start, int end){ + try { + + for(int i=start;icontrol_pressed) { break; } - m->mothurRemove(distFileName); - m->mothurRemove(namesFileName); - m->mothurRemove(listFileName); - - vector dist; //adDist - distance of sequences to centroids - vector change; //did the centroid sequence change? 0 = no; 1 = yes - vector centroids; //the representative flowgram for each cluster m - vector weight; - vector singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs) - vector nSeqsBreaks; - vector nOTUsBreaks; + vector theseFlowFileNames; theseFlowFileNames.push_back(filenames[i]); + if (large) { theseFlowFileNames = parseFlowFiles(filenames[i]); } - dist.assign(numSeqs * numOTUs, 0); - change.assign(numOTUs, 1); - centroids.assign(numOTUs, -1); - weight.assign(numOTUs, 0); - singleTau.assign(numSeqs, 1.0); - - nSeqsBreaks.assign(2, 0); - nOTUsBreaks.assign(2, 0); + if (m->control_pressed) { break; } - nSeqsBreaks[0] = 0; - nSeqsBreaks[1] = numSeqs; - nOTUsBreaks[1] = numOTUs; - - if (m->control_pressed) { break; } - - double maxDelta = 0; - int iter = 0; - - begClock = clock(); - begTime = time(NULL); + double begClock = clock(); + unsigned long long begTime; - m->mothurOut("\nDenoising flowgrams...\n"); - m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n"); - - while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){ - - if (m->control_pressed) { break; } - - double cycClock = clock(); - unsigned long long cycTime = time(NULL); - fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); - - if (m->control_pressed) { break; } + for (int g = 0; g < theseFlowFileNames.size(); g++) { - calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber); - - if (m->control_pressed) { break; } + string flowFileName = theseFlowFileNames[g]; + m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(filenames.size()) + ")\t<<<<<\n"); + m->mothurOut("Reading flowgrams...\n"); + + vector seqNameVector; + vector lengths; + vector flowDataIntI; + vector flowDataPrI; + map nameMap; + vector uniqueFlowgrams; + vector uniqueCount; + vector mapSeqToUnique; + vector mapUniqueToSeq; + vector uniqueLengths; + int numFlowCells; - maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight); + int numSeqs = getFlowData(flowFileName, seqNameVector, lengths, flowDataIntI, nameMap, numFlowCells); if (m->control_pressed) { break; } - double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight); + m->mothurOut("Identifying unique flowgrams...\n"); + int numUniques = getUniques(numSeqs, numFlowCells, uniqueFlowgrams, uniqueCount, uniqueLengths, mapSeqToUnique, mapUniqueToSeq, lengths, flowDataPrI, flowDataIntI); if (m->control_pressed) { break; } - checkCentroids(numOTUs, centroids, weight); - - if (m->control_pressed) { break; } - - calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths); - - if (m->control_pressed) { break; } - - iter++; - - m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n'); + m->mothurOut("Calculating distances between flowgrams...\n"); + string distFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.dist"; + begTime = time(NULL); + - } - - if (m->control_pressed) { break; } - - m->mothurOut("\nFinalizing...\n"); - fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); - - if (m->control_pressed) { break; } - - setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI); - - if (m->control_pressed) { break; } - - vector otuCounts(numOTUs, 0); - for(int i=0;icontrol_pressed) { break; } + flowDistParentFork(numFlowCells, distFileName, numUniques, mapUniqueToSeq, mapSeqToUnique, lengths, flowDataPrI, flowDataIntI); + + m->mothurOutEndLine(); + m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n'); + + + string namesFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.names"; + createNamesFile(numSeqs, numUniques, namesFileName, seqNameVector, mapSeqToUnique, mapUniqueToSeq); + + if (m->control_pressed) { break; } + + m->mothurOut("\nClustering flowgrams...\n"); + string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".shhh.list"; + cluster(listFileName, distFileName, namesFileName); + + if (m->control_pressed) { break; } + + vector otuData; + vector cumNumSeqs; + vector nSeqsPerOTU; + vector > aaP; //tMaster->aanP: each row is a different otu / each col contains the sequence indices + vector > aaI; //tMaster->aanI: that are in each otu - can't differentiate between aaP and aaI + vector seqNumber; //tMaster->anP: the sequence id number sorted by OTU + vector seqIndex; //tMaster->anI; the index that corresponds to seqNumber + + + int numOTUs = getOTUData(numSeqs, listFileName, otuData, cumNumSeqs, nSeqsPerOTU, aaP, aaI, seqNumber, seqIndex, nameMap); + + if (m->control_pressed) { break; } + + m->mothurRemove(distFileName); + m->mothurRemove(namesFileName); + m->mothurRemove(listFileName); + + vector dist; //adDist - distance of sequences to centroids + vector change; //did the centroid sequence change? 0 = no; 1 = yes + vector centroids; //the representative flowgram for each cluster m + vector weight; + vector singleTau; //tMaster->adTau: 1-D Tau vector (1xnumSeqs) + vector nSeqsBreaks; + vector nOTUsBreaks; + + dist.assign(numSeqs * numOTUs, 0); + change.assign(numOTUs, 1); + centroids.assign(numOTUs, -1); + weight.assign(numOTUs, 0); + singleTau.assign(numSeqs, 1.0); + + nSeqsBreaks.assign(2, 0); + nOTUsBreaks.assign(2, 0); + + nSeqsBreaks[0] = 0; + nSeqsBreaks[1] = numSeqs; + nOTUsBreaks[1] = numOTUs; + + if (m->control_pressed) { break; } + + double maxDelta = 0; + int iter = 0; + + begClock = clock(); + begTime = time(NULL); + + m->mothurOut("\nDenoising flowgrams...\n"); + m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n"); + + while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){ + + if (m->control_pressed) { break; } + + double cycClock = clock(); + unsigned long long cycTime = time(NULL); + fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); + + if (m->control_pressed) { break; } + + calcCentroidsDriver(numOTUs, cumNumSeqs, nSeqsPerOTU, seqIndex, change, centroids, singleTau, mapSeqToUnique, uniqueFlowgrams, flowDataIntI, lengths, numFlowCells, seqNumber); + + if (m->control_pressed) { break; } + + maxDelta = getNewWeights(numOTUs, cumNumSeqs, nSeqsPerOTU, singleTau, seqNumber, weight); + + if (m->control_pressed) { break; } + + double nLL = getLikelihood(numSeqs, numOTUs, nSeqsPerOTU, seqNumber, cumNumSeqs, seqIndex, dist, weight); + + if (m->control_pressed) { break; } + + checkCentroids(numOTUs, centroids, weight); + + if (m->control_pressed) { break; } + + calcNewDistances(numSeqs, numOTUs, nSeqsPerOTU, dist, weight, change, centroids, aaP, singleTau, aaI, seqNumber, seqIndex, uniqueFlowgrams, flowDataIntI, numFlowCells, lengths); + + if (m->control_pressed) { break; } + + iter++; + + m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n'); + + } + + if (m->control_pressed) { break; } + + m->mothurOut("\nFinalizing...\n"); + fill(numOTUs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, aaP, aaI); + + if (m->control_pressed) { break; } + + setOTUs(numOTUs, numSeqs, seqNumber, seqIndex, cumNumSeqs, nSeqsPerOTU, otuData, singleTau, dist, aaP, aaI); + + if (m->control_pressed) { break; } + + vector otuCounts(numOTUs, 0); + for(int i=0;icontrol_pressed) { break; } + + if ((large) && (g == 0)) { flowFileName = filenames[i]; theseFlowFileNames[0] = filenames[i]; } + string thisOutputDir = outputDir; + if (outputDir == "") { thisOutputDir = m->hasPath(flowFileName); } + string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.qual"; + string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.fasta"; + string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.names"; + string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)) + "shhh.counts"; + string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(flowFileName)); + string groupFileName = fileRoot + "shhh.groups"; + + + writeQualities(numOTUs, numFlowCells, qualityFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; } + writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, fastaFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; } + writeNames(thisCompositeNamesFileName, numOTUs, nameFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; } + writeClusters(otuCountsFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; } + writeGroups(groupFileName, fileRoot, numSeqs, seqNameVector); if (m->control_pressed) { break; } + + if (large) { + if (g > 0) { + m->appendFiles(qualityFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.qual")); + m->mothurRemove(qualityFileName); + m->appendFiles(fastaFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.fasta")); + m->mothurRemove(fastaFileName); + m->appendFiles(nameFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.names")); + m->mothurRemove(nameFileName); + m->appendFiles(otuCountsFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.counts")); + m->mothurRemove(otuCountsFileName); + m->appendFiles(groupFileName, (thisOutputDir + m->getRootName(m->getSimpleName(theseFlowFileNames[0])) + "shhh.groups")); + m->mothurRemove(groupFileName); + m->mothurRemove(theseFlowFileNames[g]); + } + } + } - writeQualities(numOTUs, numFlowCells, flowFileName, otuCounts, nSeqsPerOTU, seqNumber, singleTau, flowDataIntI, uniqueFlowgrams, cumNumSeqs, mapUniqueToSeq, seqNameVector, centroids, aaI); if (m->control_pressed) { break; } - writeSequences(thisCompositeFASTAFileName, numOTUs, numFlowCells, flowFileName, otuCounts, uniqueFlowgrams, seqNameVector, aaI, centroids);if (m->control_pressed) { break; } - writeNames(thisCompositeNamesFileName, numOTUs, flowFileName, otuCounts, seqNameVector, aaI, nSeqsPerOTU); if (m->control_pressed) { break; } - writeClusters(flowFileName, numOTUs, numFlowCells,otuCounts, centroids, uniqueFlowgrams, seqNameVector, aaI, nSeqsPerOTU, lengths, flowDataIntI); if (m->control_pressed) { break; } - writeGroups(flowFileName, numSeqs, seqNameVector); if (m->control_pressed) { break; } - m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n'); } @@ -2977,14 +3067,11 @@ void ShhherCommand::setOTUs(int numOTUs, int numSeqs, vector& seqNumber, ve } /**************************************************************************************************/ -void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string filename, vector otuCounts, vector& nSeqsPerOTU, vector& seqNumber, +void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string qualityFileName, vector otuCounts, vector& nSeqsPerOTU, vector& seqNumber, vector& singleTau, vector& flowDataIntI, vector& uniqueFlowgrams, vector& cumNumSeqs, vector& mapUniqueToSeq, vector& seqNameVector, vector& centroids, vector >& aaI){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir = m->hasPath(filename); } - string qualityFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.qual"; ofstream qualityFile; m->openOutputFile(qualityFileName, qualityFile); @@ -3088,11 +3175,9 @@ void ShhherCommand::writeQualities(int numOTUs, int numFlowCells, string filenam /**************************************************************************************************/ -void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTUs, int numFlowCells, string filename, vector otuCounts, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& centroids){ +void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTUs, int numFlowCells, string fastaFileName, vector otuCounts, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& centroids){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir = m->hasPath(filename); } - string fastaFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.fasta"; + ofstream fastaFile; m->openOutputFile(fastaFileName, fastaFile); @@ -3117,7 +3202,8 @@ void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTU } } - fastaFile << newSeq.substr(4) << endl; + if (newSeq.length() >= 4) { fastaFile << newSeq.substr(4) << endl; } + else { fastaFile << "NNNN" << endl; } } } fastaFile.close(); @@ -3136,11 +3222,9 @@ void ShhherCommand::writeSequences(string thisCompositeFASTAFileName, int numOTU /**************************************************************************************************/ -void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string filename, vector otuCounts, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU){ +void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, string nameFileName, vector otuCounts, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir = m->hasPath(filename); } - string nameFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.names"; + ofstream nameFile; m->openOutputFile(nameFileName, nameFile); @@ -3174,13 +3258,9 @@ void ShhherCommand::writeNames(string thisCompositeNamesFileName, int numOTUs, s /**************************************************************************************************/ -void ShhherCommand::writeGroups(string filename, int numSeqs, vector& seqNameVector){ +void ShhherCommand::writeGroups(string groupFileName, string fileRoot, int numSeqs, vector& seqNameVector){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir = m->hasPath(filename); } - string fileRoot = thisOutputDir + m->getRootName(m->getSimpleName(filename)); - string groupFileName = fileRoot + "shhh.groups"; - ofstream groupFile; + ofstream groupFile; m->openOutputFile(groupFileName, groupFile); for(int i=0;i& se /**************************************************************************************************/ -void ShhherCommand::writeClusters(string filename, int numOTUs, int numFlowCells, vector otuCounts, vector& centroids, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU, vector& lengths, vector& flowDataIntI){ +void ShhherCommand::writeClusters(string otuCountsFileName, int numOTUs, int numFlowCells, vector otuCounts, vector& centroids, vector& uniqueFlowgrams, vector& seqNameVector, vector >& aaI, vector& nSeqsPerOTU, vector& lengths, vector& flowDataIntI){ try { - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir = m->hasPath(filename); } - string otuCountsFileName = thisOutputDir + m->getRootName(m->getSimpleName(filename)) + "shhh.counts"; ofstream otuCountsFile; m->openOutputFile(otuCountsFileName, otuCountsFile); diff --git a/shhhercommand.h b/shhhercommand.h index c9772af..5ef43df 100644 --- a/shhhercommand.h +++ b/shhhercommand.h @@ -56,10 +56,10 @@ private: linePair(int i, int j) : start(i), end(j) {} }; - int abort; + bool abort, large; string outputDir, flowFileName, flowFilesFileName, lookupFileName, compositeFASTAFileName, compositeNamesFileName; - int processors, maxIters; + int processors, maxIters, largeSize; float cutoff, sigma, minDelta; string flowOrder; @@ -68,6 +68,7 @@ private: vector jointLookUp; vector flowFileVector; + vector parseFlowFiles(string); int driver(vector, string, string, int, int); int createProcesses(vector); int getFlowData(string, vector&, vector&, vector&, map&, int&); @@ -90,7 +91,7 @@ private: void writeQualities(int, int, string, vector, vector&, vector&, vector&, vector&, vector&, vector&, vector&, vector&, vector&, vector >&); void writeSequences(string, int, int, string, vector, vector&, vector&, vector >&, vector&); void writeNames(string, int, string, vector, vector&, vector >&, vector&); - void writeGroups(string, int, vector&); + void writeGroups(string, string, int, vector&); void writeClusters(string, int, int, vector, vector&, vector&, vector&, vector >&, vector&, vector&, vector&); void getSingleLookUp(); -- 2.39.2