X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=prcseqscommand.cpp;h=722aca537d4b88ddc38c5e1442d348c5ba11dc03;hp=cf93cb796104aa90428e5466546c9839789216f7;hb=fefd5ee1517abd3bc38b469cb2dffc85a1571c7e;hpb=b0997605981902442138b9309e9c43d95c3ba10a diff --git a/prcseqscommand.cpp b/prcseqscommand.cpp index cf93cb7..722aca5 100644 --- a/prcseqscommand.cpp +++ b/prcseqscommand.cpp @@ -312,6 +312,7 @@ int PcrSeqsCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } int start = time(NULL); + fileAligned = true; pairedOligos = false; string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(fastafile); } @@ -325,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; @@ -352,8 +353,7 @@ int PcrSeqsCommand::execute(){ if (m->control_pressed) { return 0; } set badNames; - if(processors == 1) { numFastaSeqs = driverPcr(fastafile, trimSeqFile, badSeqFile, badNames, lines[0]); } - else { numFastaSeqs = createProcesses(fastafile, trimSeqFile, badSeqFile, badNames); } + numFastaSeqs = createProcesses(fastafile, trimSeqFile, badSeqFile, badNames); if (m->control_pressed) { return 0; } @@ -429,23 +429,27 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string vector processIDS; int process = 1; int num = 0; + int pstart = -1; int pend = -1; + bool adjustNeeded = false; #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) //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){ - num = driverPcr(filename, goodFileName + toString(getpid()) + ".temp", badFileName + toString(getpid()) + ".temp", badSeqNames, lines[process]); + 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' << adjustNeeded << endl; out << num << '\t' << badSeqNames.size() << endl; for (set::iterator it = badSeqNames.begin(); it != badSeqNames.end(); it++) { out << (*it) << endl; @@ -460,7 +464,8 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string } } - num = driverPcr(filename, goodFileName, badFileName, badSeqNames, lines[0]); + 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 = ""; - if (!in.eof()) { int tempNum = 0; in >> tempNum >> numBadNames; num += tempNum; m->gobble(in); } + int tpstart = -1; bool tempAdjust = false; + + if (!in.eof()) { + 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 + } + int tempNum = 0; in >> tempNum >> numBadNames; num += tempNum; m->gobble(in); + } for (int j = 0; j < numBadNames; j++) { in >> name; m->gobble(in); badSeqNames.insert(name); @@ -485,6 +501,9 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string m->appendFiles((badFileName + toString(processIDS[i]) + ".temp"), badFileName); m->mothurRemove((badFileName + toString(processIDS[i]) + ".temp")); + + m->appendFiles((toString(processIDS[i]) + ".temp"), locationsFile); + m->mothurRemove((toString(processIDS[i]) + ".temp")); } #else @@ -498,6 +517,11 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string DWORD dwThreadIdArray[processors-1]; HANDLE hThreadArray[processors-1]; + string locationsFile = "locationsFile.txt"; + m->mothurRemove(locationsFile); + m->mothurRemove(goodFileName); + m->mothurRemove(badFileName); + //Create processor worker threads. for( int i=0; icount != pDataArray[i]->fend) { m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->fend) + " sequences assigned to it, quitting. \n"); m->control_pressed = true; } + if (pDataArray[i]->adjustNeeded) { adjustNeeded = true; } + if (pDataArray[i]->pstart != -1) { + if (pDataArray[i]->pstart != pstart) { adjustNeeded = true; } + if (pDataArray[i]->pstart < pstart) { pstart = pDataArray[i]->pstart; } + } //smallest start + for (set::iterator it = pDataArray[i]->badSeqNames.begin(); it != pDataArray[i]->badSeqNames.end(); it++) { badSeqNames.insert(*it); } CloseHandle(hThreadArray[i]); delete pDataArray[i]; @@ -536,10 +566,50 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string m->appendFiles((badFileName + toString(processIDS[i]) + ".temp"), badFileName); m->mothurRemove((badFileName + toString(processIDS[i]) + ".temp")); + + m->appendFiles((locationsFile+toString(processIDS[i]) + ".temp"), locationsFile); + m->mothurRemove((locationsFile+toString(processIDS[i]) + ".temp")); } #endif + + + 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; } @@ -550,13 +620,16 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string } //********************************************************************************************************************** -int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta, set& badSeqNames, linePair filePos){ +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); ofstream badFile; m->openOutputFile(badFasta, badFile); + + ofstream locationsFile; + m->openOutputFile(locationsName, locationsFile); ifstream inFASTA; m->openInputFile(filename, inFASTA); @@ -566,10 +639,25 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta bool done = false; int count = 0; set lengths; + 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) { @@ -577,7 +665,15 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta Sequence currSeq(inFASTA); m->gobble(inFASTA); + if (fileAligned) { //assume aligned until proven otherwise + lengths.insert(currSeq.getAligned().length()); + if (lengths.size() > 1) { fileAligned = false; } + } + string trashCode = ""; + string locationsString = ""; + int thisPStart = -1; + int thisPEnd = -1; if (currSeq.getName() != "") { if (m->debug) { m->mothurOut("[DEBUG]: seq name = " + currSeq.getName() + ".\n"); } @@ -596,13 +692,25 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta else{ //are you aligned if (aligned) { - if (!keepprimer) { - if (keepdots) { currSeq.filterToPos(mapAligned[primerEnd]); } - else { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerEnd])); } + if (!keepprimer) { + if (keepdots) { currSeq.filterToPos(mapAligned[primerEnd-1]+1); } //mapAligned[primerEnd-1] is the location of the last base in the primer. we want to trim to the space just after that. The -1 & +1 ensures if the primer is followed by gaps they are not trimmed causing an aligned sequence dataset to become unaligned. + else { + currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerEnd-1]+1)); + if (fileAligned) { + thisPStart = mapAligned[primerEnd-1]+1; //locations[0].insert(mapAligned[primerEnd-1]+1); + locationsString += currSeq.getName() + "\t" + toString(mapAligned[primerEnd-1]+1) + "\n"; + } + } } else { if (keepdots) { currSeq.filterToPos(mapAligned[primerStart]); } - else { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerStart])); } + else { + currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerStart])); + if (fileAligned) { + thisPStart = mapAligned[primerStart]; //locations[0].insert(mapAligned[primerStart]); + locationsString += currSeq.getName() + "\t" + toString(mapAligned[primerStart]) + "\n"; + } + } } isAligned(currSeq.getAligned(), mapAligned); }else { @@ -617,16 +725,28 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta int primerStart = 0; int primerEnd = 0; bool good = trim.findReverse(currSeq, primerStart, primerEnd); if(!good){ if (nomatch == "reject") { goodSeq = false; } trashCode += "r"; } - else{ - //are you aligned + else{ + //are you aligned if (aligned) { if (!keepprimer) { if (keepdots) { currSeq.filterFromPos(mapAligned[primerStart]); } - else { currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerStart])); } + else { + currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerStart])); + if (fileAligned) { + thisPEnd = mapAligned[primerStart]; //locations[1].insert(mapAligned[primerStart]); + locationsString += currSeq.getName() + "\t" + toString(mapAligned[primerStart]) + "\n"; + } + } } else { - if (keepdots) { currSeq.filterFromPos(mapAligned[primerEnd]); } - else { currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerEnd])); } + if (keepdots) { currSeq.filterFromPos(mapAligned[primerEnd-1]+1); } + else { + currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerEnd-1]+1)); + if (fileAligned) { + thisPEnd = mapAligned[primerEnd-1]+1; //locations[1].insert(mapAligned[primerEnd-1]+1); + locationsString += currSeq.getName() + "\t" + toString(mapAligned[primerEnd-1]+1) + "\n"; + } + } } } else { @@ -637,8 +757,7 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta } }else if (ecolifile != "") { //make sure the seqs are aligned - lengths.insert(currSeq.getAligned().length()); - if (lengths.size() > 1) { m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); m->control_pressed = true; break; } + if (!fileAligned) { m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); m->control_pressed = true; break; } else if (currSeq.getAligned().length() != length) { m->mothurOut("[ERROR]: seqs are not the same length as ecoli seq. When using ecoli option your sequences must be aligned and the same length as the ecoli sequence.\n"); m->control_pressed = true; break; }else { @@ -653,8 +772,7 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta } }else{ //using start and end to trim //make sure the seqs are aligned - lengths.insert(currSeq.getAligned().length()); - if (lengths.size() > 1) { m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); m->control_pressed = true; break; } + if (!fileAligned) { m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); m->control_pressed = true; break; } else { if (end != -1) { if (end > currSeq.getAligned().length()) { m->mothurOut("[ERROR]: end is longer than your sequence length, aborting.\n"); m->control_pressed = true; break; } @@ -679,7 +797,12 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta //trimming removed all bases if (currSeq.getUnaligned() == "") { goodSeq = false; } - if(goodSeq == 1) { currSeq.printSequence(goodFile); } + if(goodSeq == 1) { + currSeq.printSequence(goodFile); + if (m->debug) { m->mothurOut("[DEBUG]: " + locationsString + "\n"); } + if (thisPStart != -1) { locations.insert(thisPStart); } + if (locationsString != "") { locationsFile << locationsString; } + } else { badSeqNames.insert(currSeq.getName()); currSeq.setName(currSeq.getName() + '|' + trashCode); @@ -696,15 +819,23 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta #endif //report progress - if((count) % 100 == 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); } + if((count) % 100 == 0){ m->mothurOutJustToScreen("Processing sequence: " + toString(count)+"\n"); } } //report progress - if((count) % 100 != 0){ m->mothurOut("Processing sequence: " + toString(count)); m->mothurOutEndLine(); } + if((count) % 100 != 0){ m->mothurOutJustToScreen("Processing sequence: " + toString(count)+"\n"); } badFile.close(); goodFile.close(); inFASTA.close(); - + locationsFile.close(); + + if (m->debug) { m->mothurOut("[DEBUG]: fileAligned = " + toString(fileAligned) +'\n'); } + + if (fileAligned && !keepdots) { //print out smallest start value and largest end value + if (locations.size() > 1) { adjustNeeded = true; } + if (primers.size() != 0) { set::iterator it = locations.begin(); pstart = *it; } + } + return count; } catch(exception& e) { @@ -731,108 +862,83 @@ bool PcrSeqsCommand::isAligned(string seq, map& aligned){ 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'; } +//********************************************************************************************************************** +int PcrSeqsCommand::adjustDots(string goodFasta, string locations, int pstart, int pend){ + try { + ifstream inFasta; + m->openInputFile(goodFasta, inFasta); + + ifstream inLocations; + m->openInputFile(locations, inLocations); + + ofstream out; + m->openOutputFile(goodFasta+".temp", out); + + set lengths; + //cout << pstart << '\t' << pend << endl; + //if (pstart > pend) { //swap them + + while(!inFasta.eof()) { + if(m->control_pressed) { break; } - else if(oligo[i] == 'G'){ reverse += 'C'; } - else if(oligo[i] == 'C'){ reverse += 'G'; } + Sequence seq(inFasta); m->gobble(inFasta); - else if(oligo[i] == 'R'){ reverse += 'Y'; } - else if(oligo[i] == 'Y'){ reverse += 'R'; } + 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 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'; } + //cout << seq.getName() << '\t' << thisStart << '\t' << thisEnd << '\t' << seq.getAligned().length() << endl; + //cout << seq.getName() << '\t' << pstart << '\t' << pend << endl; - else if(oligo[i] == 'B'){ reverse += 'V'; } - else if(oligo[i] == 'V'){ reverse += 'B'; } + if (name != seq.getName()) { m->mothurOut("[ERROR]: name mismatch in pcr.seqs.\n"); } + else { + if (pstart != -1) { + if (thisStart != -1) { + if (thisStart != pstart) { + string dots = ""; + for (int i = pstart; i < thisStart; i++) { dots += "."; } + thisEnd += dots.length(); + dots += seq.getAligned(); + seq.setAligned(dots); + } + } + } + + if (pend != -1) { + if (thisEnd != -1) { + if (thisEnd != pend) { + string dots = seq.getAligned(); + for (int i = thisEnd; i < pend; i++) { dots += "."; } + seq.setAligned(dots); + } + } + } + lengths.insert(seq.getAligned().length()); + } - else if(oligo[i] == 'D'){ reverse += 'H'; } - else if(oligo[i] == 'H'){ reverse += 'D'; } - - else { reverse += 'N'; } + seq.printSequence(out); } - - - 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; - } + inFasta.close(); + inLocations.close(); + out.close(); + m->mothurRemove(locations); + m->mothurRemove(goodFasta); + m->renameFile(goodFasta+".temp", goodFasta); - return true; + //cout << "final lengths = \n"; + //for (set::iterator it = lengths.begin(); it != lengths.end(); it++) { + //cout << *it << endl; + // cout << lengths.count(*it) << endl; + // } - }catch(exception& e) { - m->errorOut(e, "PcrSeqsCommand", "readOligos"); - exit(1); - } + return 0; + } + catch(exception& e) { + m->errorOut(e, "PcrSeqsCommand", "adjustDots"); + exit(1); + } } //*************************************************************************************************************** bool PcrSeqsCommand::readEcoli(){ @@ -1088,7 +1194,7 @@ int PcrSeqsCommand::readCount(set badSeqNames){ //check for groups that have been eliminated CountTable ct; if (ct.testGroups(goodCountFile)) { - ct.readTable(goodCountFile); + ct.readTable(goodCountFile, true, false); ct.printTable(goodCountFile); } @@ -1105,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); + } +} + /**************************************************************************************/