X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=trimseqscommand.cpp;h=d72ada4cff4fb44cd557eede88d63f44a43c1cd6;hb=5c80ce8b80938d41cf6c64a017fa6fd50d45de5b;hp=b185da0905321c7980f2237a9d938bd260d9397d;hpb=c99f3846e7a7b6f06ab46508baa5409204ad6290;p=mothur.git diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp index b185da0..d72ada4 100644 --- a/trimseqscommand.cpp +++ b/trimseqscommand.cpp @@ -8,10 +8,12 @@ */ #include "trimseqscommand.h" +#include "needlemanoverlap.hpp" +#include "nast.hpp" //*************************************************************************************************************** -TrimSeqsCommand::TrimSeqsCommand(string option){ +TrimSeqsCommand::TrimSeqsCommand(string option) { try { abort = false; @@ -21,7 +23,8 @@ TrimSeqsCommand::TrimSeqsCommand(string option){ else { //valid paramters for this command - string AlignArray[] = {"fasta", "flip", "oligos", "maxambig", "maxhomop", "minlength", "maxlength", "qfile", "qthreshold", "qaverage", "allfiles", "qtrim"}; + string AlignArray[] = {"fasta", "flip", "oligos", "maxambig", "maxhomop", "minlength", "maxlength", "qfile", + "qthreshold", "qaverage", "allfiles", "qtrim","diffs", "processors", "outputdir","inputdir"}; vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); @@ -29,17 +32,54 @@ TrimSeqsCommand::TrimSeqsCommand(string option){ map parameters = parser.getParameters(); ValidParameters validParameter; + map::iterator it; //check to make sure all parameters are valid for command - for (map::iterator it = parameters.begin(); it != parameters.end(); it++) { + for (it = parameters.begin(); it != parameters.end(); it++) { if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; } } + //if the user changes the input directory command factory will send this info to us in the output parameter + string inputDir = validParameter.validFile(parameters, "inputdir", false); + if (inputDir == "not found"){ inputDir = ""; } + else { + string path; + 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; } + } + + it = parameters.find("oligos"); + //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["oligos"] = inputDir + it->second; } + } + + it = parameters.find("qfile"); + //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["qfile"] = inputDir + it->second; } + } + } + + //check for required parameters fastaFile = validParameter.validFile(parameters, "fasta", true); - if (fastaFile == "not found") { mothurOut("fasta is a required parameter for the screen.seqs command."); mothurOutEndLine(); abort = true; } + if (fastaFile == "not found") { m->mothurOut("fasta is a required parameter for the screen.seqs command."); m->mothurOutEndLine(); abort = true; } else if (fastaFile == "not open") { abort = true; } - + + //if the user changes the output directory command factory will send this info to us in the output parameter + outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ + outputDir = ""; + outputDir += hasPath(fastaFile); //if user entered a file with a path then preserve it + } //check for optional parameter and set defaults // ...at some point should added some additional type checking... @@ -65,9 +105,12 @@ TrimSeqsCommand::TrimSeqsCommand(string option){ temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found") { temp = "0"; } convert(temp, maxLength); + temp = validParameter.validFile(parameters, "diffs", false); if (temp == "not found") { temp = "0"; } + convert(temp, diffs); + temp = validParameter.validFile(parameters, "qfile", true); if (temp == "not found") { qFileName = ""; } - else if(temp == "not open") { abort = 0; } + else if(temp == "not open") { abort = true; } else { qFileName = temp; } temp = validParameter.validFile(parameters, "qthreshold", false); if (temp == "not found") { temp = "0"; } @@ -82,23 +125,26 @@ TrimSeqsCommand::TrimSeqsCommand(string option){ temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found") { temp = "F"; } allFiles = isTrue(temp); + temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; } + convert(temp, processors); + if(allFiles && oligoFile == ""){ - mothurOut("You selected allfiles, but didn't enter an oligos file. Ignoring the allfiles request."); mothurOutEndLine(); + m->mothurOut("You selected allfiles, but didn't enter an oligos file. Ignoring the allfiles request."); m->mothurOutEndLine(); } if((qAverage != 0 && qThreshold != 0) && qFileName == ""){ - mothurOut("You didn't provide a quality file name, quality criteria will be ignored."); mothurOutEndLine(); + m->mothurOut("You didn't provide a quality file name, quality criteria will be ignored."); m->mothurOutEndLine(); qAverage=0; qThreshold=0; } if(!flip && oligoFile=="" && !maxLength && !minLength && (maxAmbig==-1) && !maxHomoP && qFileName == ""){ - mothurOut("You didn't set any options... quiting command."); mothurOutEndLine(); + m->mothurOut("You didn't set any options... quiting command."); m->mothurOutEndLine(); abort = true; } } } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "TrimSeqsCommand"); + m->errorOut(e, "TrimSeqsCommand", "TrimSeqsCommand"); exit(1); } } @@ -106,30 +152,31 @@ TrimSeqsCommand::TrimSeqsCommand(string option){ void TrimSeqsCommand::help(){ try { - mothurOut("The trim.seqs command reads a fastaFile and creates .....\n"); - mothurOut("The trim.seqs command parameters are fasta, flip, oligos, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, qtrim and allfiles.\n"); - mothurOut("The fasta parameter is required.\n"); - mothurOut("The flip parameter .... The default is 0.\n"); - mothurOut("The oligos parameter .... The default is "".\n"); - mothurOut("The maxambig parameter .... The default is -1.\n"); - mothurOut("The maxhomop parameter .... The default is 0.\n"); - mothurOut("The minlength parameter .... The default is 0.\n"); - mothurOut("The maxlength parameter .... The default is 0.\n"); - mothurOut("The qfile parameter .....\n"); - mothurOut("The qthreshold parameter .... The default is 0.\n"); - mothurOut("The qaverage parameter .... The default is 0.\n"); - mothurOut("The allfiles parameter .... The default is F.\n"); - mothurOut("The qtrim parameter .... The default is F.\n"); - mothurOut("The trim.seqs command should be in the following format: \n"); - mothurOut("trim.seqs(fasta=yourFastaFile, flip=yourFlip, oligos=yourOligos, maxambig=yourMaxambig, \n"); - mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n"); - mothurOut("Example trim.seqs(fasta=abrecovery.fasta, flip=..., oligos=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n"); - mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n"); - mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/Trim.seqs .\n\n"); + m->mothurOut("The trim.seqs command reads a fastaFile and creates .....\n"); + m->mothurOut("The trim.seqs command parameters are fasta, flip, oligos, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, diffs, qtrim and allfiles.\n"); + m->mothurOut("The fasta parameter is required.\n"); + m->mothurOut("The flip parameter .... The default is 0.\n"); + m->mothurOut("The oligos parameter .... The default is "".\n"); + m->mothurOut("The maxambig parameter .... The default is -1.\n"); + m->mothurOut("The maxhomop parameter .... The default is 0.\n"); + m->mothurOut("The minlength parameter .... The default is 0.\n"); + m->mothurOut("The maxlength parameter .... The default is 0.\n"); + m->mothurOut("The diffs parameter .... The default is 0.\n"); + m->mothurOut("The qfile parameter .....\n"); + m->mothurOut("The qthreshold parameter .... The default is 0.\n"); + m->mothurOut("The qaverage parameter .... The default is 0.\n"); + m->mothurOut("The allfiles parameter .... The default is F.\n"); + m->mothurOut("The qtrim parameter .... The default is F.\n"); + m->mothurOut("The trim.seqs command should be in the following format: \n"); + m->mothurOut("trim.seqs(fasta=yourFastaFile, flip=yourFlip, oligos=yourOligos, maxambig=yourMaxambig, \n"); + m->mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n"); + m->mothurOut("Example trim.seqs(fasta=abrecovery.fasta, flip=..., oligos=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n"); + m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n"); + m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/Trim.seqs .\n\n"); } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "help"); + m->errorOut(e, "TrimSeqsCommand", "help"); exit(1); } } @@ -145,33 +192,168 @@ int TrimSeqsCommand::execute(){ try{ if (abort == true) { return 0; } + + numFPrimers = 0; //this needs to be initialized + numRPrimers = 0; + + string trimSeqFile = outputDir + getRootName(getSimpleName(fastaFile)) + "trim.fasta"; + outputNames.push_back(trimSeqFile); + string scrapSeqFile = outputDir + getRootName(getSimpleName(fastaFile)) + "scrap.fasta"; + outputNames.push_back(scrapSeqFile); + string groupFile = outputDir + getRootName(getSimpleName(fastaFile)) + "groups"; + + vector fastaFileNames; + if(oligoFile != ""){ + outputNames.push_back(groupFile); + getOligos(fastaFileNames); + } + + if(qFileName != "") { setLines(qFileName, qLines); } - ifstream inFASTA; - openInputFile(fastaFile, inFASTA); + + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + if(processors == 1){ + ifstream inFASTA; + openInputFile(fastaFile, inFASTA); + int numSeqs=count(istreambuf_iterator(inFASTA),istreambuf_iterator(), '>'); + inFASTA.close(); + + lines.push_back(new linePair(0, numSeqs)); + + driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, groupFile, fastaFileNames, lines[0], lines[0]); + + for (int j = 0; j < fastaFileNames.size(); j++) { + rename((fastaFileNames[j] + toString(getpid()) + ".temp").c_str(), fastaFileNames[j].c_str()); + } + + }else{ + setLines(fastaFile, lines); + if(qFileName == "") { qLines = lines; } + + createProcessesCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, groupFile, fastaFileNames); + + rename((trimSeqFile + toString(processIDS[0]) + ".temp").c_str(), trimSeqFile.c_str()); + rename((scrapSeqFile + toString(processIDS[0]) + ".temp").c_str(), scrapSeqFile.c_str()); + rename((groupFile + toString(processIDS[0]) + ".temp").c_str(), groupFile.c_str()); + for (int j = 0; j < fastaFileNames.size(); j++) { + rename((fastaFileNames[j] + toString(processIDS[0]) + ".temp").c_str(), fastaFileNames[j].c_str()); + } + //append files + for(int i=1;icontrol_pressed) { return 0; } + #else + ifstream inFASTA; + openInputFile(fastafileNames[s], inFASTA); + numSeqs=count(istreambuf_iterator(inFASTA),istreambuf_iterator(), '>'); + inFASTA.close(); + + lines.push_back(new linePair(0, numSeqs)); + + driverCreateSummary(fastafile, qFileName, trimSeqFile, scrapSeqFile, groupFile, fastaFileNames, lines[0], lines[0]); + + if (m->control_pressed) { return 0; } + #endif + + + for(int i=0;i'){ + inFASTA >> seqName; + outGroups << seqName << '\t' << groupVector[i] << endl; + } + while (!inFASTA.eof()) { char c = inFASTA.get(); if (c == 10 || c == 13){ break; } } + } + outGroups.close(); + inFASTA.close(); + } + + if (m->control_pressed) { + for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } + return 0; + } + + m->mothurOutEndLine(); + m->mothurOut("Output File Names: "); m->mothurOutEndLine(); + for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } + m->mothurOutEndLine(); + + return 0; + + } + catch(exception& e) { + m->errorOut(e, "TrimSeqsCommand", "execute"); + exit(1); + } +} + +/**************************************************************************************/ +int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string trimFile, string scrapFile, string groupFile, vector fastaNames, linePair* line, linePair* qline) { + try { ofstream outFASTA; - string trimSeqFile = getRootName(fastaFile) + "trim.fasta"; - openOutputFile(trimSeqFile, outFASTA); + int able = openOutputFile(trimFile, outFASTA); + + ofstream scrapFASTA; + openOutputFile(scrapFile, scrapFASTA); ofstream outGroups; vector fastaFileNames; - if(oligoFile != ""){ - string groupFile = getRootName(fastaFile) + "groups"; - openOutputFile(groupFile, outGroups); - getOligos(fastaFileNames); + if (oligoFile != "") { + openOutputFile(groupFile, outGroups); + for (int i = 0; i < fastaNames.size(); i++) { + fastaFileNames.push_back(new ofstream((fastaNames[i] + toString(getpid()) + ".temp").c_str(), ios::ate)); + } } - ofstream scrapFASTA; - string scrapSeqFile = getRootName(fastaFile) + "scrap.fasta"; - openOutputFile(scrapSeqFile, scrapFASTA); + ifstream inFASTA; + openInputFile(filename, inFASTA); ifstream qFile; if(qFileName != "") { openInputFile(qFileName, qFile); } - bool success; + qFile.seekg(qline->start); + inFASTA.seekg(line->start); + + for(int i=0;inum;i++){ + + if (m->control_pressed) { + inFASTA.close(); + outFASTA.close(); + scrapFASTA.close(); + if (oligoFile != "") { outGroups.close(); } + if(qFileName != "") { qFile.close(); } + for(int i=0;iclose(); + delete fastaFileNames[i]; + } + for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } + return 0; + } + + bool success = 1; - while(!inFASTA.eof()){ Sequence currSeq(inFASTA); + string origSeq = currSeq.getUnaligned(); if (origSeq != "") { int group; @@ -185,21 +367,22 @@ int TrimSeqsCommand::execute(){ } if(!success) { trashCode += 'q'; } } - + if(barcodes.size() != 0){ success = stripBarcode(currSeq, group); if(!success){ trashCode += 'b'; } } - + if(numFPrimers != 0){ success = stripForward(currSeq); if(!success){ trashCode += 'f'; } } - + if(numRPrimers != 0){ success = stripReverse(currSeq); if(!success){ trashCode += 'r'; } } + if(minLength > 0 || maxLength > 0){ success = cullLength(currSeq); if(!success){ trashCode += 'l'; } @@ -234,10 +417,11 @@ int TrimSeqsCommand::execute(){ } gobble(inFASTA); } + inFASTA.close(); outFASTA.close(); scrapFASTA.close(); - outGroups.close(); + if (oligoFile != "") { outGroups.close(); } if(qFileName != "") { qFile.close(); } for(int i=0;ierrorOut(e, "TrimSeqsCommand", "driverCreateTrim"); + exit(1); + } +} +/**************************************************************************************************/ +int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName, string trimFile, string scrapFile, string groupFile, vector fastaNames) { + try { +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) + int process = 0; + int exitCommand = 1; + processIDS.clear(); + + //loop through and create all the processes you want + while (process != processors) { + int pid = fork(); - while(!inFASTA.eof()){ - if(inFASTA.get() == '>'){ - inFASTA >> seqName; - outGroups << seqName << '\t' << groupVector[i] << endl; - } - while (!inFASTA.eof()) { char c = inFASTA.get(); if (c == 10 || c == 13){ break; } } - } - outGroups.close(); - inFASTA.close(); + if (pid > 0) { + processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later + process++; + }else if (pid == 0){ + driverCreateTrim(filename, qFileName, (trimFile + toString(getpid()) + ".temp"), (scrapFile + toString(getpid()) + ".temp"), (groupFile + toString(getpid()) + ".temp"), fastaNames, lines[process], qLines[process]); + exit(0); + }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); } } + //force parent to wait until all the processes are done + for (int i=0;ierrorOut(e, "TrimSeqsCommand", "createProcessesCreateTrim"); exit(1); } } +/**************************************************************************************************/ + +int TrimSeqsCommand::setLines(string filename, vector& lines) { + try { + + lines.clear(); + + vector positions; + + ifstream inFASTA; + openInputFile(filename, inFASTA); + + string input; + while(!inFASTA.eof()){ + input = getline(inFASTA); + + if (input.length() != 0) { + if(input[0] == '>'){ long int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1); } + } + } + inFASTA.close(); + + int numFastaSeqs = positions.size(); + + FILE * pFile; + long size; + + //get num bytes in file + pFile = fopen (filename.c_str(),"rb"); + if (pFile==NULL) perror ("Error opening file"); + else{ + fseek (pFile, 0, SEEK_END); + size=ftell (pFile); + fclose (pFile); + } + + int numSeqsPerProcessor = numFastaSeqs / processors; + + for (int i = 0; i < processors; i++) { + long int startPos = positions[ i * numSeqsPerProcessor ]; + if(i == processors - 1){ + numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor; + }else{ + long int myEnd = positions[ (i+1) * numSeqsPerProcessor ]; + } + lines.push_back(new linePair(startPos, numSeqsPerProcessor)); + } + + return numFastaSeqs; + } + catch(exception& e) { + m->errorOut(e, "TrimSeqsCommand", "setLines"); + exit(1); + } +} //*************************************************************************************************************** -void TrimSeqsCommand::getOligos(vector& outFASTAVec){ +void TrimSeqsCommand::getOligos(vector& outFASTAVec){ //vector& outFASTAVec try { ifstream inOligos; openInputFile(oligoFile, inOligos); @@ -311,7 +567,9 @@ void TrimSeqsCommand::getOligos(vector& outFASTAVec){ groupVector.push_back(group); if(allFiles){ - outFASTAVec.push_back(new ofstream((getRootName(fastaFile) + group + ".fasta").c_str(), ios::ate)); + //outFASTAVec.push_back(new ofstream((outputDir + getRootName(getSimpleName(fastaFile)) + group + ".fasta").c_str(), ios::ate)); + outputNames.push_back((outputDir + getRootName(getSimpleName(fastaFile)) + group + ".fasta")); + outFASTAVec.push_back((outputDir + getRootName(getSimpleName(fastaFile)) + group + ".fasta")); } } } @@ -324,12 +582,10 @@ void TrimSeqsCommand::getOligos(vector& outFASTAVec){ } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "getOligos"); + m->errorOut(e, "TrimSeqsCommand", "getOligos"); exit(1); } - } - //*************************************************************************************************************** bool TrimSeqsCommand::stripBarcode(Sequence& seq, int& group){ @@ -337,6 +593,7 @@ bool TrimSeqsCommand::stripBarcode(Sequence& seq, int& group){ string rawSequence = seq.getUnaligned(); bool success = 0; //guilty until proven innocent + //can you find the barcode for(map::iterator it=barcodes.begin();it!=barcodes.end();it++){ string oligo = it->first; if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length @@ -351,11 +608,44 @@ bool TrimSeqsCommand::stripBarcode(Sequence& seq, int& group){ break; } } + + //if you found the barcode or if you don't want to allow for diffs + if ((diffs == 0) || (success == 1)) { return success; } + + else { //try aligning and see if you can find it + //can you find the barcode + for(map::iterator it=barcodes.begin();it!=barcodes.end();it++){ + string oligo = it->first; + if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length + success = 0; + break; + } + + //use needleman to align first barcode.length()+numdiffs of sequence to each barcode + Alignment* alignment = new NeedlemanOverlap(-2.0, 1.0, -1.0, (oligo.length()+diffs+1)); + Sequence* templateSeq = new Sequence("temp", rawSequence.substr(0,(oligo.length()+diffs))); + Sequence* candidateSeq = new Sequence("temp2", oligo); + Nast nast(alignment, candidateSeq, templateSeq); + + oligo = candidateSeq->getAligned(); + cout << "barcode = " << oligo << " raw = " << rawSequence.substr(0,(oligo.length())) << endl; + delete alignment; + delete templateSeq; + delete candidateSeq; + + if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){ + group = it->second; + seq.setUnaligned(rawSequence.substr(0,oligo.length())); + success = 1; + break; + } + } + } return success; } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "stripBarcode"); + m->errorOut(e, "TrimSeqsCommand", "stripBarcode"); exit(1); } @@ -383,11 +673,43 @@ bool TrimSeqsCommand::stripForward(Sequence& seq){ } } + //if you found the primer or if you don't want to allow for diffs + if ((diffs == 0) || (success == 1)) { return success; } + + else { //try aligning and see if you can find it + //can you find the primer + for(int i=0;igetAligned(); + + delete alignment; + delete templateSeq; + delete candidateSeq; + + if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){ + seq.setUnaligned(rawSequence.substr(0,oligo.length())); + success = 1; + break; + } + } + } + return success; } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "stripForward"); + m->errorOut(e, "TrimSeqsCommand", "stripForward"); exit(1); } } @@ -417,7 +739,7 @@ bool TrimSeqsCommand::stripReverse(Sequence& seq){ } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "stripReverse"); + m->errorOut(e, "TrimSeqsCommand", "stripReverse"); exit(1); } } @@ -438,7 +760,7 @@ bool TrimSeqsCommand::cullLength(Sequence& seq){ } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "cullLength"); + m->errorOut(e, "TrimSeqsCommand", "cullLength"); exit(1); } @@ -457,7 +779,7 @@ bool TrimSeqsCommand::cullHomoP(Sequence& seq){ return success; } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "cullHomoP"); + m->errorOut(e, "TrimSeqsCommand", "cullHomoP"); exit(1); } @@ -476,7 +798,7 @@ bool TrimSeqsCommand::cullAmbigs(Sequence& seq){ return success; } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "cullAmbigs"); + m->errorOut(e, "TrimSeqsCommand", "cullAmbigs"); exit(1); } @@ -492,7 +814,7 @@ bool TrimSeqsCommand::compareDNASeq(string oligo, string seq){ for(int i=0;ierrorOut(e, "TrimSeqsCommand", "compareDNASeq"); exit(1); } @@ -526,11 +848,15 @@ bool TrimSeqsCommand::compareDNASeq(string oligo, string seq){ bool TrimSeqsCommand::stripQualThreshold(Sequence& seq, ifstream& qFile){ try { string rawSequence = seq.getUnaligned(); - int seqLength = rawSequence.length(); - string name; + int seqLength; // = rawSequence.length(); + string name, temp, temp2; - qFile >> name; - if (name.length() != 0) { if(name.substr(1) != seq.getName()) { mothurOut("sequence name mismatch btwn fasta and qual file"); mothurOutEndLine(); } } + qFile >> name >> temp; + + splitAtEquals(temp2, temp); //separates length=242, temp=length, temp2=242 + convert(temp, seqLength); //converts string to int + + if (name.length() != 0) { if(name.substr(1) != seq.getName()) { m->mothurOut("sequence name mismatch btwn fasta and qual file"); m->mothurOutEndLine(); } } while (!qFile.eof()) { char c = qFile.get(); if (c == 10 || c == 13){ break; } } int score; @@ -553,7 +879,7 @@ bool TrimSeqsCommand::stripQualThreshold(Sequence& seq, ifstream& qFile){ return 1; } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "stripQualThreshold"); + m->errorOut(e, "TrimSeqsCommand", "stripQualThreshold"); exit(1); } } @@ -568,7 +894,7 @@ bool TrimSeqsCommand::cullQualAverage(Sequence& seq, ifstream& qFile){ string name; qFile >> name; - if (name[0] == '>') { if(name.substr(1) != seq.getName()) { mothurOut("sequence name mismatch btwn fasta: " + seq.getName() + " and qual file: " + name); mothurOutEndLine(); } } + if (name[0] == '>') { if(name.substr(1) != seq.getName()) { m->mothurOut("sequence name mismatch btwn fasta: " + seq.getName() + " and qual file: " + name); m->mothurOutEndLine(); } } while (!qFile.eof()) { char c = qFile.get(); if (c == 10 || c == 13){ break; } } @@ -587,7 +913,7 @@ bool TrimSeqsCommand::cullQualAverage(Sequence& seq, ifstream& qFile){ return success; } catch(exception& e) { - errorOut(e, "TrimSeqsCommand", "cullQualAverage"); + m->errorOut(e, "TrimSeqsCommand", "cullQualAverage"); exit(1); } }