+++ /dev/null
-#ifndef Mothur_pcrseqscommand_h
-#define Mothur_pcrseqscommand_h
-
-//
-// pcrseqscommand.h
-// Mothur
-//
-// Created by Sarah Westcott on 3/14/12.
-// Copyright (c) 2012 Schloss Lab. All rights reserved.
-//
-
-
-#include "command.hpp"
-#include "sequence.hpp"
-#include "trimoligos.h"
-#include "alignment.hpp"
-#include "needlemanoverlap.hpp"
-
-class PcrSeqsCommand : public Command {
-public:
- PcrSeqsCommand(string);
- PcrSeqsCommand();
- ~PcrSeqsCommand(){}
-
- vector<string> setParameters();
- string getCommandName() { return "pcr.seqs"; }
- string getCommandCategory() { return "Sequence Processing"; }
- string getHelpString();
- string getCitation() { return "http://www.mothur.org/wiki/Pcr.seqs"; }
- string getDescription() { return "pcr.seqs"; }
-
- int execute();
- void help() { m->mothurOut(getHelpString()); }
-
-private:
-
- struct linePair {
- unsigned long long start;
- unsigned long long end;
- linePair(unsigned long long i, unsigned long long j) : start(i), end(j) {}
- linePair() {}
- };
-
- vector<linePair> lines;
- bool getOligos(vector<vector<string> >&, vector<vector<string> >&, vector<vector<string> >&);
- bool abort, keepprimer, keepdots;
- string fastafile, oligosfile, taxfile, groupfile, namefile, ecolifile, outputDir, nomatch;
- int start, end, processors, length;
-
- vector<string> revPrimer, outputNames;
- vector<string> primers;
-
- int writeAccnos(set<string>);
- int readName(set<string>&);
- int readGroup(set<string>);
- int readTax(set<string>);
- bool readOligos();
- bool readEcoli();
- int driverPcr(string, string, string, set<string>&, linePair);
- int createProcesses(string, string, string, set<string>&);
- bool findForward(Sequence&, int&, int&);
- bool findReverse(Sequence&, int&, int&);
- bool isAligned(string, map<int, int>&);
- bool compareDNASeq(string, string);
- string reverseOligo(string);
-};
-
-/**************************************************************************************************/
-//custom data structure for threads to use.
-// This is passed by void pointer so it can be any data type
-// that can be passed using a single void pointer (LPVOID).
-struct pcrData {
- string filename;
- string goodFasta, badFasta, oligosfile, ecolifile, nomatch;
- unsigned long long fstart;
- unsigned long long fend;
- int count, start, end, length;
- MothurOut* m;
- vector<string> primers;
- vector<string> revPrimer;
- set<string> badSeqNames;
- bool keepprimer, keepdots;
-
-
- pcrData(){}
- pcrData(string f, string gf, string bfn, MothurOut* mout, string ol, string ec, vector<string> pr, vector<string> rpr, string nm, bool kp, bool kd, int st, int en, int l, unsigned long long fst, unsigned long long fen) {
- filename = f;
- goodFasta = gf;
- badFasta = bfn;
- m = mout;
- oligosfile = ol;
- ecolifile = ec;
- primers = pr;
- revPrimer = rpr;
- nomatch = nm;
- keepprimer = kp;
- keepdots = kd;
- start = st;
- end = en;
- length = l;
- fstart = fst;
- fend = fen;
- count = 0;
- }
-};
-/**************************************************************************************************/
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
-#else
-static DWORD WINAPI MyPcrThreadFunction(LPVOID lpParam){
- pcrData* pDataArray;
- pDataArray = (pcrData*)lpParam;
-
- try {
- ofstream goodFile;
- pDataArray->m->openOutputFile(pDataArray->goodFasta, goodFile);
-
- ofstream badFile;
- pDataArray->m->openOutputFile(pDataArray->badFasta, badFile);
-
- ifstream inFASTA;
- pDataArray->m->openInputFile(pDataArray->filename, inFASTA);
-
- //print header if you are process 0
- if ((pDataArray->fstart == 0) || (pDataArray->fstart == 1)) {
- inFASTA.seekg(0);
- }else { //this accounts for the difference in line endings.
- inFASTA.seekg(pDataArray->fstart-1); pDataArray->m->gobble(inFASTA);
- }
-
- set<int> lengths;
- pDataArray->count = pDataArray->fend;
- for(int i = 0; i < pDataArray->fend; i++){ //end is the number of sequences to process
-
- if (pDataArray->m->control_pressed) { break; }
-
- Sequence currSeq(inFASTA); pDataArray->m->gobble(inFASTA);
-
- string trashCode = "";
- if (currSeq.getName() != "") {
-
- bool goodSeq = true;
- if (pDataArray->oligosfile != "") {
- map<int, int> mapAligned;
- //bool aligned = isAligned(currSeq.getAligned(), mapAligned);
- ///////////////////////////////////////////////////////////////
- bool aligned = false;
- string seq = currSeq.getAligned();
- int countBases = 0;
- for (int k = 0; k < seq.length(); k++) {
- if (!isalpha(seq[k])) { aligned = true; }
- else { mapAligned[countBases] = k; countBases++; } //maps location in unaligned -> location in aligned.
- } //ie. the 3rd base may be at spot 10 in the alignment
- //later when we trim we want to trim from spot 10.
- ///////////////////////////////////////////////////////////////
-
- //process primers
- if (pDataArray->primers.size() != 0) {
- int primerStart = 0; int primerEnd = 0;
- //bool good = findForward(currSeq, primerStart, primerEnd);
- ///////////////////////////////////////////////////////////////
- bool good = false;
- string rawSequence = currSeq.getUnaligned();
-
- for(int j=0;j<pDataArray->primers.size();j++){
- string oligo = pDataArray->primers[j];
-
- if (pDataArray->m->control_pressed) { primerStart = 0; primerEnd = 0; good = false; break; }
-
- if(rawSequence.length() < oligo.length()) { break; }
-
- //search for primer
- int olength = oligo.length();
- for (int l = 0; l < rawSequence.length()-olength; l++){
- if (pDataArray->m->control_pressed) { primerStart = 0; primerEnd = 0; good = false; break; }
- string rawChunk = rawSequence.substr(l, olength);
- //compareDNASeq(oligo, rawChunk)
- ////////////////////////////////////////////////////////
- bool success = 1;
- for(int k=0;k<olength;k++){
-
- if(oligo[k] != rawChunk[k]){
- if(oligo[k] == 'A' || oligo[k] == 'T' || oligo[k] == 'G' || oligo[k] == 'C') { success = 0; }
- else if((oligo[k] == 'N' || oligo[k] == 'I') && (rawChunk[k] == 'N')) { success = 0; }
- else if(oligo[k] == 'R' && (rawChunk[k] != 'A' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'Y' && (rawChunk[k] != 'C' && rawChunk[k] != 'T')) { success = 0; }
- else if(oligo[k] == 'M' && (rawChunk[k] != 'C' && rawChunk[k] != 'A')) { success = 0; }
- else if(oligo[k] == 'K' && (rawChunk[k] != 'T' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'W' && (rawChunk[k] != 'T' && rawChunk[k] != 'A')) { success = 0; }
- else if(oligo[k] == 'S' && (rawChunk[k] != 'C' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'B' && (rawChunk[k] != 'C' && rawChunk[k] != 'T' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'D' && (rawChunk[k] != 'A' && rawChunk[k] != 'T' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'H' && (rawChunk[k] != 'A' && rawChunk[k] != 'T' && rawChunk[k] != 'C')) { success = 0; }
- else if(oligo[k] == 'V' && (rawChunk[k] != 'A' && rawChunk[k] != 'C' && rawChunk[k] != 'G')) { success = 0; }
-
- if(success == 0) { break; }
- }
- else{
- success = 1;
- }
- }
-
- ////////////////////////////////////////////////////////////////////
- if(success) {
- primerStart = j;
- primerEnd = primerStart + olength;
- good = true; break;
- }
- }
- if (good) { break; }
- }
-
- if (!good) { primerStart = 0; primerEnd = 0; }
- ///////////////////////////////////////////////////////////////
-
-
- if(!good){ if (pDataArray->nomatch == "reject") { goodSeq = false; } trashCode += "f"; }
- else{
- //are you aligned
- if (aligned) {
- if (!pDataArray->keepprimer) {
- if (pDataArray->keepdots) { currSeq.filterToPos(mapAligned[primerEnd]); }
- else { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerEnd])); }
- }
- else {
- if (pDataArray->keepdots) { currSeq.filterToPos(mapAligned[primerStart]); }
- else { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerStart])); }
- }
- }else {
- if (!pDataArray->keepprimer) { currSeq.setAligned(currSeq.getUnaligned().substr(primerEnd)); }
- else { currSeq.setAligned(currSeq.getUnaligned().substr(primerStart)); }
- }
- }
- }
-
- //process reverse primers
- if (pDataArray->revPrimer.size() != 0) {
- int primerStart = 0; int primerEnd = 0;
- bool good = false;
- //findReverse(currSeq, primerStart, primerEnd);
- ///////////////////////////////////////////////////////////////
- string rawSequence = currSeq.getUnaligned();
-
- for(int j=0;j<pDataArray->revPrimer.size();j++){
- string oligo = pDataArray->revPrimer[j];
- if (pDataArray->m->control_pressed) { primerStart = 0; primerEnd = 0; good = false; break; }
- if(rawSequence.length() < oligo.length()) { break; }
-
- //search for primer
- int olength = oligo.length();
- for (int l = rawSequence.length()-olength; l >= 0; l--){
-
- string rawChunk = rawSequence.substr(l, olength);
- //compareDNASeq(oligo, rawChunk)
- ////////////////////////////////////////////////////////
- bool success = 1;
- for(int k=0;k<olength;k++){
-
- if(oligo[k] != rawChunk[k]){
- if(oligo[k] == 'A' || oligo[k] == 'T' || oligo[k] == 'G' || oligo[k] == 'C') { success = 0; }
- else if((oligo[k] == 'N' || oligo[k] == 'I') && (rawChunk[k] == 'N')) { success = 0; }
- else if(oligo[k] == 'R' && (rawChunk[k] != 'A' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'Y' && (rawChunk[k] != 'C' && rawChunk[k] != 'T')) { success = 0; }
- else if(oligo[k] == 'M' && (rawChunk[k] != 'C' && rawChunk[k] != 'A')) { success = 0; }
- else if(oligo[k] == 'K' && (rawChunk[k] != 'T' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'W' && (rawChunk[k] != 'T' && rawChunk[k] != 'A')) { success = 0; }
- else if(oligo[k] == 'S' && (rawChunk[k] != 'C' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'B' && (rawChunk[k] != 'C' && rawChunk[k] != 'T' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'D' && (rawChunk[k] != 'A' && rawChunk[k] != 'T' && rawChunk[k] != 'G')) { success = 0; }
- else if(oligo[k] == 'H' && (rawChunk[k] != 'A' && rawChunk[k] != 'T' && rawChunk[k] != 'C')) { success = 0; }
- else if(oligo[k] == 'V' && (rawChunk[k] != 'A' && rawChunk[k] != 'C' && rawChunk[k] != 'G')) { success = 0; }
-
- if(success == 0) { break; }
- }
- else{
- success = 1;
- }
- }
-
- ////////////////////////////////////////////////////////////////////
- if(success) {
- primerStart = j;
- primerEnd = primerStart + olength;
- good = true; break;
- }
- }
- if (good) { break; }
- }
-
- if (!good) { primerStart = 0; primerEnd = 0; }
-
- ///////////////////////////////////////////////////////////////
- if(!good){ if (pDataArray->nomatch == "reject") { goodSeq = false; } trashCode += "r"; }
- else{
- //are you aligned
- if (aligned) {
- if (!pDataArray->keepprimer) {
- if (pDataArray->keepdots) { currSeq.filterFromPos(mapAligned[primerStart]); }
- else { currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerStart])); }
- }
- else {
- if (pDataArray->keepdots) { currSeq.filterFromPos(mapAligned[primerEnd]); }
- else { currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerEnd])); }
- } }
- else {
- if (!pDataArray->keepprimer) { currSeq.setAligned(currSeq.getUnaligned().substr(0, primerStart)); }
- else { currSeq.setAligned(currSeq.getUnaligned().substr(0, primerEnd)); }
- }
- }
- }
- }else if (pDataArray->ecolifile != "") {
- //make sure the seqs are aligned
- lengths.insert(currSeq.getAligned().length());
- if (lengths.size() > 1) { pDataArray->m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); pDataArray->m->control_pressed = true; break; }
- else if (currSeq.getAligned().length() != pDataArray->length) {
- pDataArray->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"); pDataArray->m->control_pressed = true; break;
- }else {
- if (pDataArray->keepdots) {
- currSeq.filterToPos(pDataArray->start);
- currSeq.filterFromPos(pDataArray->end);
- }else {
- string seqString = currSeq.getAligned().substr(0, pDataArray->end);
- seqString = seqString.substr(pDataArray->start);
- currSeq.setAligned(seqString);
- }
- }
- }else{ //using start and end to trim
- //make sure the seqs are aligned
- lengths.insert(currSeq.getAligned().length());
- if (lengths.size() > 1) { pDataArray->m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); pDataArray->m->control_pressed = true; break; }
- else {
- if (pDataArray->end != -1) {
- if (pDataArray->end > currSeq.getAligned().length()) { pDataArray->m->mothurOut("[ERROR]: end is longer than your sequence length, aborting.\n"); pDataArray->m->control_pressed = true; break; }
- else {
- if (pDataArray->keepdots) { currSeq.filterFromPos(pDataArray->end); }
- else {
- string seqString = currSeq.getAligned().substr(0, pDataArray->end);
- currSeq.setAligned(seqString);
- }
- }
- }
- if (pDataArray->start != -1) {
- if (pDataArray->keepdots) { currSeq.filterToPos(pDataArray->start); }
- else {
- string seqString = currSeq.getAligned().substr(pDataArray->start);
- currSeq.setAligned(seqString);
- }
- }
-
- }
- }
-
- if(goodSeq == 1) { currSeq.printSequence(goodFile); }
- else {
- pDataArray->badSeqNames.insert(currSeq.getName());
- currSeq.setName(currSeq.getName() + '|' + trashCode);
- currSeq.printSequence(badFile);
- }
- }
-
- //report progress
- if((i+1) % 100 == 0){ pDataArray->m->mothurOut("Processing sequence: " + toString(i+1)); pDataArray->m->mothurOutEndLine(); }
- }
- //report progress
- if((pDataArray->count) % 100 != 0){ pDataArray->m->mothurOut("Thread Processing sequence: " + toString(pDataArray->count)); pDataArray->m->mothurOutEndLine(); }
-
- goodFile.close();
- inFASTA.close();
- badFile.close();
-
- return 0;
-
- }
- catch(exception& e) {
- pDataArray->m->errorOut(e, "PcrSeqsCommand", "MyPcrThreadFunction");
- exit(1);
- }
-}
-
-#endif
-
-/**************************************************************************************************/
-
-
-
-#endif
projects / mothur.git / blobdiff