X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=prcseqscommand.cpp;h=722aca537d4b88ddc38c5e1442d348c5ba11dc03;hp=d31d687c0af871ed4d5bbbde157c671900aa67d2;hb=fefd5ee1517abd3bc38b469cb2dffc85a1571c7e;hpb=1a968f34ae2d2680eaf189a197d1a21b8dfd6c03 diff --git a/prcseqscommand.cpp b/prcseqscommand.cpp index d31d687..722aca5 100644 --- a/prcseqscommand.cpp +++ b/prcseqscommand.cpp @@ -312,7 +312,7 @@ int PcrSeqsCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } int start = time(NULL); - fileAligned = true; + fileAligned = true; pairedOligos = false; string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(fastafile); } @@ -326,7 +326,7 @@ int PcrSeqsCommand::execute(){ length = 0; - if(oligosfile != ""){ readOligos(); if (m->debug) { m->mothurOut("[DEBUG]: read oligos file. numprimers = " + toString(primers.size()) + ", revprimers = " + toString(revPrimer.size()) + ".\n"); } } if (m->control_pressed) { return 0; } + if(oligosfile != ""){ readOligos(); if (m->debug) { m->mothurOut("[DEBUG]: read oligos file. numprimers = " + toString(numFPrimers) + ", revprimers = " + toString(numRPrimers) + ".\n"); } } if (m->control_pressed) { return 0; } if(ecolifile != "") { readEcoli(); } if (m->control_pressed) { return 0; } vector positions; @@ -436,20 +436,20 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string //loop through and create all the processes you want while (process != processors) { - int pid = fork(); + pid_t pid = fork(); 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){ - string locationsFile = toString(getpid()) + ".temp"; - num = driverPcr(filename, goodFileName + toString(getpid()) + ".temp", badFileName + toString(getpid()) + ".temp", locationsFile, badSeqNames, lines[process], pstart, pend, adjustNeeded); + string locationsFile = m->mothurGetpid(process) + ".temp"; + num = driverPcr(filename, goodFileName + m->mothurGetpid(process) + ".temp", badFileName + m->mothurGetpid(process) + ".temp", locationsFile, badSeqNames, lines[process], pstart, adjustNeeded); //pass numSeqs to parent ofstream out; - string tempFile = filename + toString(getpid()) + ".num.temp"; + string tempFile = filename + m->mothurGetpid(process) + ".num.temp"; m->openOutputFile(tempFile, out); - out << pstart << '\t' << pend << '\t' << adjustNeeded << endl; + out << pstart << '\t' << adjustNeeded << endl; out << num << '\t' << badSeqNames.size() << endl; for (set::iterator it = badSeqNames.begin(); it != badSeqNames.end(); it++) { out << (*it) << endl; @@ -464,8 +464,8 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string } } - string locationsFile = toString(getpid()) + ".temp"; - num = driverPcr(filename, goodFileName, badFileName, locationsFile, badSeqNames, lines[0], pstart, pend, adjustNeeded); + string locationsFile = m->mothurGetpid(process) + ".temp"; + num = driverPcr(filename, goodFileName, badFileName, locationsFile, badSeqNames, lines[0], pstart, adjustNeeded); //force parent to wait until all the processes are done for (int i=0;iopenInputFile(tempFile, in); int numBadNames = 0; string name = ""; - int tpstart = -1; int tpend = -1; bool tempAdjust = false; + int tpstart = -1; bool tempAdjust = false; if (!in.eof()) { - in >> tpstart >> tpend >> tempAdjust; m->gobble(in); + in >> tpstart >> tempAdjust; m->gobble(in); if (tempAdjust) { adjustNeeded = true; } if (tpstart != -1) { if (tpstart != pstart) { adjustNeeded = true; } if (tpstart < pstart) { pstart = tpstart; } //smallest start } - if (tpend != -1) { - if (tpend != pend) { adjustNeeded = true; } - if (tpend > pend) { pend = tpend; } //largest end - } int tempNum = 0; in >> tempNum >> numBadNames; num += tempNum; m->gobble(in); } for (int j = 0; j < numBadNames; j++) { @@ -533,7 +529,7 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string if (i!=0) {extension += toString(i) + ".temp"; processIDS.push_back(i); } // Allocate memory for thread data. - pcrData* tempPcr = new pcrData(filename, goodFileName+extension, badFileName+extension, locationsFile+extension, m, oligosfile, ecolifile, primers, revPrimer, nomatch, keepprimer, keepdots, start, end, length, pdiffs, lines[i].start, lines[i].end); + pcrData* tempPcr = new pcrData(filename, goodFileName+extension, badFileName+extension, locationsFile+extension, m, oligosfile, ecolifile, nomatch, keepprimer, keepdots, start, end, length, pdiffs, lines[i].start, lines[i].end); pDataArray.push_back(tempPcr); //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier @@ -541,7 +537,7 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string } //do your part - num = driverPcr(filename, (goodFileName+toString(processors-1)+".temp"), (badFileName+toString(processors-1)+".temp"), (locationsFile+toString(processors-1)+".temp"), badSeqNames, lines[processors-1], pstart, pend, adjustNeeded); + num = driverPcr(filename, (goodFileName+toString(processors-1)+".temp"), (badFileName+toString(processors-1)+".temp"), (locationsFile+toString(processors-1)+".temp"), badSeqNames, lines[processors-1], pstart, adjustNeeded); processIDS.push_back(processors-1); //Wait until all threads have terminated. @@ -558,11 +554,7 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string if (pDataArray[i]->pstart != pstart) { adjustNeeded = true; } if (pDataArray[i]->pstart < pstart) { pstart = pDataArray[i]->pstart; } } //smallest start - if (pDataArray[i]->pend != -1) { - if (pDataArray[i]->pend != pend) { adjustNeeded = true; } - if (pDataArray[i]->pend > pend) { pend = pDataArray[i]->pend; } - } //largest end - + for (set::iterator it = pDataArray[i]->badSeqNames.begin(); it != pDataArray[i]->badSeqNames.end(); it++) { badSeqNames.insert(*it); } CloseHandle(hThreadArray[i]); delete pDataArray[i]; @@ -580,7 +572,43 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string } #endif - if (fileAligned && adjustNeeded) { adjustDots(goodFileName, locationsFile, pstart, pend); } + + + + if (fileAligned && adjustNeeded) { + //find pend - pend is the biggest ending value, but we must account for when we adjust the start. That adjustment may make the "new" end larger then the largest end. So lets find out what that "new" end will be. + ifstream inLocations; + m->openInputFile(locationsFile, inLocations); + + while(!inLocations.eof()) { + + if (m->control_pressed) { break; } + + string name = ""; + int thisStart = -1; int thisEnd = -1; + if (numFPrimers != 0) { inLocations >> name >> thisStart; m->gobble(inLocations); } + if (numRPrimers != 0) { inLocations >> name >> thisEnd; m->gobble(inLocations); } + else { pend = -1; break; } + + int myDiff = 0; + if (pstart != -1) { + if (thisStart != -1) { + if (thisStart != pstart) { myDiff += (thisStart - pstart); } + } + } + + int myEnd = thisEnd + myDiff; + //cout << name << '\t' << thisStart << '\t' << thisEnd << " diff = " << myDiff << '\t' << myEnd << endl; + + if (thisEnd != -1) { + if (myEnd > pend) { pend = myEnd; } + } + + } + inLocations.close(); + + adjustDots(goodFileName, locationsFile, pstart, pend); + }else { m->mothurRemove(locationsFile); } return num; @@ -592,7 +620,7 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string } //********************************************************************************************************************** -int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta, string locationsName, set& badSeqNames, linePair filePos, int& pstart, int& pend, bool& adjustNeeded){ +int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta, string locationsName, set& badSeqNames, linePair filePos, int& pstart, bool& adjustNeeded){ try { ofstream goodFile; m->openOutputFile(goodFasta, goodFile); @@ -611,12 +639,25 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta bool done = false; int count = 0; set lengths; - vector< set > locations; //locations[0] = beginning locations, locations[1] = ending locations - locations.resize(2); + set locations; //locations[0] = beginning locations, //pdiffs, bdiffs, primers, barcodes, revPrimers - map faked; - TrimOligos trim(pdiffs, 0, primers, faked, revPrimer); + map primers; + map barcodes; //not used + vector revPrimer; + if (pairedOligos) { + map primerPairs = oligos.getPairedPrimers(); + for (map::iterator it = primerPairs.begin(); it != primerPairs.end(); it++) { + primers[(it->second).forward] = it->first; + revPrimer.push_back((it->second).reverse); + } + if (pdiffs != 0) { m->mothurOut("[WARNING]: Pcr.seqs is only designed to allow diffs in forward primers. Reverse primers must be an exact match.\n"); } + }else{ + primers = oligos.getPrimers(); + revPrimer = oligos.getReversePrimers(); + } + + TrimOligos trim(pdiffs, 0, primers, barcodes, revPrimer); while (!done) { @@ -759,9 +800,8 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta if(goodSeq == 1) { currSeq.printSequence(goodFile); if (m->debug) { m->mothurOut("[DEBUG]: " + locationsString + "\n"); } + if (thisPStart != -1) { locations.insert(thisPStart); } if (locationsString != "") { locationsFile << locationsString; } - if (thisPStart != -1) { locations[0].insert(thisPStart); } - if (thisPEnd != -1) { locations[1].insert(thisPEnd); } } else { badSeqNames.insert(currSeq.getName()); @@ -792,9 +832,8 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta if (m->debug) { m->mothurOut("[DEBUG]: fileAligned = " + toString(fileAligned) +'\n'); } if (fileAligned && !keepdots) { //print out smallest start value and largest end value - if ((locations[0].size() > 1) || (locations[1].size() > 1)) { adjustNeeded = true; } - if (primers.size() != 0) { set::iterator it = locations[0].begin(); pstart = *it; } - if (revPrimer.size() != 0) { set::reverse_iterator it2 = locations[1].rbegin(); pend = *it2; } + if (locations.size() > 1) { adjustNeeded = true; } + if (primers.size() != 0) { set::iterator it = locations.begin(); pstart = *it; } } return count; @@ -837,6 +876,7 @@ int PcrSeqsCommand::adjustDots(string goodFasta, string locations, int pstart, i set lengths; //cout << pstart << '\t' << pend << endl; + //if (pstart > pend) { //swap them while(!inFasta.eof()) { if(m->control_pressed) { break; } @@ -845,8 +885,10 @@ int PcrSeqsCommand::adjustDots(string goodFasta, string locations, int pstart, i string name = ""; int thisStart = -1; int thisEnd = -1; - if (primers.size() != 0) { inLocations >> name >> thisStart; m->gobble(inLocations); } - if (revPrimer.size() != 0) { inLocations >> name >> thisEnd; m->gobble(inLocations); } + if (numFPrimers != 0) { inLocations >> name >> thisStart; m->gobble(inLocations); } + if (numRPrimers != 0) { inLocations >> name >> thisEnd; m->gobble(inLocations); } + + //cout << seq.getName() << '\t' << thisStart << '\t' << thisEnd << '\t' << seq.getAligned().length() << endl; //cout << seq.getName() << '\t' << pstart << '\t' << pend << endl; @@ -887,7 +929,8 @@ int PcrSeqsCommand::adjustDots(string goodFasta, string locations, int pstart, i //cout << "final lengths = \n"; //for (set::iterator it = lengths.begin(); it != lengths.end(); it++) { - // cout << *it << endl; + //cout << *it << endl; + // cout << lengths.count(*it) << endl; // } return 0; @@ -897,109 +940,6 @@ int PcrSeqsCommand::adjustDots(string goodFasta, string locations, int pstart, i exit(1); } } -//********************************************************************/ -string PcrSeqsCommand::reverseOligo(string oligo){ - try { - string reverse = ""; - - for(int i=oligo.length()-1;i>=0;i--){ - - if(oligo[i] == 'A') { reverse += 'T'; } - else if(oligo[i] == 'T'){ reverse += 'A'; } - else if(oligo[i] == 'U'){ reverse += 'A'; } - - else if(oligo[i] == 'G'){ reverse += 'C'; } - else if(oligo[i] == 'C'){ reverse += 'G'; } - - else if(oligo[i] == 'R'){ reverse += 'Y'; } - else if(oligo[i] == 'Y'){ reverse += 'R'; } - - else if(oligo[i] == 'M'){ reverse += 'K'; } - else if(oligo[i] == 'K'){ reverse += 'M'; } - - else if(oligo[i] == 'W'){ reverse += 'W'; } - else if(oligo[i] == 'S'){ reverse += 'S'; } - - else if(oligo[i] == 'B'){ reverse += 'V'; } - else if(oligo[i] == 'V'){ reverse += 'B'; } - - else if(oligo[i] == 'D'){ reverse += 'H'; } - else if(oligo[i] == 'H'){ reverse += 'D'; } - - else { reverse += 'N'; } - } - - - return reverse; - } - catch(exception& e) { - m->errorOut(e, "PcrSeqsCommand", "reverseOligo"); - exit(1); - } -} - -//*************************************************************************************************************** -bool PcrSeqsCommand::readOligos(){ - try { - ifstream inOligos; - m->openInputFile(oligosfile, inOligos); - - string type, oligo, group; - int primerCount = 0; - - while(!inOligos.eof()){ - - inOligos >> type; - - if(type[0] == '#'){ //ignore - while (!inOligos.eof()) { char c = inOligos.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there - m->gobble(inOligos); - }else{ - m->gobble(inOligos); - //make type case insensitive - for(int i=0;i> oligo; - - for(int i=0;i> group; - }else if((type == "LINKER")||(type == "SPACER")) {;} - else{ m->mothurOut(type + " is not recognized as a valid type. Choices are forward, reverse, linker, spacer and barcode. Ignoring " + oligo + "."); m->mothurOutEndLine(); m->control_pressed = true; } - } - m->gobble(inOligos); - } - inOligos.close(); - - if ((primers.size() == 0) && (revPrimer.size() == 0)) { - m->mothurOut("[ERROR]: your oligos file does not contain valid primers or reverse primers. Please correct."); m->mothurOutEndLine(); - m->control_pressed = true; - return false; - } - - return true; - - }catch(exception& e) { - m->errorOut(e, "PcrSeqsCommand", "readOligos"); - exit(1); - } -} //*************************************************************************************************************** bool PcrSeqsCommand::readEcoli(){ try { @@ -1254,7 +1194,7 @@ int PcrSeqsCommand::readCount(set badSeqNames){ //check for groups that have been eliminated CountTable ct; if (ct.testGroups(goodCountFile)) { - ct.readTable(goodCountFile, true); + ct.readTable(goodCountFile, true, false); ct.printTable(goodCountFile); } @@ -1271,6 +1211,35 @@ int PcrSeqsCommand::readCount(set badSeqNames){ exit(1); } } +//*************************************************************************************************************** + +int PcrSeqsCommand::readOligos(){ + try { + oligos.read(oligosfile); + + if (m->control_pressed) { return false; } //error in reading oligos + + if (oligos.hasPairedBarcodes()) { + pairedOligos = true; + numFPrimers = oligos.getPairedPrimers().size(); + }else { + pairedOligos = false; + numFPrimers = oligos.getPrimers().size(); + } + numRPrimers = oligos.getReversePrimers().size(); + + if (oligos.getLinkers().size() != 0) { m->mothurOut("[WARNING]: pcr.seqs is not setup to remove linkers, ignoring.\n"); } + if (oligos.getSpacers().size() != 0) { m->mothurOut("[WARNING]: pcr.seqs is not setup to remove spacers, ignoring.\n"); } + + return true; + + } + catch(exception& e) { + m->errorOut(e, "PcrSeqsCommand", "readOligos"); + exit(1); + } +} + /**************************************************************************************/