X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=seqerrorcommand.cpp;h=5ec6cf2d1cc4778399d289ea2280c882f936b944;hb=16f9c4ab6f39769856b13e048eae2c8eaa413c02;hp=253e03e2a1d380cb4d13f360fc9a5fecaaeb0808;hpb=8977d23b7f18d981b14e59c210c97198541ca7b0;p=mothur.git diff --git a/seqerrorcommand.cpp b/seqerrorcommand.cpp index 253e03e..5ec6cf2 100644 --- a/seqerrorcommand.cpp +++ b/seqerrorcommand.cpp @@ -63,20 +63,52 @@ string SeqErrorCommand::getHelpString(){ } } //********************************************************************************************************************** +string SeqErrorCommand::getOutputFileNameTag(string type, string inputName=""){ + try { + string outputFileName = ""; + map >::iterator it; + + //is this a type this command creates + it = outputTypes.find(type); + if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); } + else { + if (type == "errorsummary") { outputFileName = "error.summary"; } + else if (type == "errorseq") { outputFileName = "error.seq"; } + else if (type == "errorquality") { outputFileName = "error.quality"; } + else if (type == "errorqualforward") { outputFileName = "error.qual.forward"; } + else if (type == "errorqualreverse") { outputFileName = "error.qual.reverse"; } + else if (type == "errorforward") { outputFileName = "error.seq.forward"; } + else if (type == "errorreverse") { outputFileName = "error.seq.reverse"; } + else if (type == "errorcount") { outputFileName = "error.count"; } + else if (type == "errormatrix") { outputFileName = "error.matrix"; } + else if (type == "errorchimera") { outputFileName = "error.chimera"; } + else if (type == "errorref-query") { outputFileName = "error.ref-query"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; } + } + return outputFileName; + } + catch(exception& e) { + m->errorOut(e, "SeqErrorCommand", "getOutputFileNameTag"); + exit(1); + } +} +//********************************************************************************************************************** SeqErrorCommand::SeqErrorCommand(){ try { abort = true; calledHelp = true; setParameters(); vector tempOutNames; - outputTypes["error.summary"] = tempOutNames; - outputTypes["error.seq"] = tempOutNames; - outputTypes["error.quality"] = tempOutNames; - outputTypes["error.qual.forward"] = tempOutNames; - outputTypes["error.qual.reverse"] = tempOutNames; - outputTypes["error.forward"] = tempOutNames; - outputTypes["error.reverse"] = tempOutNames; - outputTypes["error.count"] = tempOutNames; - outputTypes["error.matrix"] = tempOutNames; + outputTypes["errorsummary"] = tempOutNames; + outputTypes["errorseq"] = tempOutNames; + outputTypes["errorquality"] = tempOutNames; + outputTypes["errorqualforward"] = tempOutNames; + outputTypes["errorqualreverse"] = tempOutNames; + outputTypes["errorforward"] = tempOutNames; + outputTypes["errorreverse"] = tempOutNames; + outputTypes["errorcount"] = tempOutNames; + outputTypes["errormatrix"] = tempOutNames; + outputTypes["errorchimera"] = tempOutNames; + outputTypes["errorref-query"] = tempOutNames; } catch(exception& e) { m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand"); @@ -112,15 +144,17 @@ SeqErrorCommand::SeqErrorCommand(string option) { //initialize outputTypes vector tempOutNames; - outputTypes["error.summary"] = tempOutNames; - outputTypes["error.seq"] = tempOutNames; - outputTypes["error.quality"] = tempOutNames; - outputTypes["error.qual.forward"] = tempOutNames; - outputTypes["error.qual.reverse"] = tempOutNames; - outputTypes["error.forward"] = tempOutNames; - outputTypes["error.reverse"] = tempOutNames; - outputTypes["error.count"] = tempOutNames; - outputTypes["error.matrix"] = tempOutNames; + outputTypes["errorsummary"] = tempOutNames; + outputTypes["errorseq"] = tempOutNames; + outputTypes["errorquality"] = tempOutNames; + outputTypes["errorqualforward"] = tempOutNames; + outputTypes["errorqualreverse"] = tempOutNames; + outputTypes["errorforward"] = tempOutNames; + outputTypes["errorreverse"] = tempOutNames; + outputTypes["errorcount"] = tempOutNames; + outputTypes["errormatrix"] = tempOutNames; + outputTypes["errorchimera"] = tempOutNames; + outputTypes["errorref-query"] = tempOutNames; //if the user changes the input directory command factory will send this info to us in the output parameter @@ -176,7 +210,7 @@ SeqErrorCommand::SeqErrorCommand(string option) { if (queryFileName != "") { m->mothurOut("Using " + queryFileName + " as input file for the fasta parameter."); m->mothurOutEndLine(); } else { m->mothurOut("You have no current fasta file and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; } } - else if (queryFileName == "not open") { abort = true; } + else if (queryFileName == "not open") { queryFileName = ""; abort = true; } else { m->setFastaFile(queryFileName); } referenceFileName = validParameter.validFile(parameters, "reference", true); @@ -213,7 +247,7 @@ SeqErrorCommand::SeqErrorCommand(string option) { //check for optional parameter and set defaults // ...at some point should added some additional type checking... temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; } - convert(temp, threshold); + m->mothurConvert(temp, threshold); temp = validParameter.validFile(parameters, "save", false); if (temp == "not found"){ temp = "f"; } save = m->isTrue(temp); @@ -242,10 +276,15 @@ SeqErrorCommand::SeqErrorCommand(string option) { temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } m->setProcessors(temp); - convert(temp, processors); + m->mothurConvert(temp, processors); substitutionMatrix.resize(6); for(int i=0;i<6;i++){ substitutionMatrix[i].resize(6,0); } + + if ((namesFileName == "") && (queryFileName != "")){ + vector files; files.push_back(queryFileName); + parser.getNameFile(files); + } } } catch(exception& e) { @@ -264,14 +303,15 @@ int SeqErrorCommand::execute(){ totalBases = 0; totalMatches = 0; - string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary"; - outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName); + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + string errorSummaryFileName = fileNameRoot + getOutputFileNameTag("errorsummary"); + outputNames.push_back(errorSummaryFileName); outputTypes["errorsummary"].push_back(errorSummaryFileName); - string errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq"; - outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName); + string errorSeqFileName = fileNameRoot + getOutputFileNameTag("errorseq"); + outputNames.push_back(errorSeqFileName); outputTypes["errorseq"].push_back(errorSeqFileName); - string errorChimeraFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.chimera"; - outputNames.push_back(errorChimeraFileName); outputTypes["error.chimera"].push_back(errorChimeraFileName); + string errorChimeraFileName = fileNameRoot + getOutputFileNameTag("errorchimera"); + outputNames.push_back(errorChimeraFileName); outputTypes["errorchimera"].push_back(errorChimeraFileName); getReferences(); //read in reference sequences - make sure there's no ambiguous bases @@ -293,7 +333,7 @@ int SeqErrorCommand::execute(){ if(qualFileName == "") { qLines = lines; rLines = lines; } //fills with duds int numSeqs = 0; -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) if(processors == 1){ numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]); }else{ @@ -312,10 +352,10 @@ int SeqErrorCommand::execute(){ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } - string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count"; + string errorCountFileName = fileNameRoot + getOutputFileNameTag("errorcount"); ofstream errorCountFile; m->openOutputFile(errorCountFileName, errorCountFile); - outputNames.push_back(errorCountFileName); outputTypes["error.count"].push_back(errorCountFileName); + outputNames.push_back(errorCountFileName); outputTypes["errorcount"].push_back(errorCountFileName); m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n"); m->mothurOut("Errors\tSequences\n"); errorCountFile << "Errors\tSequences\n"; @@ -329,10 +369,10 @@ int SeqErrorCommand::execute(){ printSubMatrix(); - string megAlignmentFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.ref-query"; + string megAlignmentFileName = fileNameRoot + getOutputFileNameTag("errorref-query"); ofstream megAlignmentFile; m->openOutputFile(megAlignmentFileName, megAlignmentFile); - outputNames.push_back(megAlignmentFileName); outputTypes["error.ref-query"].push_back(megAlignmentFileName); + outputNames.push_back(megAlignmentFileName); outputTypes["errorref-query"].push_back(megAlignmentFileName); for(int i=0;i >::iterator it; int num = 0; -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) //loop through and create all the processes you want while (process != processors) { @@ -677,7 +717,8 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, if(numParentSeqs > 1 && ignoreChimeras == 1) { ignoreSeq = 1; } else { ignoreSeq = 0; } - Compare minCompare = getErrors(query, referenceSeqs[closestRefIndex]); + Compare minCompare; + getErrors(query, referenceSeqs[closestRefIndex], minCompare); if(namesFileName != ""){ it = weights.find(query.getName()); @@ -734,14 +775,14 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, index++; - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) unsigned long long pos = queryFile.tellg(); if ((pos == -1) || (pos >= line.end)) { break; } #else if (queryFile.eof()) { break; } #endif - if(index % 100 == 0){ m->mothurOut(toString(index) + '\n'); } + if(index % 100 == 0){ m->mothurOut(toString(index)); m->mothurOutEndLine(); } } queryFile.close(); if(qFileName != "" && rFileName != ""){ reportFile.close(); qualFile.close(); } @@ -749,7 +790,7 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, errorSeqFile.close(); //report progress - if(index % 100 != 0){ m->mothurOut(toString(index) + '\n'); } + if(index % 100 != 0){ m->mothurOut(toString(index)); m->mothurOutEndLine(); } return index; } @@ -780,12 +821,16 @@ void SeqErrorCommand::getReferences(){ // // int endPos = rdb->referenceSeqs[i].getEndPos(); // if(endPos < minEndPos) { minEndPos = endPos; } + if (rdb->referenceSeqs[i].getNumBases() == 0) { + m->mothurOut("[WARNING]: " + rdb->referenceSeqs[i].getName() + " is blank, ignoring.");m->mothurOutEndLine(); + }else { + referenceSeqs.push_back(rdb->referenceSeqs[i]); + } - referenceSeqs.push_back(rdb->referenceSeqs[i]); } referenceFileName = rdb->getSavedReference(); - m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences.");m->mothurOutEndLine(); + m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(referenceSeqs.size()) + " sequences.");m->mothurOutEndLine(); }else { int start = time(NULL); @@ -803,9 +848,12 @@ void SeqErrorCommand::getReferences(){ // // int endPos = currentSeq.getEndPos(); // if(endPos < minEndPos) { minEndPos = endPos; } - referenceSeqs.push_back(currentSeq); - - if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); } + if (currentSeq.getNumBases() == 0) { + m->mothurOut("[WARNING]: " + currentSeq.getName() + " is blank, ignoring.");m->mothurOutEndLine(); + }else { + referenceSeqs.push_back(currentSeq); + if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); } + } m->gobble(referenceFile); } @@ -819,7 +867,7 @@ void SeqErrorCommand::getReferences(){ for(int i=0;imothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n"); @@ -834,7 +882,7 @@ void SeqErrorCommand::getReferences(){ //*************************************************************************************************************** -Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){ +int SeqErrorCommand::getErrors(Sequence query, Sequence reference, Compare& errors){ try { if(query.getAlignLength() != reference.getAlignLength()){ m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n"); @@ -845,7 +893,7 @@ Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){ string r = reference.getAligned(); int started = 0; - Compare errors; + //Compare errors; for(int i=0;ierrorOut(e, "SeqErrorCommand", "getErrors"); @@ -1039,10 +1087,11 @@ void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs, ofstream& void SeqErrorCommand::printSubMatrix(){ try { - string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + string subMatrixFileName = fileNameRoot + getOutputFileNameTag("errormatrix"); ofstream subMatrixFile; m->openOutputFile(subMatrixFileName, subMatrixFile); - outputNames.push_back(subMatrixFileName); outputTypes["error.matrix"].push_back(subMatrixFileName); + outputNames.push_back(subMatrixFileName); outputTypes["errormatrix"].push_back(subMatrixFileName); vector bases(6); bases[0] = "A"; bases[1] = "T"; @@ -1085,10 +1134,11 @@ void SeqErrorCommand::printSubMatrix(){ void SeqErrorCommand::printErrorFRFile(map > errorForward, map > errorReverse){ try{ - string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + string errorForwardFileName = fileNameRoot + getOutputFileNameTag("errorforward"); ofstream errorForwardFile; m->openOutputFile(errorForwardFileName, errorForwardFile); - outputNames.push_back(errorForwardFileName); outputTypes["error.forward"].push_back(errorForwardFileName); + outputNames.push_back(errorForwardFileName); outputTypes["errorforward"].push_back(errorForwardFileName); errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl; for(int i=0;i > errorForward, map } errorForwardFile.close(); - string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse"; + string errorReverseFileName = fileNameRoot + getOutputFileNameTag("errorreverse"); ofstream errorReverseFile; m->openOutputFile(errorReverseFileName, errorReverseFile); - outputNames.push_back(errorReverseFileName); outputTypes["error.reverse"].push_back(errorReverseFileName); + outputNames.push_back(errorReverseFileName); outputTypes["errorreverse"].push_back(errorReverseFileName); errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl; for(int i=0;i > errorForward, map void SeqErrorCommand::printErrorQuality(map > qScoreErrorMap){ try{ - - string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + string errorQualityFileName = fileNameRoot + getOutputFileNameTag("errorquality"); ofstream errorQualityFile; m->openOutputFile(errorQualityFileName, errorQualityFile); - outputNames.push_back(errorQualityFileName); outputTypes["error.quality"].push_back(errorQualityFileName); + outputNames.push_back(errorQualityFileName); outputTypes["errorquality"].push_back(errorQualityFileName); errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl; for(int i=0;i<41;i++){ @@ -1165,11 +1215,11 @@ void SeqErrorCommand::printQualityFR(vector > qualForwardMap, vector } } } - - string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + string qualityForwardFileName = fileNameRoot + getOutputFileNameTag("errorqualforward"); ofstream qualityForwardFile; m->openOutputFile(qualityForwardFileName, qualityForwardFile); - outputNames.push_back(qualityForwardFileName); outputTypes["error.qual.forward"].push_back(qualityForwardFileName); + outputNames.push_back(qualityForwardFileName); outputTypes["errorqualforward"].push_back(qualityForwardFileName); for(int i=0;i > qualForwardMap, vector qualityForwardFile.close(); - string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse"; + string qualityReverseFileName = fileNameRoot + getOutputFileNameTag("errorqualreverse"); ofstream qualityReverseFile; m->openOutputFile(qualityReverseFileName, qualityReverseFile); - outputNames.push_back(qualityReverseFileName); outputTypes["error.qual.reverse"].push_back(qualityReverseFileName); + outputNames.push_back(qualityReverseFileName); outputTypes["errorqualreverse"].push_back(qualityReverseFileName); for(int i=0;i > qualForwardMap, vector int SeqErrorCommand::setLines(string filename, string qfilename, string rfilename, vector& fastaFilePos, vector& qfileFilePos, vector& rfileFilePos) { try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) //set file positions for fasta file fastaFilePos = m->divideFile(filename, processors);